View source for GEOS-Chem Adjoint User's Guide

On this page we describe how to download and use the GEOS-Chem adjoint model.

= Getting started =

== Brief overview ==

The GEOS-Chem adjoint model is an adjoint model derived from the GEOS-Chem CTM. Although in terms of GEOS-Chem activities, it is one of many projects in terms of code itself, it is a super-structure, containing both, the forward GEOS-Chem and its derivative adjoint code. Great effort is made to keep the adjoint current with updates in the GEOS-Chem, which have to be implemented manually. Currently, there is not an adjoint equivalent of every part of GEOS-Chem.

The FORTRAN code that is GEOS-Chem adjoint can perform a number of calculations: sensitivity calculations (most efficient if ''y'' in  <span class="math">∂<em>y</em> / ∂<em>x</em></span> is a scalar, and ''x'' is a 2-, 3- or 4-dimensional field). The adjoint code can also be used for inverse problems, although some code development might be required to interface with observational datasets. Currently, some observation operators are available for several species.

Original work on the adjoint of GEOS-Chem began in 2003, focusing on the adjoint of the offline aerosol simulation. By 2005, the adjoint was expanded to include a tagged CO simulation and a full chemistry simulation ''(Kopacz, et al. 2009a,b; Henze, et al. 2007,2009)''; an adjoint of GEOS-Chem v7 was also developed in the following years'' (Zhang et al. 2009; Singh et al. 2009)''. Each of these branches of the adjoint code were been constructed in a hybrid fashion using a combination of automatic differentiation software (TAMC, KPP) and manual coding of both discrete and continuous adjoints. They shared many common elements yet had unique features for different applications. During the summer of 2009, the existing branches were merged and updated to bring the adjoint into alignment with the latest release of GEOS-Chem, v8-02-01. This merged adjoint model is now the standard adjoint code into which all further development efforts will be placed.

The adjoint model is maintained by a group of its users. The Adjoint Model Scientist is Prof. Daven K. Henze at University of Colorado, who is also in charge of the repository and code distribution. Questions regarding this manual and code in general can be directed to him (daven.henze@colorado.edu), as well as to Monika Kopacz (mkopacz@princeton.edu) and Kumaresh Singh (kumaresh@vt.edu). As the adjoint code is still in the early phase of development and testing, if you wish you publish your work done with the standard adjoint code described in this manual, we request that you offer co-authorship to the following group of core developers:

Monika Kopacz (mkopacz@princeton.edu) 

Daven K. Henze (daven.henze@colorado.edu)

Kumaresh Singh (kumaresh@vt.edu)

Changsub Shim (changsub.shim@jpl.nasa.gov)

== Recent and ongoing updates ==

See the wiki <tt>(http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_Adjoint)</tt> for a complete list of features that are implemented and/or in the process of being updated to the GC v8 adjoint.

== Obtaining the adjoint model ==

