GEOS-Chem Adjoint v32

From Geos-chem
Revision as of 19:49, 18 January 2012 by Daven (Talk | contribs) (Cleanup and merge CALC_APRIORI (adj32_017))

Jump to: navigation, search



Under development. Expected release date 02/01/2012

What's new in this version

Updates to forward model

Update tagged CO simulation (adj32_017)

tagged_co_mod.f has been updated to distribute emissions through the boundary layer and removing obsolete COSF, TIME_STEPS, GET_CO_CH4, SET_CO_CH4, and GET_OPT_CO.

Update several files to v9 code to support full-chem nested simulation (adj32_015)

Updates to the following files were made by Zhe Jiang to support the NESTED_NA and NESTED_SD simulation:

Update DEAD dust emissions for nested simulations and MERRA (v9-01-01)

When the forward model was updated to support MERRA, the dust emissions got retuned MERRA_implementation_details#Dust_emissions. But also the emissions for GEOS-5 were affected, and support was added for nested grid.

The dust_dead_mod.f file has been updated to match that in v9-01-01 (dust_dead_mod.f,v 1.3 2010/02/02 16:57:53 bmy). Common blocks CMN_SIZE and CMN_GCTM are reinstated for backwards compatibility.

--Daven 19:58, 6 January 2012 (EST)

Updated TOMS O3

toms_mod.f was updated to v9-01-02. This uses the newer TOMS inventory TOMS_200906 and allows for rescaling of years up to 2008 (previous version stopped at 2007).


cac_mod.f was updated to v9-01-02.


emissions_mod.f was updated to support GRID05x0666 for CAC, NEI05 and EMEP emissions following v9-01-02.


file_mod.f was updated to include IU_BC_05x06.


CMN_SIZE was updated to include size of NA and SD nested domain.


commsoil.h was updated to include size of NA and SD nested domain.


define.h was updated to include the NESTED_SD option.


tpcore_bc_mod.f has been updated to include

  • OPEN_BC_FILE_05x0666
  • ITS_TIME_FOR_BC_05x0666INIT_TPCORE_BC_05x0666

and the BC_05x06 array.


regrid_1x1_mod has been updated to support the NESTED_SD domain.

Update lightning_nox_mod.f

See GEOS-Chem_Adjoint_v32#Updated_lightning_parameterization_and_fix_for_cloud-top-height_algorithm_.28v9-01-01.29


acetone_mod.f has been updated to support the nested domains. It now matches v9-01-02, with the additional update for the adjoint that we do not multiply and divide by ACETONE, thus making OCEAN_SINK_ACET self-adjoint.


input_mod.f has been updated in several places to support nested domains.


transport_mod.f has been updated to support the NESTED_SD domain.


co2_mod.f has been updated to support nested domains.

Bug fixes in forward model

Ionic strength below zero in rpmares (adj32_001)

On prospero, we often get the message: Ionic strength below zero...negative concentrations.

Despite the setting of negative roots to zero, sometimes the code still ends up with negative CRUTES. Trap for it here to prevent mysterious crashing later. Also, add the following

          !IF ( CRUTES(1) .ne. ABS(CRUTES(1) ) ) CRUTES(1)=1d9
          !IF ( CRUTES(2) .ne. ABS(CRUTES(2) ) ) CRUTES(2)=1d9
          !IF ( CRUTES(3) .ne. ABS(CRUTES(3) ) ) CRUTES(3)=1d9

which for some reason (I have no idea why) prevents the code from returning a negative CRUTES on prospero.

--Daven 14:58, 6 January 2012 (EST)

Updated lightning parameterization and fix for cloud-top-height algorithm (v9-01-01)

Updates to a6_read_mod.f and lightning_nox_mod.f. Details in Lee's description:

--Daven 15:18, 6 January 2012 (EST)

Fix for GEIA emissions scaling factor over Botswana (v9-01-02h)

Updates to scale_anthro_mod.f. See GEOS-Chem_v9-01-02#Fix_for_GEIA_emissions_scaling_factor_over_Botswana

--Daven 16:16, 6 January 2012 (EST)

TINY parameter in convection_mod.f (v8-03-02)

Update to convection_mod.f GEOS-Chem_v8-03-02#TINY_parameter_in_convection_mod.f

--Daven 16:16, 6 January 2012 (EST)

Updated to REGRID_05x0666_NESTED for non-standard domains (v?-?-?)

Zhe Jiang wrote

