GEOS-Chem Adjoint v32

From Geos-chem
Revision as of 05:23, 13 January 2012 by Daven (Talk | contribs) (Updates to adjoint model)

Jump to: navigation, search


Under development. Expected release date 02/01/2012

What's new in this version

Updates to forward model

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 buy). Common blocks CMN_SIZE and CMN_GCTM are reinstated for backwards compatibility.

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

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)

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)

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.

  • Dust adjoint (UNL, SNU)
  • CH4 adjoint (Wecht, Harvard)
  • Inverse Hessian approximation (Tang)
  • 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)?
  • IMPROVE observation operators (UCLA)?
  • ETOH simulations (Millet, UMN)