Model packages including source code and run directories are available through the GIT repository. This [http://wiki.seas.harvard.edu/geos-chem/index.php/Quick_Start_Guide. Quick Start Guide] goes through the process of how to create an account get the latest version of the code.

More information on GIT is at '' [http://git-scm.com/documentation GIT Documentation] '' and '' [http://www.kernel.org/pub/software/scm/git/docs/git.html GIT Man Page] '' . Please review this online documentation prior to using GIT. Once you have familiarized yourself with the sorfware you may proceed.

==== Initial download ====

To download the most current copy of the code, enter the following command

<tt>git clone http://adjoint.colorado.edu:8080/gcadj_std.git</tt>

This will download a project directory ''“gcadj_std”'' containing the source code and run directories.

==== Tracking subsequent changes by yourself and others ====

Perhaps you already have a version of model. To determine the status of your existing project vs the current repository version, enter your local copy of your project directory and type

<tt>git status</tt>

To determine the difference between your local copy and the current repository copy, type

<tt>git diff --word-diff=color origin/master</tt>

The ''—word-diff=color ''option makes the output with colors so it’s easy to read, specifiying ''origin/master'' takes the difference between your copy and the newest repository copy. Without ''origin/master'', you will see the difference between your copy and the version that was in the repository as of the last time you checked out the code.

The commands above are without arguments and will thus apply to all files in the project. To use them for only one specific file, for example ''geos_chem_adj_mod.f'', type:

<tt>git diff --word-diff=color origin/master -- geos_chem_adj_mod.f</tt>

To replace your local copy of ''geos_chem_adj_mod.f'' with the newest version from the repository,

<tt>git checkout origin/master -- geos_chem_adj_mod.f</tt>

If you have also modified your local copy, you can merge your changes with those made to the repository version

<tt>git checkout --patch origin/master -- geos_chem_adj_mod.f</tt>

If you want your updates to be added to the main repository please create a [http://wiki.seas.harvard.edu/geos-chem/index.php/Using_Git_with_GEOS-Chem#Sharing_your_revisions_with_others_.28and_vice_versa.29 ''PATCH''] and contact the [http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_Adjoint#Contact_information Support Team]. Daven Henze will handle incorporating changes into the code repository. Please do not distribute copies of the code at this time; you may refer such requests to Daven.

== Additional files for analysis ==

Some IDL and MATLAB scripts for plotting results of finite difference tests are here:

<tt>http://adjoint.colorado.edu/~daven/gcadj_std/tools.tar.gz</tt>

== Benchmark simulations ==

==== Setting input files ====

You will need to make the following changes to the input files, as they depend upon your particular filesystem:

* ''run'' script
** Set <tt>DRUN</tt> and <tt>DSAVE</tt> in the ''run'' script, as dictated by your filesystem.
** Depending upon your computer system, you may or may not need to include coping of you project directory to the $DSAVE filesystem after each iteration. If you won’t lose access to the local filesystem where you model is running after it has completed executing, then you don't need to set <tt>DSAVE</tt>.
** If you are running on multiple CPU cores, be mindful of the line  <tt>export OMP_NUM_THREADS=24</tt>  and adjust as appropriate for your system.
* Set data folder locations in ''input.geos'' following standard forward model procedures.
* Set NetCDF locations in ''Makefile''
* Set HDF locations in ''Makefile'' (if your simulation requires these)

==== Setting source code ====

The source code is set to run the ''geos5'' benchmark. To use the ''geos4'' benchmark,

* change the preprocessor flag for meteorology from GEOS_5 to GEOS_4 in ''code/define.h''
* disable the IN_CLOUD_OD preprocessor flag in ''code/define.h''.
* enable the PSEUDO_OBS preprocessor flag in ''code/adjoint/define_adj.h''.

==== Running the tests ====

Now you are ready to execute the run script, with a unix command like:

<tt>./run &gt; log.my_benchmark_test &amp;</tt>

If all goes well, the output should finish with:

<pre> ------------------------------------------------
 G E O S   C H E M   A D J O I N T   E X I T E D
            N O R M A L L Y
 ------------------------------------------------</pre>
=== The geos5 full chemistry finite difference test ===

[sec:geos5<sub>b</sub>ench] The ''geos5'' benchmark simulation is a full chemistry sensitivity test. It checks the sensitivity of O<math>_x</math> with respect to NO<math>_x</math> using adjoints and finite difference calculations. The results are in the ''diagadj/*.fdglob.*'' file.

To speed up the evaluation of this test, chemistry is only calculated in the <tt> LFD</tt> level.

This simulation will run for 1 complete iteration, and then only the forward part of the 2nd and 3rd iterations, after which the log file will state something like:

<pre>  Global validation test for values &gt;    1.00000000000000
  MAX of global 2nd order ADJ / FD  =    3.411938
  MIN of global 2nd order ADJ / FD  =  -6.9405900
  Number of places where ratio off by   0.100000000000000    =      88
 </pre>
That the min and max of the ADJ / FD ratio is much greater or less than 1.d0 may or may not be significant. In general, it shouldn’t be off by more than x10, and the number of places where the ratio is off by more than 10% shouldn’t be much more than 100. However, to tell whether ADJ / FD ratios that deviate from 1.0 are owing to errors in the adjoint code or to limitations of finite difference sensitivities, the results will have to be plotted.

The key results are saved to the ''*.fdglob.*'' file. The contents can be analyzed using GAMAP. You can load a file and view results such as finite sensitivities, adjoint sensitivities, and the ratio of the two.

<pre>IDL&gt; gamap, file='gctm.fdglob.20050701.0500' </pre>
Next generate a data set for a scatter plot of adjoint vs finite difference values using (thanks to Kevin Wecht for providing this IDL script):

<pre>IDL&gt;plot_fdglob</pre>
Alternatively, if you prefer to generate figures in MATLAB, scripts for the following scheme are also provided. First output the sensitivities to some text files.

<pre>IDL&gt; fd_stats</pre>
and then plot the results in MATLAB using

<pre>&gt;&gt;fd_vs_adj</pre>
Sample output is shown here:

<tt>http://spot.colorado.edu/~henzed/GC_adj/fd_validation.pdf</tt>

More about setting up and designing your own finite difference validation tests is in Sect. [sec:fd<sub>c</sub>hecklist].

=== The geos4 tagged CO optimization test ===

The optimization benchmark simulation attempts to optimize initial concentrations of CO, starting with domain wide linear scaling factors of 0.5 and pseudo observations that were generated with scaling factors of 1.00. The simulation is only 1 day long and should start to converge after a few iterations.

Output from our benchmark runs are in the directory ''../OptData/''. To quickly see how the cost function converges,

<tt>grep '' OptData/cfn*</tt>

The cost function should have reduced by about an order of magnitude. Note: not all value of the cost function listed here correspond to accepted iterations of the optimization procedure. Some are ‘function evaluations’ that correspond to the optimization performing searches in various directions before finding the optimal path towards a minimum. To see which values correspond to accepted iterations,

<tt>grep 'iterate' log</tt>

To see the inverse modeling solution after the 6th iteration, look at the optimized scaling factors in the ''gctm.sf.06'' file using gamap,

<pre>IDL&gt; gamap, file='gctm.sf.06', 'IJ-ICS-$', yrange=[0,1] </pre>
Scaling factors should be close to 1.0 in locations where CO concentrations are significant.

= Directories and input/output files =

== Directory structure ==

[sec:dir] The adjoint code package ''gcadj_std'' contains the following directories:

l p8cml

''code/'' &amp;Code directory, which contains all unmodified GEOS-Chem files, a Makefile and a few subdirectories relevant to the adjoint code, listed below. This is also where all the object files are placed.<br />
''code/modified/'' &amp;Subdirectory that contains all (forward) GEOS-Chem files that have been modified for the adjoint.<br />
''code/adjoint/'' &amp;Subdirectory that contains all adjoint specific files.<br />
''code/obs_operators/'' &amp;Subdirectory contains all files relevant to observation operators.<br />
''code/new/'' &amp; Subdirectory contains all new files (new to forward AND adjoint).<br />
''runs/'' &amp; Run directory with subdirectories for each met field type.<br />
''runs/v8-02-1/geos*/'' &amp; Run directories with subdirectories for output files.<br />


''runs/../geos*/adjtmp/'' &amp; Adjoint temporary file directory (''gctm.chk.*; gctm.obs.*; gctm.adj.*'' ), the name can be changed in ''input.gcadj''.<br />
''runs/../geos*/tmp/'' &amp; Forward model temporary file directory (unzipped met fields).<br />
''runs/../geos*/OptData/'' &amp; Results from each iteration (''gctm.gdt.*; cfn.*; gctm.sf.*, fwd_dat.*.tar; gctm.arr''), the name can be changed in ''input.gcadj''.<br />
''runs/../geos*/diagadj/'' &amp; Adjoint diagnostic files (''*.fd.*;*.fdglob.*; aero.ave.*; satave.*; jsave.*; gctm.iteration, emis.adj.*''), the name can be changed in ''input.gcadj''.<br />


== Input files ==

=== ''input.gcadj'': the main adjoint input file ===

''input.gcadj'' file contains the following options:

<pre>01: %%% ADJOINT SIMULATION MENU %%%
02: Do adjoint run    LADJ         : T
03: Selecet one simulation type    :---
04: Invese problem    L4DVAR       : T
05: Kalman filter     L3DVAR       : F
06: Sensitivity       LSENS        : F
07:  =&gt; spot finite diff  FD_SPOT  : F
08:  =&gt; global finite diff FD_GLOB : F</pre>
''Description''

l p12cm l

01: header line &amp; –<br />
02: ''LADJ'' &amp; Global switch for adjoint option. If set to FALSE, it will overwrite all other options and make the run a forward mode only run.<br />
03: line&amp; If LADJ = T, need to pick one of the following options. 3DVAR not yet supported.<br />
04: ''L4DVAR'' &amp; Switch for 4d-var runs.<br />
05: ''L3DVAR'' &amp; Switch for Kalman filter.<br />
06: ''LSENS'' &amp; Switch for sensitivity runs. If performing a finite difference test, pick one option from below:<br />
07: ''FD_SPOT'' &amp; Switch for spot finite difference test.<br />
08: ''FD_GLOB'' &amp; Switc for global finite difference test).<br />


<pre>01: %%% FORWARD MODEL OPTIONS %%%
02: adjoint chemistry LADJ_CHEM    : T
03: aerosol thermo    LAERO_THEM   : T</pre>
''Description''<br />


l p12cm l

01: header line &amp; –<br />
02: ''LADJ_CHEM'' &amp; Switch for adjoint chemistry option. If set to FALSE, it will turn off adjoint chemistry. Make sure that LCHEM is set to the same value.<br />
03: ''LAERO_THEM'' &amp; Switch for aerosol thermodynamics option (applies to forward and adjoint).<br />


<pre>01: %%% ADJOINT MODEL OPTIONS %%%
02: Include a priori term APSRC    : F
03: compute bkgrnd error cov       : F
04: Compute DFP inverse Hessian    : F
05: Compute BFGS inverse Hessian   : F
06: Include rxn rate sensitivities : F
07: Delete chk files LDEL_CHKPT    : T
08: Scale up and FILL adj transport: F</pre>
''Description''<br />


l p12cm l

01: header line &amp; –<br />
02: ''APSRC'' &amp; Switch for calculating the a priori term in the cost function. Valid option for 4d-var runs.<br />
03: ''COV'' &amp; Switch for computing the full error covariance matrices (not yet implemented).<br />
04: ''LINVH'' &amp; Switch for computing an approximate inverse Hessian matrix using the DFP method.<br />
04: ''LBFGS'' &amp; Switch for computing an approximate inverse Hessian matrix using the L-BFGS method.<br />
06: Rxn sensitivities &amp; Switch for storing sensitivities wrt reaction rates.<br />
07: LDEL_CHKPT &amp; Delete checkpoint files after they are used in adj run. Set to F to reuse them for multiple adj runs.<br />
08: LFILL_ADJ &amp; Scale up adjoints and then use the LFILL option in tpcore for advection.<br />


<pre>01: %%% DIRECTORIES %%%
02: Optimization output            : OptData/
03: Temporary adjoint dir adjtmp   : adjtmp/
04: Diagnostics ouptut             : diagadj/</pre>
''Description''<br />


l p12cm l

01: header line &amp; –<br />
02: OptData. &amp; Specify output directory where essential output will go, typically set to OptData.<br />
03: adjtmp. &amp; Specify output directory where nonessential output will go, typically set to adjtmp.<br />
04: diagadj &amp; Specify output directory where diagnostic output will go, typically set to diagadj.<br />


<pre>01: %%% CONTROL VARIABLE MENU %%%
02: Initial conditions LICS        : F
03: ... OR emissions   LADJ_EMS    : T
04:  =&gt; strat prod/loss LADJ_STRAT : T
05: &gt;------------------------------&lt; 
06:           FOR LICS             :
07: NSOPT: number of tracers opt   : 1
08:   =&gt; opt these tracers------&gt;  : TRC# trc_name SF_DEFAULT REG_PARAM ERROR
09: Tracer  #1                     : 1    NOx      1          1          1
10: &gt;------------------------------&lt;
11:           FOR LADJ_EMS         :
12: NNEMS: ems groups implemented  : 33
13: Emission entries ------------&gt; : EMS# ems_name       opt  SF_DEFAULT REG_PARAM ERROR
14: Emission #1                    : 1    IDADJ_ENH3_an  T    1          1         1
15: Emission #2                    : 2    IDADJ_ENH3_na  T    1          1         1
16: Emission #3                    : 3    IDADJ_ENH3_bb  T    1          1         1
17: Emission #4                    : 4    IDADJ_ENH3_bf  T    1          1         1
18: Emission #5                    : 5    IDADJ_ESO2_an1 T    1          1         1
19: Emission #6                    : 6    IDADJ_ESO2_an2 T    1          1         1
20: Emission #7                    : 7    IDADJ_ESO2_bf  T    1          1         1
21: Emission #8                    : 8    IDADJ_ESO2_bb  T    1          1         1
22: Emission #9                    : 9    IDADJ_ESO2_sh  T    1          1         1
23: Emission #10                   : 10   IDADJ_EBCPI_an T    1          1         1
24: Emission #11                   : 11   IDADJ_EBCPO_an T    1          1         1
25: Emission #12                   : 12   IDADJ_EOCPI_an T    1          1         1
26: Emission #13                   : 13   IDADJ_EOCPO_an T    1          1         1
27: Emission #14                   : 14   IDADJ_EBCPI_bf T    1          1         1
28: Emission #15                   : 15   IDADJ_EBCPO_bf T    1          1         1
29: Emission #16                   : 16   IDADJ_EOCPI_bf T    1          1         1
30: Emission #17                   : 17   IDADJ_EOCPO_bf T    1          1         1
31: Emission #18                   : 18   IDADJ_EBCPI_bb T    1          1         1
32: Emission #19                   : 19   IDADJ_EBCPO_bb T    1          1         1
33: Emission #20                   : 20   IDADJ_EOCPI_bb T    1          1         1
34: Emission #21                   : 21   IDADJ_EOCPO_bb T    1          1         1
35: Emission #22                   : 22   IDADJ_ENOX_so  F    1          1         1
36: Emission #23                   : 23   IDADJ_ENOX_li  F    1          1         1
37: Emission #24                   : 24   IDADJ_ENOX_ac  F    1          1         1
38: Emission #25                   : 25   IDADJ_ENOX_an  F    1          1         1
39: Emission #26                   : 26   IDADJ_ENOX_bf  F    1          1         1
40: Emission #27                   : 27   IDADJ_ENOX_bb  F    1          1         1
41: Emission #28                   : 28   IDADJ_ECO_an   F    1          1         1
42: Emission #29                   : 29   IDADJ_ECO_bf   F    1          1         1
43: Emission #30                   : 30   IDADJ_ECO_bb   F    1          1         1
44: Emission #31                   : 31   IDADJ_EISOP_an F    1          1         1
45: Emission #32                   : 32   IDADJ_EISOP_bf F    1          1         1
46: Emission #33                   : 33   IDADJ_EISOP_bb F    1          1         1
47: Number emis time group MMSCL   : 1
48: &gt;------------------------------&lt;
49:         FOR LADJ_STRAT         :
50: NSTPL: strat prod &amp; loss trcs  : 24
51: Strat prod &amp; loss trc entries  : ID# trc_name    opt   SF_DEFALUT  REG_PARAM ERROR
52: Tracer  #1                     : 1    NOx_p      T     1           1         1
53: Tracer  #2                     : 2    Ox_p       T     1           1         1
54: Tracer  #3                     : 3    PAN_p      T     1           1         1
55: Tracer  #4                     : 4    CO_p       T     1           1         1
56: Tracer  #5                     : 5    ALK4_p     T     1           1         1
57: Tracer  #6                     : 6    ISOP_p     T     1           1         1
58: Tracer  #7                     : 7    HNO3_p     T     1           1         1
59: Tracer  #8                     : 8    H2O2_p     T     1           1         1
60: Tracer  #9                     : 9    ACET_p     T     1           1         1
61: Tracer  #10                    : 10   MEK_p      T     1           1         1
62: Tracer  #11                    : 11   ALD2_p     T     1           1         1
63: Tracer  #12                    : 12   RCHO_p     T     1           1         1
64: Tracer  #13                    : 13   MVK_p      T     1           1         1
65: Tracer  #14                    : 14   MACR_p     T     1           1         1
66: Tracer  #15                    : 15   PMN_p      T     1           1         1
67: Tracer  #16                    : 16   PPN_p      T     1           1         1
68: Tracer  #17                    : 17   R4N2_p     T     1           1         1
69: Tracer  #18                    : 18   PRPE_p     T     1           1         1
70: Tracer  #19                    : 19   C3H8_p     T     1           1         1
71: Tracer  #20                    : 20   CH2O_p     T     1           1         1
72: Tracer  #21                    : 21   C2H6_p     T     1           1         1
73: Tracer  #22                    : 22   N2O5_p     T     1           1         1
74: Tracer  #23                    : 23   HNO4_p     T     1           1         1
75: Tracer  #24                    : 24   MP_p       T     1           1         1
76: Tracer  #1                     : 1    NOx_l      T     1           1         1
77: Tracer  #2                     : 2    Ox_l       T     1           1         1
78: Tracer  #3                     : 3    PAN_l      T     1           1         1
79: Tracer  #4                     : 4    CO_l       T     1           1         1
80: Tracer  #5                     : 5    ALK4_l     T     1           1         1
81: Tracer  #6                     : 6    ISOP_l     T     1           1         1
82: Tracer  #7                     : 7    HNO3_l     T     1           1         1
83: Tracer  #8                     : 8    H2O2_l     T     1           1         1
84: Tracer  #9                     : 9    ACET_l     T     1           1         1
85: Tracer  #10                    : 10   MEK_l      T     1           1         1
86: Tracer  #11                    : 11   ALD2_l     T     1           1         1
87: Tracer  #12                    : 12   RCHO_l     T     1           1         1
88: Tracer  #13                    : 13   MVK_l      T     1           1         1
89: Tracer  #14                    : 14   MACR_l     T     1           1         1
90: Tracer  #15                    : 15   PMN_l      T     1           1         1
91: Tracer  #16                    : 16   PPN_l      T     1           1         1
92: Tracer  #17                    : 17   R4N2_l     T     1           1         1
93: Tracer  #18                    : 18   PRPE_l     T     1           1         1
94: Tracer  #19                    : 19   C3H8_l     T     1           1         1
95: Tracer  #20                    : 20   CH2O_l     T     1           1         1
96: Tracer  #21                    : 21   C2H6_l     T     1           1         1
97: Tracer  #22                    : 22   N2O5_l     T     1           1         1
98: Tracer  #23                    : 23   HNO4_l     T     1           1         1
99: Tracer  #24                    : 24   MP_l       T     1           1         1</pre>
''Description''<br />


l p12cm l

01: header line &amp; Need to pick one of the three possible sets of control parameters: initial conditions, emissions, or stratospheric production and loss rates. Emissions must be turned on to select stratospheric fluxes.<br />
02: LICS &amp; Tracer initial conditions as control parameters.<br />
03: LADJ_EMS &amp; Emissions as control parameters.<br />
04: LADJ_STRAT &amp; Stratospheric production and loss rates as control parameters.<br />
05: spacer line &amp; -<br />
06: FOR LICS &amp; Specify which tracers to allow as control parameters. Note: the range of possible tracers is defined in input.geos. The adjoint will always include adjoints of all tracers. So here we just need to list which of these tracers will be optimized. All and only those tracers listed below will be optimized. This LICS section of the MENU will be ignored if LADJ_EMS = T.<br />
07: NSOPT &amp; Total number of tracers to optimize listed in the submenu below.<br />
08: subheader &amp; -<br />
09: Tracers #1 &amp; List the corresponding tracer number (TRC#) and name (trc_name) from input.goes. Here you can also specify a global default scaling factor (SF_DEFAULT) for the first iteration, a regularization parameter (REG_PARAM) and an error (ERROR). The latter two only have an effect if LAPSRC = T.<br />
…&amp; add more lines like 08 if you want to make the initial conditions for more than one tracer active. …<br />
10: spacer line &amp; -<br />
11: FOR LADJ_EMS &amp; Specify all of the emissions adjoints that are currently implemented. This LADJ_EMS section of the menu will be ingnored if LICS = T.<br />
12: NNEMS &amp; Total number of active emissions groups listed in the submenu below.<br />
13: subheader &amp; -<br />
14: Emission #1 &amp; List the emission number (EMS#) and name (ems_name). Names must begin with IDADJ_E…. Select wether this emissions group is to be optimized (opt). Here you can also specify a global default scaling factor (SF_DEFAULT) for the first iteration, a regularization parameter (REG_PARAM), and an error. The latter two only have an impact of APSRC = T.<br />
15…46 &amp; additional emission group definitions and options.<br />
47: MMSCL &amp; Number of emissions sub-scaling groups. MMSCL <math>></math> 1 not yet supported.<br />
48: spacer line &amp; -<br />


l p12cm l

49: FOR LADJ_STRAT &amp; Specify all of the stratospheric production and loss rate adjoints that are currently implemented. This LADJ_STRAT section of the menu will be only considered if both LADJ_EMS = T and LADJ_STRAT = T.<br />
50: NSTPL &amp; Total number of active stratospheric tracers that have production and loss rates listed in the menu below. NOTE: each tracer has production and loss rates.<br />
51: subheader &amp; -<br />
52: Tracer #1 &amp; List the corresponding tracer number (ID#) and name (trc_name). Names must end with …_p (production) or…_l (loss). Select whether this stratospheric flux is to be optimized (opt), although optimization of the stratospheric fluxes has yet to be fully tested. Here you can also specify a global default scaling factor (SF_DEFAULT) for the first iteration, a regularization parameter (REG_PARAM), and an error. For now, only SF_DEFAULT is supported.<br />
53…99 &amp; additional stratospheric flux definitions and options.<br />


<pre>01: %%% OBSERVATION MENU %%%
02: %%%  for PSUEDO_OBS  %%%
03: %%%   or LSENSE      %%%
04: Observation frequency OBS_FREQ : 60
05: Limit number of observations?  : F
06:   =&gt;  max number of obs NSPAN  : 1
07: COST FUNCTION options for LSENS:---
08:   =&gt;  tracer kg/box            : T
09:   =&gt;  tracer ug/m3             : F
10:   =&gt;  tracer ppb               : F
11:   =&gt;  tracer ppm free trop     : F
12:   =&gt;  species ppb w/averaging  : F
13:   =&gt;  tracer ug/m3 pop weight  : F
14: &gt;------------------------------&lt; 
15: NOBS: number of tracers to obs : 2
16:  =&gt; obs these tracers------&gt;  : TRC# tracer_name
17: Tracer  #1                     : 34   BCPI
18: Tracer  #2                     : 35   OCPI
19: &gt;------------------------------&lt;
20: NOBS_CSPEC: # of species to obs: 0
21:  =&gt; obs these species------&gt;  :species_name
22: Species #1                     : O3</pre>
''Description''<br />


l p12cm l

01: header line &amp; The options pretty much pertain only to sensitivity calculations or pseudo observation tests.<br />
02: header line &amp; The only exception to that is if you have an observation operator specific to a chemical species.<br />
03: header line &amp; which is not a tracer (e.g., O<math>_3</math>), in which case you need to specify it in the observed species submenu.<br />
04: ''OBS_FREQ'' &amp; Frequency (in min) of checking and assimilating observations, both pseudo and real, typically 60.<br />
05 LMAX_OBS: &amp; Set this if you wish to limit the number of observations. For example, if you want the cost function to be evaluated only during the final day of your simulation, and OBS_FREQ = 60, then set LMAX_OBS = T and NSPAN = 24. Setting FD_GLOB will trigger LMAX_OBS = T and NSPAN = 1.<br />
06: NSPAN &amp; If LMAX_OBS = T, then use this to set the number of times the cost function is evaluated. Setting FD_GLOB will trigger LMAX_OBS = T and NSPAN = 1.<br />
07: subheader &amp; Below are some options for evaluating the cost function during a sensitivity run. Note the distinction between tracers (STT) and species (CSPEC). Some of these options include the WEIGHT array, which allows for spatial masking. Check the respective code segments for details.<br />
08: LKGBOX &amp; Evaluate the cost function for tracer concentrations in units of kg/box. Note: FD simulations will default to this option.<br />
09: LUGM3 &amp; Evaluate the cost function for tracer concentrations in units of ug/m3.<br />
10: LSTT_PPB &amp; Evaluate the cost function for tracer concentrations in units of ppb.<br />
11: LSTT_TROP_PPM &amp; Evaluate the cost function for tracer concentrations only in the free troposphere in units of ppm.<br />
12: LCSPEC_PPB &amp; Evaluate the cost function for species concentrations in units of ppb, averaged over the range NSPAN. There are also hardwired options within CALC_ADJ_FORCE_FOR_SENS that can be used to specify a sub-domain over which to average: LMIN, LMAX, JMIN, JMAX, IMIN, IMAX.<br />
13: LPOP_UGM3 &amp; Domain-wide average population weighted aerosol concentrations<br />
14: spacer line &amp; This next section is only important if you are using a cost function that involves tracers (STT).<br />
15: NOBS &amp; The number of tracers involved in your cost function. It must match the number of tracers listed below, or if it is zero the section below will be ignored.<br />


l p12cm l

16: subheader &amp; -<br />
17: Tracer #1 &amp; List the ID number of the tracer from input.geos (TRC#) and its name (tracer_name). Make as many entries in this section as necessary.<br />
18- …: &amp; …<br />
19: spacer &amp; -<br />
20: NOBS_CSPEC &amp; The number of species involved in your cost function. It must match the number of species listed below, or if it is zero the section below will be ignored.<br />
21: subheader &amp; -<br />
22: Species #1 &amp; Enter the names of the species to be observed (species_name) They don’t need to be ordered or numbered according to their definition in CSPEC, just make sure that the name is exactly as it is listed in globchem.dat so that the code can match the name and find the corresponding an ID index in CSPEC. Make as many entries in this section as necessary.<br />
23-…: &amp; …<br />


<pre>01: %%% FINITE DIFFERENCE MENU %%%
02: fd perturbation        FD_DIFF : 0.1
03: Numerator of derivative to test:---
04:  =&gt; longitude degree   LONFD   : 32
05:  =&gt; latitude  degree   LATFD   : 21
06:  =&gt; OR pick box by grid index? : T
07:    =&gt; longidute index  IFD     : 21
08:    =&gt; latitude index   JFD     : 31
09:  =&gt; altidude index     LFD     : 1
10:  =&gt; tracer (STT TRC#)  NFD     : 2
11: Denomenator of deriv.   to test:
12:  =&gt; w/LEMS: emis group MFD     : 1
13:  =&gt; w/LEMS: sector     EMSFD   : 9
14:  =&gt; w/LICS: tracer     ICSFD   : 1
15:  =&gt; w/LSTR: tracer     STRFD   : 1</pre>
''Description'' Selections in this menu will apply if we are performing a finite difference test, i.e. LSENS = T and either FD_GLOB or FD_SPOT = T in the ADJOINT SIMULATION MENU.

l p12cm l

01: header line &amp; –<br />
02: ''FD_DIFF'' &amp; The size of the finite difference perturbation that is applied to the control parameters (LICS or LADJ_EMS), depending on which ones are selected for the finite difference test.<br />
03: spacer line &amp; The numerator of the derivative to test is selected in the following submenu. Note: for debugging, it can be useful to point these indices to troublesome values and then turn on LPRINTFD in the DIAGNOSTICS MENU.<br />
04: ''LONFD'' &amp; Longitude of the gridbox in the SPOT finite difference test and debugging (converted to index online).<br />
05: ''LATFD'' &amp; Latitude of the gridbox in the SPOT finite difference test and debugging (converted to index online).<br />
06: specify box &amp; This flag, if set to TRUE, will set the finite difference box as specified by the gridbox indecies, and not lat/lon indecies.<br />
07: ''IFD'' &amp; Longitude index of the gridbox in the SPOT finite difference test and debugging.<br />
08: ''JFD'' &amp; Latitude index of the gridbox in the SPOT finite difference test and debugging.<br />
09: ''LFD'' &amp; Vertical level index of the gridbox in the SPOT and GLOB finite difference test and debugging.<br />
10: ''NFD'' &amp; Tracer (STT TRC#) index of the SPOT and GLOB finite difference test and in debugging. This value will override whatever tracer is listed in the OBSERVATION MENU.<br />
11: spacer line &amp; The denomenator of the derivative to test (for both SPOT and GLOB) is selected in the following submenu.<br />
12: ''MFD'' &amp; Emission temporal group index of the finite difference test. Only matters for LADJ_EMS.<br />
13: ''EMSFD'' &amp; Emission sector index (full chem only) of the finite difference test. Only matters for LADJ_EMS.<br />
14: ''ICSFD'' &amp; Initial conditions type of the finite difference test. Only matters for LICS.<br />
15: ''STRFD'' &amp; Stratospheric tracer index (full chem only) of the finite difference test. Only matters for LADJ_STRAT. NOTE: Finite difference test for the stratospheric fluxes is available for one of either production or loss rates. By default it is for loss rates. You can change it to production rate by modifying the code (inverse_mod.f).<br />


<pre>01: %%% DIAGNOSTICS MENU %%%
02: General                        : T
03:  =&gt; print debug  LPRINTFD      : T
04:  =&gt; jsave, jsave2              : F
05:  =&gt; adjoint traj LADJ_TRAJ     : T
06:    =&gt; w.r.t. scale factors?    : T
07:  =&gt; save iteration diags LITR  : T
08:  =&gt; sense w.r.t absolute emis  : F
09: CO satellite diganostics       : F
10:  =&gt; H(model)                   : F
11:  =&gt; h(obs)                     : F
12:  =&gt; H(model)-h(obs)            : F
13:  =&gt; adjoint forcing            : F
14:  =&gt; model bias                 : F
15:  =&gt; observation count          : F
16:  =&gt; DOFs                       : F
17: TES NH3 diagnostics            :---
18:  =&gt; BLVMR                      : F</pre>
''Description''<br />


l p12cm l

01: header line &amp; –<br />
02: general switch &amp; Global diagnostic switch. Needs to be set to TRUE for any of the diagnostics to be saved or printed.<br />
03: ''LPRINTFD'' &amp; Print (to log file) debugging messages and tracer values in (IFD,JFD,LFD,NFD) gridbox.<br />
04: jsave, jsave2 &amp; Switch to save jsave and jsave2 values (debugging/finite difference output).<br />
05: ''LADJ_TRAJ'' &amp; Switch to store adjoint trajectory. Note: turning this on can generate a lot of output for fullchemistry simulations. You may want to edit the source code (routine MAKE_ADJ_FILE) so that only a subset of the adjoint species are written out, or a subset of vertical levels.<br />
06: ''LADJ_TRAJ'' &amp; Switch to save adjoint trajectories as sensitivities w.r.t scaling factors.<br />
07: ''LITR'' &amp; Save iteration diagnostic file ''gctm.iteration'' to the diagnostic directory. This file contains information related to convergence of the optimization routine.<br />
08: ''LEMS_ABS'' &amp; Switch to save sensitivities w.r.t. absolute values of emissions (as opposed to scaling factors). These will be in the ems.adj.* files in the diagnostic directory. Currently only supported for SO<math>_2</math>, NH<math>_3</math>, BC and OC emissions.<br />
09: CO diags &amp; Global switch for CO satellite diagnostics. To be used only when MOPITT, AIRS and/or SCIAMACHY observations are used. Corresponding code is very easy to adapt to other observational operators and saves information in the same space for model and observations, making model-data comparisons very easy.<br />
10: ''H(model)'' &amp; Switch to save model convolved by satellite averaging kernels, currently as a column value.<br />
11: ''h(obs)'' &amp; Switch to save satellite observations averaged on GEOS-Chem grid (1h resolution).<br />
12: ''H(model)-h(obs)'' &amp; Switch to save model-satellite differences in satellite/model space.<br />
13: adj forcing &amp; Switch to save adjoint forcing.<br />
14: model bias &amp; Switch to save model bias (model-obs)/obs. This is very useful for computing Relative Residual Error (RRE) values for optimization studies.<br />
15: ''OBS_COUNT'' &amp; Switch to save the number of observations per gridbox in a given simulation.<br />
16: ''DOF'' &amp; Switch to save an average degrees of freedom (DOF) per gridbox.<br />
17: TES NH<math>_3</math> diags &amp; switches for TES NH<math>_3</math> assimilation<br />
18: BLVMR &amp; Switch to read in BLVMR values from TES and add GEOS-Chem BLVMR values to the TES data files.<br />


=== ''define_adj.h'': observation selection ===

The following are observation options, all set in ''define_adj.h''. The corresponding files are in ''code/obs/''. Some are not yet fully implemented, as indicated by “'''Placeholder''' ”.

Note: No observation operators are needed for basic sensitivity runs. Basic sensitivity runs are those where the cost function is calculated directly within subroutine ''CALC_ADJ_FORCE_FOR_SENS''.

''CO''

l p12cm l

<tt>MOPITT_V3_CO_OBS</tt> &amp; Assimilate MOPITT CO (column) v3 observations to perform an optimization problem; data in hdf-eos4 file format.<br />
<tt>MOPITT_V4_CO_OBS</tt> &amp; Assimilate MOPITT CO (column) v3 observations to perform an optimization problem; data in hdf-eos4 file format.<br />
<tt>AIRS_CO_OBS</tt> &amp; Assimilate CO (column) observations from AIRS v5, data in hdf-eos4 file format.<br />
<tt>SCIA_BRE_CO_OBS</tt> &amp; Assimilate CO (column) observations from SCIAMACHY (Bremen retrieval only), data in ASCII format.<br />


''Aerosols''

l p12cm l

<tt>TES_NH3_OBS</tt> &amp;<br />
<tt>SCIA_DAL_SO2_OBS</tt> &amp; '''Placeholder''' for SO<math>_2</math> from SCIA.<br />
<tt>PM_ATTAINMENT</tt> &amp; '''Placeholder''' for aerosol attainment sensitivities; no observations (pseudo or real) required.<br />
<tt>IMPROVE_SO4_NIT_OBS</tt> &amp; '''Placeholder''' for aerosol observations from IMPROVE data files.<br />
<tt>CASTNET_NH4_OBS</tt> &amp; '''Placeholder''' for aerosol observations from CASTNET data files.<br />


''Ozone''

l p12cm l

<tt>SOMO35_ATTAINMENT</tt> &amp; '''Placeholder''' for ozone attainment sensitivities; no observations (pseudo or real) required.<br />
<tt>TES_O3_OBS</tt> &amp; TES O<math>_3</math> data in netcdf format from JPL.<br />


''NO2''

l p12cm l

<tt>SCIA_KNMI_NO2_OBS</tt> &amp; '''Placeholder''' for SCIAMACHY or GOME NO<math>_2</math> data from KNMI hdf files.<br />
<tt>SCIA_DAL_NO2_OBS</tt> &amp; '''Placeholder''' for …<br />


''other options''

l p12cm l

<tt>PSEUDO_OBS</tt> &amp; Use pseudo observations (generated by the model).<br />
<tt>LOG_OPT</tt> &amp; Use log-scaling factors.<br />
<tt>LIDORT</tt> &amp; Compile LIDORT code for radiative forcing caculations.<br />
<tt>LBFGS</tt> &amp; Calculate inverse Hessian matrix using L-BFGS method.<br />


=== Other input files ===

All the files listed below, if needed, should be place in ''run/../geos*/'' directory.

* ''run'' script, see Sec. [sec:run<sub>s</sub>cript].
* Observational error file(s): for each dataset used. Currently the code expects it in the following format:<br />
<tt>RRE_YYYYMMairsGlobal.bpch</tt><br />
<tt>RRE_seasonMay1mopittGlobal.bpch</tt><br />
<tt>RRE_seasonMay1sciabrGlobal.bpch</tt><br />
So for example for AIRS data, in May 2004, the file name is RRE_200405airsGlobal.bpch. AIRS data uses monthly errors, while MOPITT and SCIAMACHY (due to scarcity of data), use seasonal errors. Currently there is 1 file per month for AIRS and 1 file each for MOPITTT and SCIA for the whole year. Current setting spans May 1, 2004 to May 1, 2005 NOTE: There is an option to specify a fixed error across all times and gridboxes and not use the files. This is useful for saving model and data to compute Relative Residual Error (RRE) quantities for later use in an inversion.
* ''iter.txt'': This input file contains the iteration number of the next execution of the model. It starts with 1 and the code updates it. This way the same executable can be used for multiple runs (different dates and different iterations).
* MOPITT files: ''mopitt_v3_apriori.dat'' (a priori profile)
* SCIAMACHY files: ''ak_co_wfmdscia_V2.dat'' (averaging kernels, a priori profile), ''SCIA_pressure.dat'' (pressure levels for SCIA Bremen)

== Output files ==

See Sect. [sec:dir] for locations. All are binary bunch files unless otherwise noted. Note: files that don’t have an iteration token in their name will be overwritten or removed at each iteration. Exceptions are ''aero.ave*'' and ''satave*'' files, which are lumped into ''fwd_dat.*.tar'' files after each iteration, see the run script.

=== Essential output files ===

lp8 cml

''adjtmp/gctm*.chk.*''&amp; Checkpoint files. Generated during the forward run; deleted after they are used in the backward run (although the <tt>L_DEL_CHECKPT</tt> flag in ''CMN_ADJ'' allows you to not delete them if desired). See ''checkpt_mod.f'' for more details on content.<br />
''OptData/gctm.gdt.NN'' &amp; Gradients of active parameters at each iteration NN (<tt>IJ-GDE-$</tt>). These gradients are semi-normalized. To include fully-normalized gradients (<tt>IJ-GDEN$'</tt> for diagnostic purposes), set <tt>L_WRITE_GDEN = .TRUE.</tt> at the top of the rountine <tt>MAKE_GDT_FILE</tt>.<br />
''OptData/cfn.NN'' &amp; Cost function at each iteration NN. ASCII<br />


=== Nonessential output files ===

lp8 cm l

''OptData/gctm.ics.NN'' &amp; Scaling factors at each iteration NN. The scaled emissions themselves at each iteration are also included in this file for diagnostic purposes. Hence, the a priori estimates of emissions will be found in ''gctm.ics.01'', optimized emissions found in ''gctm.ics.10'', etc. The scaling factors are <tt>IJ-EMS-$</tt> and the emissions themselves are <tt>IJ-EM0-$</tt>.<br />
''OptData/gctm.obs.NN'' &amp; Satellite observations in the CO version. Files contain 3 data structures/arrays. Tracer 1 corresponds to MOPITT obs, 2 to SCIA, 3 is AIRS. The arrays are (IIPAR, JJPAR, NDAYS), where NDAYS is number of simulation days.<br />


''OptData/gctm.model.NN'' &amp; Model columns corresponding to the satellite data in time, space and retrieval processing to observations in gctm.obs.NN. These arrays correspond to the satellite ones and are updated at each iteration of the optimization.<br />


''OptData/gctm.costf.NN'' &amp; These files contain 2 arrays; 1: Cumulative cost function (summed over all observation types) with dimensions (IIPAR,JJPAR,NDAYS), 2: observation count in each gridbox with dimensions (IIPAR,JJPAR).<br />


''OptData/gctm.modbias.NN'' &amp; Contains <tt>(model-obs)/model</tt> at (IIPAR,JJPAR,NDAYS) resolution. Useful for computing Relative Residual Errors (RRE).<br />


''OptData/gctm.forcing.NN'' &amp; Contains <tt>2*(model-obs)/err^2</tt>, also known as adjoint forcing; a diagnostic.<br />


''OptData/gctm.gdta.NN'' &amp; Gradients of all parameters at each iteration NN. Only generated if <tt>CALL MAKE_GDT_ALL_FILE</tt> in ''inverse_driver.f''.<br />


''OptData/gctm.arr'' &amp; Integrated reaction rate constant sensitivities.<br />


lp8cm l

''diagadj/gctm.adj.YYYYMMDD.hhmmss'' &amp; Adjoint state variables (<math>\lambda_c</math>). Purely diagnostic; only written if <tt>LADJ_TRAJ = T</tt> (in ''input.gcadj'') and <tt>ITS_TIME_FOR_OBS</tt>. The number of these files kept can be controlled using the <tt>REMOVE_ADJ_FILE</tt> routine (defunct, need to fix!) in ''geos_chem_adj_mod.f''and the <tt>N_ADJ_KEEP</tt> parameter in ''input.gcadj''. Can be saved as sensitivities with respect to concentrations or concentration scaling factors, see <tt>MAKE_ADJ_FILE</tt> in ''geos_chem_adj_mod.f''.<br />


''adjtmp/gctm.obs.YYYYMMDD.hhmmss'' &amp; Pseudo observation file.<br />


''diagadj/gctm.fd.YYYYMMDD.hhmmss'' &amp; Diagnostic file. Used for process specific finite difference tests.<br />


''diagadj/gctm.fdglob.YYYYMMDD.hhmmss'' &amp; Diagnostic file. Used for process specific second order finite difference tests. Only generated if FD_GLOB is set to true.<br />


''diagadj/aero.ave.YYYYMMDD'' &amp; Aerosol data and corresponding model predictions. Generated when running with <tt>IMPROVE_OBS</tt> or <tt>PM_ATTAINMENT</tt> options.<br />


''diagadj/jsave.YYYYMMDD'' &amp; Contributions to cost function. Generated when running with <tt>IMPROVE_OBS</tt> or <tt>ATTAINMENT</tt> options.<br />


''diagadj/satave.bpch'' &amp; Satellite data and corresponding model predictions for <tt>NO2_SAT_OBS</tt>.<br />


''diagadj/ems.adj.*'' &amp; Emissions sensitivities on per-kg basis (as opposed to scaling factors). Currently only implemented for SO<math>_2</math>, NH<math>_3</math>, OC and BC emissions. Enabled with the ''LEMS_ABS'' switch.<br />


''diagadj/gctm.iteration'' &amp; ASCII. Statistics from the optimization procedure. Enabled with the ''LITR'' switch.<br />


''runs/../geos*/FWD_met'' &amp; Elements of the forward meteorology in the FD cell.<br />


''runs/../geos*/BACKWD_met'' &amp; Elements of the backward meteorology in FD cell. Should agree exactly with contents of ''FWD_met''.<br />


== The ''run'' script ==

[sec:run<sub>s</sub>cript] Each run of the adjoint can be done with the same executable. Thus the simplest run script is sufficient for 1 iteration. The number in ITER updates ''N_CALC_STOP'' variable in ''inverse_driver.f'' and determines the number of iterations in the next execution of the code.

'''NOTE:''' The ''run'' script uses bash like syntax, we expect it to work if you have bash on your system. If you experience problems please write to the development team in other to adress the problem.

To use, you have to set several variables in the script, following the instructions therein.

Change often (almost every run):

<pre>X
XSTOP
RNAME
TYPE
Makefile</pre>
Change rarely (only when migrating to a new filesystem):

<pre>IFORT_OPT
RECOMPILE
SAVE
DSAVE
ARCHIVE
DARCHIVE
DRUNDIR
DPACK</pre>
Unlikely to need to be changed from the defaults:

<pre>DCODE
DRUN</pre>
The frequently changed variables are <tt>X</tt>, <tt>XSTOP</tt>, <tt>RNAME</tt> and <tt>TYPE</tt>. Set <tt>RNAME</tt> to be a descriptive name of your current calculation (such as <tt>ADJv23_optimize_NH3_emissions</tt>). Set <tt>TYPE</tt> acording to the type of simulation operator, the model supports 4 types of simulation:

<pre>DEFAULT    : As it says is the default type that most of the people run.
HDF        : Set this type when satellite HDF data is going to be used.
LIDORT     : Set this type when running a LIDORT simulation.
SAT_NETCDF  : Set this type when satellite NetCDF data is goin to be used.</pre>
Set the start (or current) iteration number, <tt>X</tt>. Set <tt>XSTOP</tt> to the final iteration number. For example, if you’ve already computed five iterations, then set <tt>X=6</tt> and <tt>XSTOP=9</tt> to compute three more. At each iteration, the current value of <tt>X</tt> is assigned to the variable <tt>N_CALC_STOP</tt> in ''inverse_driver.f'' and the code is executed. <tt>X = 0</tt> is used for generating pseudo observations. For example, if <tt>X=0</tt> and <tt>XSTOP =4</tt>:

l | l | l | l | l

<br />
<tt>X=0</tt> &amp; <tt>X=1</tt> &amp; <tt>X=2</tt>&amp; <tt>X=3</tt> &amp; <tt>X=4</tt><br />
<math>\sigma</math>=1&amp;<math>\sigma</math>=1+<math>\delta \sigma</math> &amp; read *.01 &amp; read *.01 &amp;read *.01<br />
<tt>DO_GEOS_CHEM</tt>&amp;<tt>DO_GEOS_CHEM</tt> &amp; update <math>\sigma</math> &amp; update <math>\sigma</math> &amp; update <math>\sigma</math><br />
make *obs* &amp;<tt>DO_ADJOINT</tt> &amp; <tt>DO_GEOS_CHEM</tt> &amp; read *.02 &amp; read *.02<br />
&amp;write *.01 &amp; <tt>DO_ADJOINT</tt> &amp; update <math>\sigma</math> &amp; update <math>\sigma</math><br />
&amp; &amp; write *.02 &amp; <tt>DO_GEOS_CHEM</tt> &amp; read *.03<br />
&amp; &amp; &amp; <tt>DO_ADJOINT</tt> &amp; update <math>\sigma</math><br />
&amp; &amp; &amp; write *.03 &amp; <tt>DO_GEOS_CHEM</tt><br />
&amp; &amp; &amp; &amp; <tt>DO_ADJOINT</tt><br />
&amp; &amp;&amp;&amp; write *.04<br />


A reason to do only a single iteration per execution is because of code fragments like this:

<tt>LOGICAL</tt>, <tt>SAVE   :: FIRST = .TRUE.</tt>

<pre>IF ( FIRST ) THEN
   CALL INIT_ARRAYS
   FIRST = .FALSE.
ENDIF</pre>
By exiting the code between each iteration, these logical control switches get automatically reset to <tt>FIRST = .TRUE.</tt>. Doesn’t this repeated reading / writing / optimizing take a lot of time? No, not with respect to the expense of the forward and backward model calculations.

For sensitivity calculations, set <tt>X</tt> and <tt>XSTOP</tt> to 1

For global finite difference calculations, set <tt>X=1</tt> and <tt>XSTOP=3</tt>.

For spot finite difference calculations, set <tt>X=1</tt> and <tt>XSTOP=2</tt>.

For generating pseudo observations using ''PSEUDO_OBS'', set <tt>X=0</tt> and <tt>XSTOP</tt> equal to the total number of iterations you wish to perform. Subsequent tests using the same set of pseudo observations can start from <tt>X=1</tt>.

The run script will compile the ''geos'' executable,

<pre> ##################################################################
 # Compile geos, move it to the run directory and execute
 ##################################################################
        
 cd $DRUN/$DPACK/$DCODE
 if  [ -n $(grep 'Compute BFGS inverse Hessian: T' $DRUNDIR/input.gcadj) ]; then
        ./objects.sh $TYPE
 else
        ./objects.sh $TYPE LBFGS
  fi
        
  if   [ $TYPE = 'DEFAULT' ]; then
       IFORT_OPT=&quot;$IFORT_OPT &quot;
  elif [ $TYPE = 'HDF' ]; then
        IFORT_OPT=&quot;$IFORT_OPT HDF=yes&quot; 
  elif [ $TYPE = 'SAT_NETCDF' ]; then
        IFORT_OPT=&quot;$IFORT_OPT SAT_NETCDF=yes&quot;
  elif [ $TYPE = 'LIDORT' ]; then
        IFORT_OPT=&quot;$IFORT_OPT LIDORT=yes&quot;
  fi
                
  if  [ -n $(grep 'Compute BFGS inverse Hessian: T' $DRUNDIR/input.gcadj) ]; then
      make $IFORT_OPT 
 else
      make LBFGS='yes' $IFORT_OPT
 fi
        
 if [ $RECOMPILE = 'YES' ]; then    
     mv -f geos ../$DRUNDIR/
 else
      cp -f geos ../$DRUNDIR/
 fi
        
 cd ../$DRUNDIR/
 time ./geos</pre>
If you are running an interactive job and wish to save a log file that doesn’t get overwritten each iteration, than you may replace

<pre>        time ./geos</pre>
with

<pre>        time ./geos &gt; log.${X}</pre>
The run script will optionally save a copy of the entire package and move it to <tt>DSAVE</tt>. This may be necessary if you are using <tt>DRUN</tt> on the local storage of a compute node to which you no longer have access after your log has completed. The enable or disable this feature, set the variable <tt>SAVE</tt>.

Lastly, users running parallel jobs with OpenMP on multi-core compute nodes should be mindful to set the following to be appropriate for their system:<br />


<pre>#  Set number of threads
export OMP_NUM_THREADS=24</pre>
= 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

<tt>resolution:</tt> &amp; 4x5 and 2x2.5, nested CO capabilities<br />
<tt>simulation types:</tt> &amp; full chemistry, tagged CO, tagged O3 and offline CO<math>_2</math>.<br />
<tt>meteorology:</tt> &amp; GEOS3, GEOS4, GEOS5<br />
<tt>chemistry:</tt> &amp; via KPP Rosenbrock solver<br />


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 <tt>'make help'</tt> inside the <tt>code/</tt> 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:fd<sub>c</sub>hecklist] 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:4dvar<sub>c</sub>hecklist] 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 <tt>LOG_OPT</tt>) 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, <tt>NCOEFF</tt> in ''gckpp_adj_Global.f90''. <tt>NCOEFF</tt> is a the size of the vector <tt>JCOEFF</tt>, 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 <tt>JCOEFF</tt>. Define the value of JCOEFF(new index) to be the index of the reaction rate in the KPP generated reaction rate constant vector <tt>RCONST</tt>. 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 <tt>RCONST2RKPP</tt> 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 <tt>ros_DadjInt</tt> 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 (Henze, Hakami, and Seinfeld 2007). 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> (<tt>DJDR</tt>) and <math>f_p^T(T_i,C_i)</math> (<tt>DFDR</tt>) 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'':

<pre>{40}  ALD2 + NO3 = HNO3 + MCO3:                 RKPP(JLOOP,54);</pre>
If <tt>NCOEFF</tt> was previously 8, then we would set <tt>NCOEFF = 9</tt> in ''gckpp_adj_Global.f90''. In ''gckpp_adj_Util.f90'', in the subroutine <tt>INIT_KPP</tt> we would add the following lines:

<tt>      JCOEFF(9) = 40   ! ALD2 + NO3 --&gt; HNO3 + MCO3</tt>

<tt>      RCONST2RRATE(40) = 54</tt>

As mentioned above, in the routine <tt>ros_DadjInt</tt> in ''gckpp_adj_Integrator.f90'', these sensitivities can be calculated using the general formula, or one could implement the sparse formula:

<pre> 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)</pre>
==== Saving reaction rate sensitivities ====

The subroutine <tt>SAVE_INST_ARR</tt> can be used to print out reaction rate sensitivities in a particular cell at every time step. This routine also calls the routine <tt>MAKE_ARR_FILE</tt> 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 <tt>ARE_FLAGS_VALID</tt> 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 (Henze, Hakami, and Seinfeld 2007).

For example, to check the adjoint of the chemistry only in a single vertical level, set the following:

in ''run'':

<tt>X=1</tt>

<tt>XSTOP=3</tt>

<tt>TYPE=DEFAULT</tt>

in ''input.gcadj'':

<tt>LSENS</tt> = T

<tt>FD_GLOB</tt> = T

<tt>LICS</tt> or <tt>LADJ_EMS</tt> or <tt>LADJ_EMS</tt> and <tt>LADJ_STRAT</tt> = T

select a dependent species, <tt> NFD</tt>

select a vertical level, <tt> LFD</tt>

select a control parameter, <tt>ICSFD</tt> or <tt>EMSFD</tt> and/or <tt>STRFD</tt>

<tt>IFD</tt> and <tt>JFD</tt> don’t matter as we’re doing a domain wide test

set the finite difference perturbation (<math>\delta\sigma</math>), <tt> FD_DIFF = 1.d-1</tt>

Set the <tt>OBS_FREQ</tt> to 60, but set LMAX_OBS = T and NSPAN = 1.

Set <tt>LAERO_THERM = .FALSE.</tt>

in ''input.geos'':

Turn off convection ( <tt>LCONV = F</tt> )

Turn off turbulent mixing ( <tt>LTURB = F</tt> )

Turn off wet deposition ( <tt>LWETD = F</tt> )

Leave on dry deposition ( <tt>LDRYD = T</tt> )

Can leave on transport ( <tt>LTRAN = T</tt> ), which will be overridden if <tt>FD_GLOB = T</tt>.

in ''chemistry_mod.f'', can save time by computing gas-phase chemistry only in certain cells by placing an <tt>IF</tt> statement around <tt>CALL INTEGRATE_ADJ</tt> such as:

<pre>       IF ( L == LFD ) THEN

           CALL INTEGRATE_ADJ(...

       ENDIF</pre>
Analysis of global FD tests is described in Sect. [sec:geos5<sub>b</sub>ench].

== 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'':
** <tt>X=1</tt>
** <tt>XSTOP=2</tt>
** <tt>TYPE=DEFAULT</tt>
* in ''input.gcadj'':
** <tt>LSENS</tt>
** Make sure that <tt>FD_GLOB</tt> is set to FALSE and that <tt> FD_SPOT</tt> is set to TRUE.
** active variables are <tt>LICS</tt> or <tt>LADJ_EMS</tt> or <tt>LADJ_EMS</tt> and <tt>LADJ_STRAT</tt>
** select a dependent species, <tt> NFD</tt>
** select particular control parameter <tt>IFD</tt> <tt>JFD</tt> <tt>LFD</tt> and <tt>ICSFD</tt> or <tt>EMSFD</tt> and/or <tt>STRFD</tt>
** set the finite difference perturbation (<math>\delta\sigma</math>), <tt> FD_DIFF = 1.d-1</tt>
** probably want to turn on additional diagnostic output, so set <tt>L_PRINTFD = .TRUE.</tt>
** Make sure that <tt>OBS_FREQ</tt> 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 <tt>ADJ / FD</tt> 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 <tt>ADJ / FD</tt> can not be expected to be unity. Benchmark tests are given in (Henze, Hakami, and Seinfeld 2007).

= Generating forward and reverse code for GEOS-Chem with the Kinetic PreProcessor (KPP) =

(Damian et al. 2002; Sandu, Daescu, and Carmichael 2003; [CSL BIBLIOGRAPHIC DATA ERROR: reference &quot;Adj1.b&quot; not found.])

== KPP input files ==

== Post processing KPP generated code ==

=== Interfacing KPP code with GEOS-Chem ===

=== Implementing OpenMP Parallelization in KPP generated code ===

== Performing global benchmarks of new chemical solvers ==

Anon. “Adj1.b not found!”

Damian, V., A. Sandu, M. Damian, F. Potra, and G. R. Carmichael. 2002. “The kinetic preprocessor KPP - a software environment for solving chemical kinetics.” ''Comput. Chem. Eng.'' 26: 1567–1579.

Henze, D. K., A. Hakami, and J. H. Seinfeld. 2007. “Development of the adjoint of GEOS-Chem.” ''Atmos. Chem. Phys.'' 7: 2413–2433.

Henze, D. K., J. H. Seinfeld, and D. Shindell. 2009. “Inverse modeling and mapping U.S. air quality influences of inorganic PM$_2.5$ precursor emissions using the adjoint of GEOS-Chem.” ''Atmos. Chem. Phys.'' 9: 5877–5903.

Kopacz, M., D. J. Jacob, J. A. Fisher, J. A. Logan, L. Zhang, I. A. Megretskaia, B. M. Yantosca, et al. 2009. “Global estimates of CO sources with high resolution by adjoint inversion of multiple satellite datasets (MOPITT, AIRS, SCIAMACHY, TES) .” ''Atmos. Chem. Phys. Discuss.'' submitted.

Kopacz, M., D. Jacob, D. K. Henze, C. L. Heald, David G. Streets, and Qiang Zhang. 2009. “A comparison of analytical and adjoint Bayesian inversion methods for constraining Asian sources of CO using satellite (MOPITT) measurements of CO columns.” ''J. Geophys. Res.-Atmos.'' 114. doi:0.1029/2007JD009264.

Sandu, A., D. N. Daescu, and G. R. Carmichael. 2003. “Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part I - theory and software tools.” ''Atmos. Environ.'' 37: 5083–5096.

Singh, K., P. Eller, A. Sandu, D. K. Henze, K. Bowman, M. Kopacz, and M. Lee. 2009. “Towards the construction of a standard GEOS-Chem adjoint model.” ''ACM High Performance Computing Conference''.

Zhang, L., D. J. Jacob, M. Kopacz, D. K. Henze, K. Singh, and D. A. Jaffe. 2009. “Intercontinental source attribution of ozone pollution at western US sites using an adjoint method.” ''Geophys. Res. Lett.'' 36. doi:10.1029/2009gl037950.

= 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, <tt>NCOEFF</tt> in ''gckpp_adj_Global.f90''. <tt>NCOEFF</tt> is a the size of the vector <tt>JCOEFF</tt>, 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 <tt>JCOEFF</tt>. Define the value of JCOEFF(new index) to be the index of the reaction rate in the KPP generated reaction rate constant vector <tt>RCONST</tt>. 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 <tt>RCONST2RKPP</tt> 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 <tt>ros_DadjInt</tt> 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> (<tt>DJDR</tt>) and <math>f_p^T(T_i,C_i)</math> (<tt>DFDR</tt>) 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'':

<pre>{40}  ALD2 + NO3 = HNO3 + MCO3:                 RKPP(JLOOP,54);</pre>
If <tt>NCOEFF</tt> was previously 8, then we would set <tt>NCOEFF = 9</tt> in ''gckpp_adj_Global.f90''. In ''gckpp_adj_Util.f90'', in the subroutine <tt>INIT_KPP</tt> we would add the following lines:

<tt>      JCOEFF(9) = 40   ! ALD2 + NO3 --&gt; HNO3 + MCO3</tt>

<tt>      RCONST2RRATE(40) = 54</tt>

As mentioned above, in the routine <tt>ros_DadjInt</tt> in ''gckpp_adj_Integrator.f90'', these sensitivities can be calculated using the general formula, or one could implement the sparse formula:

<pre> 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)</pre>
==== Saving reaction rate sensitivities ====

The subroutine <tt>SAVE_INST_ARR</tt> can be used to print out reaction rate sensitivities in a particular cell at every time step. This routine also calls the routine <tt>MAKE_ARR_FILE</tt> 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 <tt>ARE_FLAGS_VALID</tt> 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'':

<tt>X=1</tt>

<tt>XSTOP=3</tt>

<tt>TYPE=DEFAULT</tt>

in ''input.gcadj'':

<tt>LSENS</tt> = T

<tt>FD_GLOB</tt> = T

<tt>LICS</tt> or <tt>LADJ_EMS</tt> or <tt>LADJ_EMS</tt> and <tt>LADJ_STRAT</tt> = T

select a dependent species, <tt> NFD</tt>

select a vertical level, <tt> LFD</tt>

select a control parameter, <tt>ICSFD</tt> or <tt>EMSFD</tt> and/or <tt>STRFD</tt>

<tt>IFD</tt> and <tt>JFD</tt> don’t matter as we’re doing a domain wide test

set the finite difference perturbation (<math>\delta\sigma</math>), <tt> FD_DIFF = 1.d-1</tt>

Set the <tt>OBS_FREQ</tt> to 60, but set LMAX_OBS = T and NSPAN = 1.

Set <tt>LAERO_THERM = .FALSE.</tt>

in ''input.geos'':

Turn off convection ( <tt>LCONV = F</tt> )

Turn off turbulent mixing ( <tt>LTURB = F</tt> )

Turn off wet deposition ( <tt>LWETD = F</tt> )

Leave on dry deposition ( <tt>LDRYD = T</tt> )

Can leave on transport ( <tt>LTRAN = T</tt> ), which will be overridden if <tt>FD_GLOB = T</tt>.

in ''chemistry_mod.f'', can save time by computing gas-phase chemistry only in certain cells by placing an <tt>IF</tt> statement around <tt>CALL INTEGRATE_ADJ</tt> such as:

<pre>       IF ( L == LFD ) THEN

           CALL INTEGRATE_ADJ(...

       ENDIF</pre>
Analysis of global FD tests is described in Sect. [sec:geos5<sub>b</sub>ench].

== 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'':
** <tt>X=1</tt>
** <tt>XSTOP=2</tt>
** <tt>TYPE=DEFAULT</tt>
* in ''input.gcadj'':
** <tt>LSENS</tt>
** Make sure that <tt>FD_GLOB</tt> is set to FALSE and that <tt> FD_SPOT</tt> is set to TRUE.
** active variables are <tt>LICS</tt> or <tt>LADJ_EMS</tt> or <tt>LADJ_EMS</tt> and <tt>LADJ_STRAT</tt>
** select a dependent species, <tt> NFD</tt>
** select particular control parameter <tt>IFD</tt> <tt>JFD</tt> <tt>LFD</tt> and <tt>ICSFD</tt> or <tt>EMSFD</tt> and/or <tt>STRFD</tt>
** set the finite difference perturbation (<math>\delta\sigma</math>), <tt> FD_DIFF = 1.d-1</tt>
** probably want to turn on additional diagnostic output, so set <tt>L_PRINTFD = .TRUE.</tt>
** Make sure that <tt>OBS_FREQ</tt> 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 <tt>ADJ / FD</tt> 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 <tt>ADJ / FD</tt> can not be expected to be unity. Benchmark tests are given in .

= Generating forward and reverse code for GEOS-Chem with the Kinetic PreProcessor (KPP) =


== KPP input files ==

== Post processing KPP generated code ==

=== Interfacing KPP code with GEOS-Chem ===

=== Implementing OpenMP Parallelization in KPP generated code ===

== Performing global benchmarks of new chemical solvers ==


=References=

*Daescu, D. N., A. Sandu, and G. R. Carmichael (2003), Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: II - numerical validation and applications, Atmos. Environ., 37(36), 5097–5114.

*Damian, V., A. Sandu, M. Damian, F. Potra, and G. R. Carmichael. 2002. “The kinetic preprocessor KPP - a software environment for solving chemical kinetics.” ''Comput. Chem. Eng.'' 26: 1567–1579.

*Henze, D. K., A. Hakami, and J. H. Seinfeld. 2007. “Development of the adjoint of GEOS-Chem.” ''Atmos. Chem. Phys.'' 7: 2413–2433.

*Henze, D. K., J. H. Seinfeld, and D. Shindell. 2009. “Inverse modeling and mapping U.S. air quality influences of inorganic PM$_2.5$ precursor emissions using the adjoint of GEOS-Chem.” ''Atmos. Chem. Phys.'' 9: 5877–5903.

*Kopacz, M., D. J. Jacob, J. A. Fisher, J. A. Logan, L. Zhang, I. A. Megretskaia, B. M. Yantosca, et al. 2009. “Global estimates of CO sources with high resolution by adjoint inversion of multiple satellite datasets (MOPITT, AIRS, SCIAMACHY, TES) .” ''Atmos. Chem. Phys. Discuss.'' submitted.

*Kopacz, M., D. Jacob, D. K. Henze, C. L. Heald, David G. Streets, and Qiang Zhang. 2009. “A comparison of analytical and adjoint Bayesian inversion methods for constraining Asian sources of CO using satellite (MOPITT) measurements of CO columns.” ''J. Geophys. Res.-Atmos.'' 114. doi:0.1029/2007JD009264.

*Sandu, A., D. N. Daescu, and G. R. Carmichael. 2003. “Direct and adjoint sensitivity analysis of chemical kinetic systems with KPP: Part I - theory and software tools.” ''Atmos. Environ.'' 37: 5083–5096.

*Singh, K., P. Eller, A. Sandu, D. K. Henze, K. Bowman, M. Kopacz, and M. Lee. 2009. “Towards the construction of a standard GEOS-Chem adjoint model.” ''ACM High Performance Computing Conference''.

*Zhang, L., D. J. Jacob, M. Kopacz, D. K. Henze, K. Singh, and D. A. Jaffe. 2009. “Intercontinental source attribution of ozone pollution at western US sites using an adjoint method.” ''Geophys. Res. Lett.'' 36. doi:10.1029/2009gl037950.