I revised the SUBROUTINE REGRID_05x0666_NESTED.
The problem in the original code is that the Longitude of LL (Left-Low of nested window) has to be even number and the latitude of LL has to be integer number. It has no influence on our standard NA/CH/EU nested setting.
I think it would be safer to fix this problem because we perhaps need to define new nested window in future.

This subroutine has now been replaced with Zhe's update. Not sure if there is a corresponding update in the v9 forward model yet.

--Daven 19:06, 6 January 2012 (EST)

Fix variable declarations in UPBDFLX_NOY (adj32_008)

Thomas Walker noted that having DTDYN, AIRDENS, PNOY declared as REAL*4 can lead to overflow in the adjoint model. The upbdflx_mod.f has now been updated to use REAL*8 instead.

--Daven 18:00, 8 January 2012 (EST)

Updates to time_mod.f for tagged-CO (ad33_009)

Zhe Jiang wrote

I currently found it perhaps can only affect tagged-CO inversion. Thus, it can not be said as a BUG.
For example, if we do our inversion for June 1-15, when the backward running, because June 15 is not the last day of the calendar month, these emissions will not be read. However, we don't need to read them since the model will use the value in the forward run. For the tagged-CO backward run, the CH4 concentration will only be read if June 15 is the last day of the calendar year. Thus, in the original code, the CH4 concentration, and thus gradient on CH4 oxidation, is ZERO.
      IF ( ITS_A_NEW_YEAR() ) THEN
         CALL GET_GLOBAL_CH4( GET_YEAR(), .TRUE., A3090S, A0030S, A0030N, A3090N )
In my previous code, I saved the CH4 concentration in the forward run and then read them in the backward run. After I modified time_mod.f, I have removed the code for the CH4 saving/reading. Anyway, the effect is trivial. I can solve the problem in the tagged_CO inversion by modifying the previous code.

To address these issues, subroutines ITS_A_NEW_MONTH and ITS_A_NEW_YEAR have been updated in time_mod.f.

--Daven 18:18, 8 January 2012 (EST)

NEI05 ID names for non-full chemistry simulation (adj32_015)

Explanation from Zhe...

 IF ( .not. ITS_A_FULLCHEM_SIM() )  THEN    !Zhe (16/3/2011)
      ! Zhe (16/3/2011), (dkh,  
     $                     IDTACET, IDTALK4, IDTC2H6, IDTC3H8,
     $                     IDTOCPI, IDTBCPI,
     $                     IDTALD2, IDTCH2O, IDTPRPE, IDTMEK
     $     /)   

         IDTNOX  = 1
         IDTCO   = 4
         IDTSO2  = 26
         IDTSO4  = 27
         IDTNH3  = 30
         IDTACET = 9
         IDTALK4 = 5
         IDTC2H6 = 21
         IDTC3H8 = 19
         IDTOCPI = 35
         IDTBCPI = 34
         IDTALD2 = 11
         IDTCH2O = 20
         IDTPRPE = 18
         IDTMEK  = 10

--Daven 19:46, 16 January 2012 (EST)

Fixes for minor issues affecting nested-grid simulations (v9_01_02e)

Updates to nei2005_anthro_mod.f and input_mod.f, see GEOS-Chem_v9-01-02#Fixes_for_minor_issues_affecting_nested-grid_simulations

--Daven 19:51, 16 January 2012 (EST)

Updates to adjoint model

Updates to aerosol direct radiative forcing adjoint

Several updates in lidort_mod.f and mie_mod.f:

  • From Farhan Akhtar, now include primary OC in the LIDORT radiative forcing calculation.
  • Now switch back to using analytic derivatives from mie table rather than FD

--Daven 18:50, 8 January 2012 (EST)

Dust adjoint and MODIS AOD observation operator (adj32_011)

From Richard Xiaoguang Xu at UNL, the dust adjoint includes sensitivity with respect to dust emissions from the DEAD module. Also included in this update is an observation operator for assimilating AOD observations from MODIS using the UNL retrieval.

--Daven 00:23, 13 January 2012 (EST)

Inverse Hessian approximation (adj32_012)

The inverse Hessian approximation (e.g., Henze et al., 2009) using the DFP formula is now implemented with the LINVH flag in input.gcadj. Instructions for how to perform this calculation are included in the UPDATE_HESSIAN subroutine in the new inv_hessian_mod.f file.

--Daven 00:23, 13 January 2012 (EST)

Offline aerosol adjoint. (adj32_013)

Yuhao Mao and Qinbin Li at UCLA have implemented the offline aerosol simulation adjoint for black carbon. This update also included an observation operator for measurements of these species at IMPROVE sites, IMPROVE_BC_OC_OBS and improve_bc_mod.f.

