GEOS-Chem Adjoint User's Guide
Running the adjoint code
Selecting the adjoint model operational mode
Selection of adjoint model operation mode is done via an input file input.gcadj in the run directory and a preprocessor file define_adj.h in the code/adjoint/ directory.
Active variable selection
There are two options for active variables (the ones with respect to which we compute sensitivities): initial conditions (LICS) and emissions (LADJ_EMS), which could include chemical production sources. The active variable tracers (LICS) need to be listed in input.gcadj and their order and numbering needs to correspond to the tracers listed in the input.geos file. For LADJ_EMS, the emissions also need to be listed in input.gcadj.
Forward model settings
Not all (forward) GEOS-Chem options are accommodated in the adjoint model. Currently, the forward model settings that are supported are as follows:
lp12cml
resolution: & 4x5 and 2x2.5, nested CO capabilities
simulation types: & full chemistry, tagged CO, tagged O3 and offline CO<math>_2</math>.
meteorology: & GEOS3, GEOS4, GEOS5
chemistry: & via KPP Rosenbrock solver
The forward model needs to be set to one of the supported options, unless we set LADJ to FALSE (in input.gcadj), which will run the code in forward model only mode. Check the GEOS-Chem adjoint wiki for up-to-date list of current supported features of the forward model.
Adjoint code checklist
Generally and ideally all settings should be controlled via the input files (input.gcadj and define_adj.h). However, a few hardwired options remain and are listed below. Each needs to be checked before a new simulation.
Subroutine INIT_WEIGHT in adj_arrays_mod.f:
This subroutine is used to specify the spatial domain over which observations are included in the cost function for the PSEUDO_OBS 4D-Var tests as well as for the following sensitivity run cost function options: LKGBOX, LUGM3, LSTT_PPB, and LSTT_TROP_PPM. If you’re computing a sensitivity of specific gridbox or region, you need to make sure that the WEIGHT array is initialized accordingly.
Subroutines CALC_ADJ_FORCE_FOR_SENS
If you are using a cost function option that depends upon species concentrations (i.e., CSPEC), you can adjust the LMIN, LMAX, etc, parameters at the top of this routine to limit the spatial coverage of the cost function.
Makefile
Makefile settings have to correspond to settings in define_adj.h file for observational operators, i.e. if we’re using hdf or netcdf files, we need to include appropriate libraries. Alternatively, we need to make sure that the files requiring those libraries are not being compiled when they are not needed. For a list of all options supported type 'make help' inside the code/ directory.
subroutines SET_SF or SET_LOG_SF in inverse_mod.f.
The SF_DEFAULT values in input.gcadj can be used to set global values for the scaling factors during iteration X=1. Any non-global settings of initial guesses require changing the code in these subroutines.
Finite difference test checklist
[sec:fdchecklist] The finite difference test is either done 1 gridbox at a time (FD_SPOT is TRUE) or globally at once (FD_GLOB is TRUE). Finite difference test compares an adjoint gradient to its finite difference approximation. The finite difference perturbation, as well as the gridbox for the spot test can be set in input.gcadj. Current preferred method of testing adjoint code is to perform a global FD test.
Checklist:
Comment out all observation operators (including pseudo observations) in define_adj.h.
Set LSENS to TRUE, and both L4DVAR and L3DVAR to FALSE.
Set FD_GLOB or FD_SPOT to TRUE (not both).
Set LADJ_EMS or LICS to TRUE.
Set finite perturbation, FD_DIFF (e.g. to something like 0.1).
Select a numerator and denominator for the test in the FD MENU.
Decide whether to specify finite difference box with index (Specify box is TRUE) or lon/lat values.
Run from iteration 1 to iteration 3 for a global test . The first iteration performs the base case forward model run and the adjoint run. For iterations 2 and 3, the forward model will be run for positive and negative perturbations of the control variables being tested.
For a spot test, run from iteration 1 to 2. Both iterations include a forward model run and adjoint model run. The second run is evaluated with a positive perturbation. Finite difference sensitivities are compared to the average of the adjoints from the first and second run.
A few things to note about FD tests:
- First guesses settings SF_DEFAULT do not apply here.
- The default cost function setting for FD_GLOB is LKGBOX. It applies to the entire layer LFD.
- FD_SPOT can use any of the STT-based cost function unit options. It applies to grid cell IFD,JFD,LFD.
- If (FD_GLOB is TRUE), transport will automatically be turned off, overwriting settings in input.geos).
- Currently FD_GLOB is only setup to work with tracers, not yet species.
- For FD_SPOT, species-based cost functions can be implemented by using the CSPEC_PPB option and filling in the CSPEC OBS MENU. NFD will be ignored.
- Setting FD_GLOB will trigger LMAX_OBS = T and NSPAN = 1.
- Setting NFD will override any tracer that is listed in the OBSERVATION MENU.
Relevant diagnostic output:
diagadj/gctm.fd.YYYYMMDD.hhmmss Diagnostic file. Used for process specific finite difference tests. It contains 6 data blocks. The first two are gradients, the third one is the ratio of the gradients.
diagadj/gctm.fdglob.YYYYMMDD.hhmmss Diagnostic file. Used for process specific second order finite difference tests. Only generated if FD_GLOB is set to true. It contains 6 data blocks. The first two are gradients, the third one is the ratio of the gradients
Sensitivity (non-finite difference checklist)
The sensitivity run local and global sensitivity calculation. LSENS controls this simulation type. The options here are sensitivity of a predefined gridbox concentration, average concentration etc. (in CALC_ADJ_FORCING in geos_chem_adj_mod.f).
Note that if LSENS is TRUE, both L4DVAR and L3DVAR have to be FALSE.
Comment out all observation operators (including pseudo observations) in define_adj.h.
Set LSENS to TRUE, and both L4DVAR and L3DVAR to FALSE.
Set FD_GLOB and FD_SPOT to FALSE.
Set LADJ_EMS or LICS to TRUE.
Specify as control variables the tracers (for LICS) or emissions LADJ_EMS for which you would like to calculate sensitivities. For LADJ_EMS, set opt=T, otherwise the gradients get set to zero.
Specify observations / cost function options.
- Set OBS_FREQ to desired frequency of numerator calculation. Note that even if you want a monthly mean concentration, you still have to set OBS_FREQ to something like 60 (min).
- Use LMAX_OBS and NSPAN to evaluate your cost function over a limited time range.
- Select an option for the cost function evaluation, and recognize whether it depends upon tracers (STT) or species (CSPEC)
- List the tracers or species to be included in the cost function in input.gcadj.
- Modify the spatial domain of the cost function using WEIGHT for tracers or LMIN, LMAX etc. for species.
Run for iteration 1 only. This iteration performs the base case forward model run and the adjoint run.
4D-var checklist
[sec:4dvarchecklist] 4D-var run is an inverse model run, which takes pseudo or real observations to estimate initial conditions, emissions, chemical sources and/or reaction rates. Observations must be set in define_adj.h. All other flags are in input.gcadj. This simulation might require special data libraries (e.g. ncdf, hdf), data, and additional input files (see section 3.3). Supported active/control variables are LICS, LADJ_EMS (someday, both).
For inverse model tests using pseudo observations, we usually set the “initial guess” of scaling factors to be some value other than one (or zero for LOG_OPT) and then try to converge to values of one (zero). For real data assimilation problems or attainment studies, we begin with our best estimate, i.e. scaling factors equal to one (zero), and converge to values that lead to best agreement with observations.
Checklist:
Select desired observations in define_adj.h
Decide whether to optimize scaling factors (default) or their log, also in define_adj.h.
Depending on whether we are using observations requiring additional libraries (e.g, hdf or netcdf), modify the run script to use the appropriate Makefile.
Set a desired number of iterations in the run script if using a multiple iteration script.
Set L4DVAR to TRUE., check that LSENS and L3DVAR are FALSE
Select control parameters LADJ_EMS or LICS.
Set APSRC to TRUE if including a priori cost function term. If TRUE, make sure that REG_PARAM and ERROR are intentionally defined in the CONTROL VARIABLE MENU.
Set initial conditions for ICS_SF and EMS_SF. This can be done globally using SF_DEFAULT, or on a regionally specific basis by modifying code in inverse_mod.f.
Set observation frequency.
List optimized species and/or emissions.
If your observation operator depends upon a chemical species in CSPEC, such as O<math>_3</math> or NO<math>_2</math>, list this species as an observed species in the OBSERVATION MENU.
Specify for which gridbox or lat/lon you want debugging print statements in the finite difference menu.
Select diagnostic output (in the diagnostics menu) such as LITR.
3D-var checklist
3D-var run is a Kalman filter forward estimation. The code exists outside the v8 framework. Contact Kumaresh Singh if interested in details.
Coding, debugging and testing
Sensitivity with respect to reaction rate coefficients: adding new reactions
[sec:arr]
Currently, sensitivities with respect to the following reaction rate constants are included:
To have the model calculate sensitivity with respect to additional rate constants, the following needs to be updated:
- Update the total number of active reaction rate constants, NCOEFF in gckpp_adj_Global.f90. NCOEFF is a the size of the vector JCOEFF, which holds the indicies of the active reaction rate constants.
- For the new active rate constants, make an entry in gckpp_adj_Util.f90 for JCOEFF. Define the value of JCOEFF(new index) to be the index of the reaction rate in the KPP generated reaction rate constant vector RCONST. This number can be found in the KPP input file gckpp_adj.eqn on left side of each reaction definition. The corresponding SMVGEAR index is on the right side.
- Add an entry in the mapping array RCONST2RKPP for the new reaction in gckpp_adj_Util.f90.
- Make sure that the sensitivity with respect to this rate constant is actually being calculated in the subroutine ros_DadjInt in gckpp_adj_Integrator.f90. The code contains the general formula for calculating discrete sensitivities with respect to reaction rate constants, see Appendix A of . The general formula calculates all elements of <math>(J_p(t^n,c^n) \times k_i)^T \cdot u_i</math> and <math>f_p^T(T_i,C_i)\cdot u_i</math>. However, many of the elements of <math>J_p(t^n,c^n) \times k_i)^T </math> (DJDR) and <math>f_p^T(T_i,C_i)</math> (DFDR) are actually zero. Automatically exploiting such sparsity will be a future feature of KPP. Until then, we have to enter the non zero elements manually if we want only the sparse calculation.
For example, if we would like to calculate the sensitivity of the cost function with respect to the rate constant for the reaction ALD<math>_2</math> + NO<math>_3 </math><math>\rightarrow</math> HNO<math>_3</math> + MCO<math>_3</math>, we would find the following entry for this reaction in the KPP input file gckpp_adj.eqn:
{40} ALD2 + NO3 = HNO3 + MCO3: RKPP(JLOOP,54);
If NCOEFF was previously 8, then we would set NCOEFF = 9 in gckpp_adj_Global.f90. In gckpp_adj_Util.f90, in the subroutine INIT_KPP we would add the following lines:
JCOEFF(9) = 40 ! ALD2 + NO3 --> HNO3 + MCO3
RCONST2RRATE(40) = 54
As mentioned above, in the routine ros_DadjInt in gckpp_adj_Integrator.f90, these sensitivities can be calculated using the general formula, or one could implement the sparse formula:
CASE(9) ! f_r(y) * u j = ind_ALD2 V_P(vpstart+icoeff,m) = DFDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) j = ind_NO3 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DFDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) j = ind_HNO3 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DFDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) j = ind_MCO3 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DFDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) ! ( J_r x k )^T * u j = ind_ALD2 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) j = ind_NO3 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) j = ind_HNO3 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) j = ind_MCO3 V_P(vpstart+icoeff,m) = V_P(vpstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m)
Saving reaction rate sensitivities
The subroutine SAVE_INST_ARR can be used to print out reaction rate sensitivities in a particular cell at every time step. This routine also calls the routine MAKE_ARR_FILE which writes the integrated sensitivities globally to a binary punchfile, gctm.arr.
Troubleshooting and debugging
If you find any incompatible set of switches in input.gcadj that crash the code, please add corresponding error checks in ARE_FLAGS_VALID in input_adj_mod.f.
Validating code
[sec:val]
The adjoint of GEOS-Chem is designed to be validated in several manners.
Global tests of a subset of the adjoint model
By turning off multidirectional transport related processes, the adjoint model sensitivities can be compared to finite difference sensitivities on a global scale, see Fig. 1 and 3 of .
For example, to check the adjoint of the chemistry only in a single vertical level, set the following:
in run:
X=1
XSTOP=3
TYPE=DEFAULT
in input.gcadj:
LSENS = T
FD_GLOB = T
LICS or LADJ_EMS or LADJ_EMS and LADJ_STRAT = T
select a dependent species, NFD
select a vertical level, LFD
select a control parameter, ICSFD or EMSFD and/or STRFD
IFD and JFD don’t matter as we’re doing a domain wide test
set the finite difference perturbation (<math>\delta\sigma</math>), FD_DIFF = 1.d-1
Set the OBS_FREQ to 60, but set LMAX_OBS = T and NSPAN = 1.
Set LAERO_THERM = .FALSE.
in input.geos:
Turn off convection ( LCONV = F )
Turn off turbulent mixing ( LTURB = F )
Turn off wet deposition ( LWETD = F )
Leave on dry deposition ( LDRYD = T )
Can leave on transport ( LTRAN = T ), which will be overridden if FD_GLOB = T.
in chemistry_mod.f, can save time by computing gas-phase chemistry only in certain cells by placing an IF statement around CALL INTEGRATE_ADJ such as:
IF ( L == LFD ) THEN
CALL INTEGRATE_ADJ(...
ENDIF
Analysis of global FD tests is described in Sect. [sec:geos5bench].
Spot tests of full adjoint model
Alternatively, we can compare adjoint gradients to finite difference gradients for control parameters one location at a time, but with all model processes turned on.
- in run:
- X=1
- XSTOP=2
- TYPE=DEFAULT
- in input.gcadj:
- LSENS
- Make sure that FD_GLOB is set to FALSE and that FD_SPOT is set to TRUE.
- active variables are LICS or LADJ_EMS or LADJ_EMS and LADJ_STRAT
- select a dependent species, NFD
- select particular control parameter IFD JFD LFD and ICSFD or EMSFD and/or STRFD
- set the finite difference perturbation (<math>\delta\sigma</math>), FD_DIFF = 1.d-1
- probably want to turn on additional diagnostic output, so set L_PRINTFD = .TRUE.
- Make sure that OBS_FREQ is set long enough so that the cost function is only evaluated once during the simulation.
- in input.geos:
- Turn on all desired processes.
The adjoint and finite difference gradients and the ratio ADJ / FD will be written to standard output at the end of the run. The reported adjoint gradient is actually the average of the values at <math>\sigma = 1 </math> and <math>\sigma = 1 + \delta\sigma</math>. Since such tests involve the continuous adjoint of advection, the ratio ADJ / FD can not be expected to be unity. Benchmark tests are given in .