diff --git a/code/obs_operators/geocape_ch4_mod.f b/code/obs_operators/geocape_ch4_mod.f index d2a19ad..6f3d05b 100644 --- a/code/obs_operators/geocape_ch4_mod.f +++ b/code/obs_operators/geocape_ch4_mod.f @@ -1012,16 +1012,16 @@ ! FILENAME_OBS = '/home/kjw/GEOS-Chem/gcadj_std/runs/v8-02-01/' // ! & 'ch4/geocape/' // GET_RES_EXT() // '/adjtmp/' // ! & 'gctm.obs.YYYYMMDD.hhmm' - FILENAME_OBS = 'gctm.obs.YYYYMMDD.hhmm' + FILENAME_OBS = 'gctm.obs.YYYYMMDD.hhmm' ! 观测文件, 居然是直接指定的么 CALL EXPAND_DATE( FILENAME_OBS, GET_NYMD(), GET_NHMS() ) FILENAME_OBS = TRIM( ADJTMP_DIR ) // TRIM( FILENAME_OBS ) XTAU = GET_TAU() CALL READ_BPCH2( TRIM(FILENAME_OBS), 'IJ-OBS-$', 1, & XTAU, IIPAR, JJPAR, - & LLPAR, DUMMY_TRUE, QUIET=.TRUE.) - GC_CH4_TRUE_ARRAY(:,:,:) = DUMMY_TRUE(:,:,:) + & LLPAR, DUMMY_TRUE, QUIET=.TRUE.) ! 读取观测数据 + GC_CH4_TRUE_ARRAY(:,:,:) = DUMMY_TRUE(:,:,:) ! 将观测数据赋值 - ! Convert from [kg] --> [v/v] + ! Convert from [kg] --> [v/v] 逐个网格转换单位 DO II=1,IIPAR DO JJ=1,JJPAR DO LG=1,LLPAR @@ -1031,23 +1031,23 @@ ENDDO ENDDO - ! Read a priori vertical profiles from file + ! Read a priori vertical profiles from file 读取先验廓线 FILENAME = '/home/kjw/new_satellites/geocape/data/' // & 'geocape_prior.' // GET_RES_EXT() // '.bpch' XTAU = GET_TAU0( 1, 1, 1985 ) CALL READ_BPCH2( FILENAME, 'IJ-AVG-$', 1, & XTAU, IGLOB, JGLOB, & LLGEOCAPE, DUMMY_PRIOR, QUIET=.TRUE. ) - CH4_PRIOR(:,:,:) = DUMMY_PRIOR(:,:,:) + CH4_PRIOR(:,:,:) = DUMMY_PRIOR(:,:,:) ! 给先验廓线赋值 - ! Select arbitrary II, JJ and NT value + ! Select arbitrary II, JJ and NT value 为什么这里的数值都直接给定了 II=40 JJ=JJPAR-10 NT=100 - ! Initialize variables + ! Initialize variables 初始化相关的数组 GC_PCENTER(:) = 0d0 GC_PEDGE(:) = 0d0 GC_AD(:) = 0d0 @@ -1079,12 +1079,12 @@ ENDDO - ! Number of vertical levels to use in these observations + ! Number of vertical levels to use in these observations 处理观测和模式垂直层上的问题 ! Chop off lowermost levels if ! GEOS-Chem surface pressure < GEOCAPE pressure levels nlev = count( PRESSURE_EDGE .LT. GC_PEDGE(1) ) IF ( nlev .LT. 13 ) nlev = nlev + 1 - lind = LLGEOCAPE + 1 - nlev ! minimum vertical index on GEOCAPE grid + lind = LLGEOCAPE + 1 - nlev ! minimum vertical index on GEOCAPE grid 垂直层相关 ! Get interpolation matrix that maps GEOS-Chem to GEOCAPE grid GRIDMAP(1:LLPAR, 1:LLGEOCAPE) = @@ -1092,11 +1092,11 @@ ! Get GEOS-Chem column from "truth" run to make pseudo-observations GC_CH4_NATIVE_OB(:) = 0d0 - GC_CH4_NATIVE_OB(:) = GC_CH4_TRUE_ARRAY(II,JJ,:) + GC_CH4_NATIVE_OB(:) = GC_CH4_TRUE_ARRAY(II,JJ,:) ! 给观测数据赋值 ! Interpolate GEOS-Chem CH4 column and observation to GEOCAPE grid - ! Column in [v/v] + ! Column in [v/v] 在垂直上进行循环, 进行插值 DO LN = lind, LLGEOCAPE GC_CH4_onGEOCAPE(LN) = 0d0 GC_CH4_onGEOCAPE_OB(LN) = 0d0 @@ -1122,8 +1122,8 @@ ! A = GEOCAPE averaging kernel !-------------------------------------------------------------- - ! x_m - x_a for model and "observation" - ! [v/v] --> ln( v/v ) happens here + ! x_m - x_a for model and "observation" 以下部分似乎是在计算模式柱浓度 + ! [v/v] --> ln( v/v ) happens here 这部分是在观测垂直层上计算的 DO LN = lind, LLGEOCAPE GC_CH4_onGEOCAPE(LN) =MAX(GC_CH4_onGEOCAPE(LN), 1d-10) GC_CH4_onGEOCAPE_OB(LN)=MAX(GC_CH4_onGEOCAPE_OB(LN),1d-10) @@ -1153,7 +1153,7 @@ ! For safety, initialize these up to LLGEOCAPE - ! Add random error to this observation + ! Add random error to this observation 给观测增加随机误差 CH4_HAT_werr(:) = CH4_HAT(:) DO LN = lind, LLGEOCAPE @@ -1171,59 +1171,59 @@ ! J = [ model - obs ]^T S_{obs}^{-1} [ model - obs ] !-------------------------------------------------------------- - ! Calculate difference between modeled and observed profile + ! Calculate difference between modeled and observed profile 计算差分部分 DO LN = lind, LLGEOCAPE DIFF(LN) = CH4_HAT_werr(LN) - CH4_HAT_OB(LN) ENDDO - - ! Calculate adjoint forcing: 2 * DIFF^T * S_{obs}^{-1} - ! and cost function: DIFF^T * S_{obs}^{-1} * DIFF - DO LN = lind, LLGEOCAPE - DO LLN = lind, LLGEOCAPE - FORCE(LN) = FORCE(LN) + + ! 因为卫星垂直上各层的误差存在相关性, 所以最终每层都要处理与其他各层之间的关联 + ! Calculate adjoint forcing: 2 * DIFF^T * S_{obs}^{-1} 计算伴随强迫 + ! and cost function: DIFF^T * S_{obs}^{-1} * DIFF 代价函数就是距离 + DO LN = lind, LLGEOCAPE ! 在垂直上进行循环 + DO LLN = lind, LLGEOCAPE ! 继续在垂直上进行循环? + FORCE(LN) = FORCE(LN) + ! 伴随强迫的每一层都加上了该层与其他层之间的协方差 & 2d0 * OBSERROR_INV(LN,LLN) * DIFF(LLN) ENDDO NEW_COST(LN) = NEW_COST(LN) + 0.5*DIFF(LN)*FORCE(LN) - ENDDO + ENDDO ! 也就是说, 这部分是将有限差分作为切线性算子的替代, 计算相关参数 !-------------------------------------------------------------- - ! Begin adjoint calculations + ! Begin adjoint calculations 开始计算伴随部分 !-------------------------------------------------------------- - + ! 整理以下这部分就是: FORCE = DIFF = CH4_HAT_werr_ADJ = CH4_HAT_ADJ ! The adjoint forcing is 2 * S_{obs}^{-1} * DIFF = FORCE - DIFF_ADJ(:) = FORCE(:) + DIFF_ADJ(:) = FORCE(:) ! 差分的伴随就是伴随强迫 ! Adjoint of GEOS-Chem - Observation difference - CH4_HAT_werr_ADJ(:) = DIFF_ADJ(:) + CH4_HAT_werr_ADJ(:) = DIFF_ADJ(:) ! 观测模拟差的伴随 - ! Adjoint of adding random error to observation + ! Adjoint of adding random error to observation 随机观测误差的伴随 DO LN=lind,LLGEOCAPE - CH4_HAT_ADJ(LN) = 0d0 + CH4_HAT_ADJ(LN) = 0d0 ! 初始化为 0 - DO LLN=lind,LLGEOCAPE + DO LLN=lind,LLGEOCAPE ! 考虑协方差, 计算最终误差 CH4_HAT_ADJ(LN) = CH4_HAT_ADJ(LN) + & CH4_HAT_werr_ADJ(LLN) * RANDNUM(NT) * OBSERROR(LLN,LN) ENDDO ENDDO - CH4_HAT_ADJ(:) = CH4_HAT_werr_ADJ(:) + CH4_HAT_ADJ(:) = CH4_HAT_werr_ADJ(:) ! 那你这上面计算的有什么意义 - ! Adjoint of GEOCAPE observation operator + ! Adjoint of GEOCAPE observation operator 观测算子的伴随 CH4_PERT_ADJ(:) = CH4_HAT_ADJ(:) DO LN=lind,LLGEOCAPE - CH4_PERT_ADJ(LN) = 0D0 + CH4_PERT_ADJ(LN) = 0D0 ! 又给初始化了, 那上面的赋值有啥意义 - DO LLN=lind,LLGEOCAPE + DO LLN=lind,LLGEOCAPE ! 反正和协方差计算相关, 累积 CH4_PERT_ADJ(LN) = CH4_PERT_ADJ(LN) + & AVGKERNEL(LLN,LN) * CH4_HAT_ADJ(LLN) ENDDO ENDDO - ! Adjoint of x_m - x_a - DO LN = lind, LLGEOCAPE - ! fwd code: - !GC_CH4(LN) = MAX(GC_CH4(LN), 1d-10) + ! Adjoint of x_m - x_a 状态变量的伴随? + DO LN = lind, LLGEOCAPE ! 对于观测垂直层进行循环 + ! fwd code: 看上去这边是针对正向模拟写的, 因为正向模拟这边存在范围限制 + !GC_CH4(LN) = MAX(GC_CH4(LN), 1d-10) 限制了最小值为 1e-10 !CH4_PERT(LN) = LOG(GC_CH4(LN)) - LOG(PRIOR(LN)) ! adj code: IF ( GC_CH4_onGEOCAPE(LN) > 1d-10 ) THEN @@ -1235,28 +1235,28 @@ ENDDO - ! Adjoint of interpolation - DO LN=lind,LLGEOCAPE - DO LG=1,LLPAR + ! Adjoint of interpolation 插值还需要伴随的么 + DO LN=lind,LLGEOCAPE ! LN 是观测的垂直层 + DO LG=1,LLPAR ! LG 是模式的垂直层 GC_CH4_NATIVE_ADJ(LG) = GC_CH4_NATIVE_ADJ(LG) + - & GRIDMAP(LG,LN) * GC_CH4_onGEOCAPE_ADJ(LN) + & GRIDMAP(LG,LN) * GC_CH4_onGEOCAPE_ADJ(LN) ! 从卫星垂直层插值到模式垂直层 ENDDO ENDDO - ! Adjoint of unit conversion + ! Adjoint of unit conversion 单位转换的伴随, 大概就是从原始值到体积比的换算 DO LG=1,LLPAR GC_CH4_NATIVE_ADJ(LG) = GC_CH4_NATIVE_ADJ(LG) & * XNUMOL(1) / ( XNUMOLAIR * GC_AD(LG) ) ENDDO - ! Pass adjoing forcing back to adjoint tracer array - DO LG=1,LLPAR + ! Pass adjoing forcing back to adjoint tracer array 将伴随强迫返回到大数组上 + DO LG=1,LLPAR ! 最终的伴随强迫是在模式的垂直层上 ADJ(LG) = GC_CH4_NATIVE_ADJ(LG) * CHK_STT(II,JJ,LG,1) ENDDO - ! Update cost function + ! Update cost function 更新代价函数 COST_FUNC_A = COST_FUNC_A + SUM(NEW_COST(:))