The offline aerosol adjoint will also work for the recently implemented dust aerosol adjoint (it does not yet include support for sulfate aerosol, sea salt or SOA). For example, you may now run with input.geos

%%% TRACER MENU %%%     : 
Type of simulation      : 10
Number of Tracers       : 4
Tracer Entries -------> : TR#   Name  g/mole   Tracer Members; () = emitted
Tracer #1               :  1    BCPI    12.0
Tracer #2               :  2    OCPI    12.0
Tracer #3               :  3    BCPO    12.0
Tracer #4               :  4    OCPO    12.0

and input.gcadj

          FOR LADJ_EMS         :
NNEMS: ems groups implemented  : 12
Emission entries ------------> : EMS# ems_name       opt  SF_DEFAULT REG_PARAM ERROR
Emission #1                    : 1    IDADJ_EBCPI_an T    1          1         1
Emission #2                    : 2    IDADJ_EBCPO_an T    1          1         1
Emission #3                    : 3    IDADJ_EOCPI_an T    1          1         1
Emission #4                    : 4    IDADJ_EOCPO_an T    1          1         1
Emission #5                    : 5    IDADJ_EBCPI_bf T    1          1         1
Emission #6                    : 6    IDADJ_EBCPO_bf T    1          1         1
Emission #7                    : 7    IDADJ_EOCPI_bf T    1          1         1
Emission #8                    : 8    IDADJ_EOCPO_bf T    1          1         1
Emission #9                    : 9    IDADJ_EBCPI_bb T    1          1         1
Emission #10                   : 10   IDADJ_EBCPO_bb T    1          1         1
Emission #11                   : 11   IDADJ_EOCPI_bb T    1          1         1
Emission #12                   : 12   IDADJ_EOCPO_bb T    1          1         1

--Daven 19:27, 13 January 2012 (EST)

Augmented list of supported VOC sensitivities (adj32_014)

From Kateryna Lapina, sensitivities of anthropogenic, biofuel and biomass burning emissions of all VOCs (ALK4, MEK, ALD2, PRPE, C3H8, CH2O C2H6). To implement, simply list these in input.gcadj as:

Emission #34                   : 34   IDADJ_EALK4_an T    1          1         1
Emission #35                   : 35   IDADJ_EALK4_bf T    1          1         1
Emission #36                   : 36   IDADJ_EALK4_bb T    1          1         1
Emission #37                   : 37   IDADJ_EACET_an T    1          1         1
Emission #38                   : 38   IDADJ_EACET_bf T    1          1         1
Emission #39                   : 39   IDADJ_EACET_bb T    1          1         1
Emission #40                   : 40   IDADJ_EMEK_an  T    1          1         1
Emission #41                   : 41   IDADJ_EMEK_bf  T    1          1         1
Emission #42                   : 42   IDADJ_EMEK_bb  T    1          1         1
Emission #43                   : 43   IDADJ_EALD2_an T    1          1         1
Emission #44                   : 44   IDADJ_EALD2_bf T    1          1         1
Emission #45                   : 45   IDADJ_EALD2_bb T    1          1         1
Emission #46                   : 46   IDADJ_EPRPE_an T    1          1         1
Emission #47                   : 47   IDADJ_EPRPE_bf T    1          1         1
Emission #48                   : 48   IDADJ_EPRPE_bb T    1          1         1
Emission #49                   : 49   IDADJ_EC3H8_an T    1          1         1
Emission #50                   : 50   IDADJ_EC3H8_bf T    1          1         1
Emission #51                   : 51   IDADJ_EC3H8_bb T    1          1         1
Emission #52                   : 52   IDADJ_ECH2O_an T    1          1         1
Emission #53                   : 53   IDADJ_ECH2O_bf T    1          1         1
Emission #54                   : 54   IDADJ_ECH2O_bb T    1          1         1
Emission #55                   : 55   IDADJ_EC2H6_an T    1          1         1
Emission #56                   : 56   IDADJ_EC2H6_bf T    1          1         1
Emission #57                   : 57   IDADJ_EC2H6_bb T    1          1         1

where in the example above emissions 1-33 were the previously defined standard set.

--Daven 19:27, 13 January 2012 (EST)

Remove bin_data.f and adj_bin_data.f

From Zhe Jiang, remove the obsolete files bin_data.f and adj_bin_data.f, and remove their listings in Makefile.ifort.hdf

--Daven 20:27, 16 January 2012 (EST)

Several updates to support nested full chemistry adjoint (adj32_015)

Several updates from Zhe Jiang support running the adjoint full-chemistry simulation over the NESTED_NA and NESTED_SD domain. The latter is a new domain, like NESTED_NA, but slightly less padding around NA.

Support nested in emissions_adj.f

emissions_adj.f has been updated to support nested emissions from the CAC, NEI05 and EMEP inventories.

Cleanup code in adj_arrays_mod.f

Remove obsolete variable NOR in adj_arrays_mod.f

MOPITT v5 observation operator (adj32_016)

Zhe Jiang added support for the MOPITT v5 products.

  • Add MOPITT_V5_CO_OBS pre-processor flag to define_adj.h
  • Updates to adj_arrays_mod.f, mopitt_obs_mod.f , HdfVdModule.f90, input_adj_mod.f.

--Daven 01:45, 17 January 2012 (EST)

Updates for tagged CO simulation (adj32_017)

Include CH4 gradient

Update inverse_mod.f so that oxidation of CH4 is incorporated in the gradient

Cleanup and small fixes to tagged_co_adj_mod.f

Remove COSF and TIME_STEPS in tagged_co_adj_mod.f.

--Daven 01:45, 17 January 2012 (EST)

Cleanup and merge CALC_APRIORI

Updates from Zhe Jiang for calculation of the penalty term:

  • Now valid for full chemistry and tagged CO simulations
  • Add subroutine READ_APERROR

Enforce factor of 0.5 in cost function definition (adj32_017)

A consequence of the merging of the full-chem and CO routines in CALC_APRIORI is that the definition of the cost function needs to be made consistent. Originally, the CO code omitted the factors of 1/2 in both terms of the cost function. To correct for this, the cost function (and adjoint forcing) has been updated in the following CO observation operators to match the definition in the full-chem simulation:

  • airs_co_obs_mod.f
  • mopitt_obs_mod.f
  • sciabr_co_obs_mod

--Daven 14:39, 18 January 2012 (EST)

Bug fixes in adjoint model

Mover error checking for O3 assimilation (adj32_002)

Nicolas Bousserez wrote:

the code crashes in "input_adj_mod.f" when doing the checking of observation settings (see line 2027 in my code). It is apparently due to "IDCSPEC_ADJ", which is not allocated/defined yet at the time the program executes those lines. I'm skipping these checks for now then. Note: The code crashes also when using "TES_O3_OBS" instead of "SCIA_DAL_NO2_OBS", so it is not related to what I implemented apparently.

The following block of code has been moved from input_adj_mod.f to subroutine INIT_CSPEC_ADJ in adj_arrays_mod.f, as it can only be called after the CSPEC variables have been initialized.

#elif defined ( TES_O3_OBS ) 
     ! Since the O3 obs operators will pass adjoints back 
     ! to CSPEC via CSPEC_AFTER_CHEM_ADJ, we need to make sure that 
     ! these species are listed as observed species
     FOUND = .FALSE.
     DO N = 1, NOBS_CSPEC

        IF ( TRIM( NAMEGAS( IDCSPEC_ADJ(N) ) ) == 'O3' ) THEN
           FOUND = .TRUE.

     IF ( .not. FOUND ) THEN

        CALL ERROR_STOP( ' Need to list O3 as observed species',
    &                    ' input_adj_mod ' )

--Daven 16:38, 6 January 2012 (EST)

Retire CSPEC_NO2_ADJ (adj32_003)

Nicolas Bousserez wrote:

In "chemdr_adj.F" I removed this part:
#if defined( SCIA_KNMI_NO2_OBS ) || defined( SCIA_DAL_NO2_OBS )
! Apply forcing from satellite observations
CSPEC_NO2_ADJ(:) = 0d0
because it seems like I'm already updating the adjoint of NO2 in "sciadal_no2_obs_mod.f" through CSPEC_AFTER_CHEM_ADJ. I would like to be sure that it is correct though.

Daven Henze wrote:

Yes, that is correct. I didn't want to take it out incase you had built your obs operator using CSPEC_NO2_ADJ. But since you used the more generic CSPEC_AFTER_CHEM_ADJ, then we can delete this section entirely and use the section higher up in this file ( IF (LCSPEC_OBS) ...) to apply the forcing from your obs operator to CSPEC_ADJ.

Now use CSPEC_AFTER_CHEM_ADJ instead of CSPEC_NO2_ADJ in chemdr_adj.f.

--Daven 18:04, 6 January 2012 (EST)

Fix to NaNs in STT_ADJ everywhere after adjoint of partitioning (adj32_004)

The forward gas-phase chemical solver may sometimes need to repeat the internal integration of the reaction rate ODEs within a single grid cell. When this happens, we now only redo the model calculation during the forward model run, not the adjoint, as this can lead to a corrupt CSPEC_ADJ array globally (not sure how that happens, but it does), eventually triggering a check for NaNs after the partitioning adjoint. See updates to chemistry_mod.f and gckpp_adj_Integrator.f90. Thanks to Richard Xu for helping pinpoint this bug.

--Daven 18:03, 6 January 2012 (EST)

Fix to NaNs in sulfate aerosol adjoint tracers (adj32_005)

Owing to cancelation of errors, it's possible to seemingly randomly end up with NaNs coming out of CHEM_SO2_ADJ . This can happen if L3 ends up equaling 1d0 exactly. We now trap for this and skip in the rate event this happens.

--Daven 18:03, 6 January 2012 (EST)

Remove superfluous debugging which had LAVHRRLAI set to FALSE

Noted by Fabien. Upated in input_mod.f.

--Daven 18:07, 6 January 2012 (EST)

Fix for reseting deposition adjoints (adj32_006)

Fabien Paulot wrote:

It seems it's "a bug" in the code. To keep track of the dry deposition flux of HNO3 I use dryhno3 (in cspec array). The CSPEC_ADJ(DRYHNO3,:) should be reset to zero at every time step but it is not. That's because you are saving in chk_cspec before it is reset to zero. Easy fix, just force to zero in chemdr_adj before the forcing is applied et voila, it works (the adjoint does not diverge any more after one hour). It should not affect anything as long as there is no forcing on dryhno3 (or any dry dep counters for that matter).

We add the following to reset the dry dep adjoints in chemdr_adj.f:

      DO N = 1,NUMDEP
         NK = NTDEP(N)
         IF (NK.NE.0) THEN
            JJ = IRM(NPRODLO+1,NK,NCS)
            IF (JJ.GT.0) THEN
               CSPEC_ADJ(:,JJ) = 0.0D0

--Daven 18:34, 6 January 2012 (EST)

Update comments for STT_ADJ unit conversion

Martin Keller wrote:

I'm currently guessing that CONVERT_UNITS(1,etc.) does something like:
STT:     kg -> v/v

STT_ADJ: J/[v/v] -> J/kg
and similarly CONVERT_UNITS(2,etc.)
STT:     v/v -> kg

STT_ADJ: J/kg -> J/[v/v]
From a dimensional analysis this would seem plausible. It would also explain some of the errors I'm currently fighting with. If this is correct, then a short comment should perhaps be added to the code.

The comments in geos_chem_adj_mod.f have been updated surrounding calls to CONVERT_UNITS.

--Daven 18:43, 6 January 2012 (EST)

Only write upbd checkpoint files if necessary (adj32_007)

Thanks to Farhan Akhtar for pointing this out. We now check to make sure writing these checkpoint files is necessary (i.e., it is an adjoint run, and not just a pseudo obs calculation) in linoz_mod.f:

      ! Now see if checkpointing is necessary (fa, dkh, 01/06/12, adj32_007) 
      IF ( LADJ .and. N_CALC > 0 ) THEN

--Daven 18:54, 6 January 2012 (EST)

Fix for reading TES diagnostic flags from input file (adj32_010)

The problem was that if the CO diagnostics (LDCOSAT) were set to false, then the line counting would be off for any subsequent lines read in input.gcadj. The file input_adj_mod.f has been updated to work properly regardless of the value of LDCOSAT.

--Daven 18:59, 8 January 2012 (EST)

Outstanding issues not yet resolved in v32

This is just a quick list of the features that we're working on putting into v32. As they are integrated and validated, we'll move these up to section 2 and expand the descriptions.

  • CH4 adjoint (Wecht, Harvard)
  • Strat fluxes and chemistry (Lee, Colorado)
  • NH3 updates (Zhu, Colorado)
  • O3 IRK (Henze, Colorado)
  • Population weighting (Vogel, Colorado)
  • NO2/SO2 obs operators (Bousserez, Dalhousie)
  • ANISO (Capps, Georgia Tech)

Things to potentially add to the pipeline for the v32 release:

  • Nested grid full chemistry (Jiang, Toronto)?
  • Reaction rates / deposition / updated chemistry (Paulot, Harvard)?
  • ETOH simulations (Millet, UMN)