GEOS-Chem Adjoint v32: Difference between revisions

From Geos-chem
Jump to navigation Jump to search
 
(43 intermediate revisions by 2 users not shown)
Line 1: Line 1:
== Overview ==
== Overview ==
Under development.  Expected release date 02/01/2012
Under development.  Expected release date 02/15/2012


*Previous version: [[GEOS-Chem_Adjoint_v31]]
*Previous version: [[GEOS-Chem_Adjoint_v31]]
Line 102: Line 102:
===== nei05 and input_mod.f =====
===== nei05 and input_mod.f =====
<tt>nei2005_anthro_mod.f</tt> and <tt>input_mod.f</tt> have been updates, see [[GEOS-Chem_Adjoint_v32#Fixes_for_minor_issues_affecting_nested-grid_simulations_.28v9_01_02e.29]]
<tt>nei2005_anthro_mod.f</tt> and <tt>input_mod.f</tt> have been updates, see [[GEOS-Chem_Adjoint_v32#Fixes_for_minor_issues_affecting_nested-grid_simulations_.28v9_01_02e.29]]
==== Update CH4 simulation (adj32_023) ====
[[GEOS-Chem_Adjoint_v32#Include_adjoint_of_CH4_simulation_.28adj32_023.29]]
==== Update stratospheric chemistry and prod/loss rates  (adj32_025) ====
[[GEOS-Chem_Adjoint_v32#Add_strat_production_and_loss_rates_as_active_control_parameters_.28adj32_025.29]]


=== Bug fixes in forward model ===
=== Bug fixes in forward model ===
Line 107: Line 113:
----
----


==== Ionic strength below zero in rpmares (adj32_001) ====
==== Prevent overflow / underflow in rpmares (adj32_001) ====
 
===== Ionic strength below zero =====
On prospero, we often get the message: <tt>Ionic strength below zero...negative concentrations</tt>.
On prospero, we often get the message: <tt>Ionic strength below zero...negative concentrations</tt>.


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
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        


          !IF ( CRUTES(1) .ne. ABS(CRUTES(1) ) ) CRUTES(1)=1d9
        ! Trap for negative CRUTES (dkh, 01/06/12, adj32_001))
          !IF ( CRUTES(2) .ne. ABS(CRUTES(2) ) ) CRUTES(2)=1d9
        IF ( CRUTES(1) < 0d0 ) THEN
          !IF ( CRUTES(3) .ne. ABS(CRUTES(3) ) ) CRUTES(3)=1d9
            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
            CRUTES(1) = MIN( CRUTES(1), CRUTES(2), CRUTES(3) )
            IF ( CRUTES(1) < 0d0 ) THEN
              print*, ' CRUTES = ', CRUTES        , NR
              CALL ERROR_STOP( ' checking for neg val failed',
    &                      ' rpmares_mod.f ' )
            ENDIF
        ENDIF


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


            ! reinstate (dkh, 02/12/12, adj32_001)
            IF (HPLUS <= 0d0) GOTO 100
--[[User:Daven|Daven]] 14:58, 6 January 2012 (EST)
--[[User:Daven|Daven]] 14:58, 6 January 2012 (EST)
===== RPMARES: not safe to divide...exit =====
When running with LADJ = F, this error would arise.  The reason is not entirely clear, but with the exact same inputs, it sometimes goes away by adding a print statement here or there, so is likely owing to overflow / underflow within RPMARES somewhere.  The RPMARES_FORADJ routine is not subject to this instability, so now we move the LADJ flag within RPMARES_FORADJ and always use this routine.
--[[User:Daven|Daven]] 12:32, 12 February 2012 (EST)
===== Compiler level fix =====
'''''[mailto:fabienpaulot@gmail.com Fabien Paulot]''''' wrote
:Hi Daven,
:I am sorry to bug you again about that issue but this time I may have a solution. It seems the problem with RPMARES is related with the compilation option.
:I don't get the error ("not safe to divide") when I use the following option
:FFLAGS = -cpp -w -O2 -auto -noalign -convert big_endian -fp-model precise and remove the -mp option from the F90 compiler options.
:(Note switching to O2 is not the reason why it's working but Bob suggested not to use O3).
We now implement the -fp-model precise flag and do not use the -mp flag.  The latter was originally implemented to allow KPP code to run on 32-bit systems. 
This compiler level fix does address the two issues with rpmares noted above on prospero.  But we will leave those patches in the code for now, as they don't hurt, and we'll see how this plays out on other platforms.
--[[User:Daven|Daven]] 00:34, 1 March 2012 (EST)


==== Updated lightning parameterization and fix for cloud-top-height algorithm (v9-01-01) ====
==== Updated lightning parameterization and fix for cloud-top-height algorithm (v9-01-01) ====
Line 191: Line 236:
==== Fixes for tagged CO EMEP emissions (adj32_019) ====
==== Fixes for tagged CO EMEP emissions (adj32_019) ====
To allow for tagged-CO simulation to use the EMEP emissions, use hardwired tracer numbers in emep_mod.f  This should eventually be migrated into the standard forward model code.  
To allow for tagged-CO simulation to use the EMEP emissions, use hardwired tracer numbers in emep_mod.f  This should eventually be migrated into the standard forward model code.  
--[[User:Daven|Daven]] 15:23, 9 February 2012 (EST)
==== Fix for nested boundary condition (adj32_021) ====
In <tt>tpcore_bc_mod.f</tt>, subroutine DO_WINDOW_TPCORE_BC_05x0666, the OpenMP statements were causing problems, so remove these for now.
--[[User:Daven|Daven]] 15:23, 9 February 2012 (EST)


==== Fix for nested boundary condition (adj32_022) ====
==== Fix for nested boundary condition (adj32_022) ====
In subroutine DO_WINDOW_TPCORE_BC_05x0666, the OpenMP statements were causing problems, so remove these for now.
Updated to<tt>transport_mod.f</tt> and <tt>diag3.f</tt>.  This issue has been corrected in the standard forward model with the following note:
 


      ! Note: the mass flux diagnostic arrays (MASSFLEW, MASSFLNS and MASSFLUP)
      ! are incremented upside-down (level 1 = top of the atmosphere).
      ! The levels order is reversed only when written out in diag3.f
      ! (ccc, 3/8/10),


--[[User:Daven|Daven]] 15:23, 9 February 2012 (EST)


--[[User:Daven|Daven]] 19:51, 16 January 2012 (EST)
--[[User:Daven|Daven]] 19:51, 16 January 2012 (EST)
Line 337: Line 393:
*sciabr_co_obs_mod
*sciabr_co_obs_mod


--[[User:Daven|Daven]] 15:01, 18 January 2012 (EST)
==== Include adjoint of CH4 simulation (adj32_023) ====
Updates from Kevin Wecht.  Modifications in the following files:
* checkpt_mod.f
* adj_arrays_mod.f
* chemistry_adj_mod.f
* cleanup_adj.f
* define_adj.h
* emissions_adj_mod.f
* geos_chem_adj_mod.f
* input_adj_mod.f
* a3_read_mod.f
* dao_mod.f
* geos_chem_mod.f
* global_ch4_mod.f
* global_oh_mod.f
* CMN_DIAG and diag3.f, diag_mod.f, initialize.f, ndxx_setup.f
* julday_mod.f (fwd model bug fix)
New files
* netcdf_util_mod.f
* global_ch4_adj_mod.f
* tes_ch4_mod.f
* scia_ch4_mod.f
* leo_ch4_mod.f
* mem_ch4_mod.f
*geocape_ch4_mod.f
===== Add support for ND19, ND59 and ND60 =====
Modified CMN_DIAG and diag3.f, diag_mod.f, initialize.f, ndxx_setup.f
===== Add support for OH_LOSS =====
Modified diag_oh_mod.f
--[[User:Daven|Daven]] 13:26, 12 February 2012 (EST)
==== Add population weighted cost function option (adj32_024) ====
Update from Steven Vogel adds a new option for cost function evaluation, LPOP_UGM3, and a new file <tt>population_mod.f</tt>. 
--[[User:Daven|Daven]] 01:10, 14 February 2012 (EST)
==== Add strat production and loss rates as active control parameters (adj32_025) ====
Hyung-Min Lee has developed the adjoint with respect to production and loss control parameters, LADJ_STRAT, based on the new stratospheric fluxes from Lee Murray [[Stratospheric_chemistry]]. Also fix debug printing for JLOP.
Updated files include:
*logical_adj_mod.f
*adj_arrays_mod.f
*checkpoint_mod.f
*checkpt_mod.f
*chemdr_adj.f
*geos_chem_adj_mod.f
*input_adj_mod.f
*inverse_driver.f
*inverse_mod.f
*linoz_adj_mod.f (completely replaced)
*gamap_mod.f
*partition_mod.f
*tacer_mod.f
*strat_chem_adj_mod.f (added)
fwd:
*logical_mod.f
modified:
*chemdr.f
*geos_chem_mod.f (remove UPBDFLX_NOY and DO_UPBDFLX)
*input_mod.f
*linoz_mod.f (replace w/v9)
*strat_chem_mod.f (add from v9)
*Makefile.ifort.* (because new strat chem uses netcdf)
--[[User:Daven|Daven]] 15:22, 14 February 2012 (EST)
==== Calculate array sizes dynamically in CINSPECT (adj32_026) ====
The arrays VAR_FD and RCONST_FD have been moved to adj_arrays_mod.f and are now dynamically allocated.  The following error trap (which would frequently cause problems for new users) is no longer relevant:


--[[User:Daven|Daven]] 15:01, 18 January 2012 (EST)
      IF ( NCHEM > 1000 ) THEN
        CALL ERROR_STOP( 'Need to boost NCHEM max', 'chemistry_mod.f')
      ENDIF
 
--[[User:Daven|Daven]] 16:22, 23 February 2012 (EST)


=== Bug fixes in adjoint model ===
=== Bug fixes in adjoint model ===
Line 483: Line 624:
F_OF_PBL and ED_CO were declared thread PRIVATE in several parallel do loops.  
F_OF_PBL and ED_CO were declared thread PRIVATE in several parallel do loops.  


--[[User:Daven|Daven]] 18:59, 8 January 2012 (EST)
==== Bug fix for adjoints in nested cushions (adj32_020) ====
In <tt>transport_mod.f</tt>, now call DO_WINDOW_TPCORE_BC_ADJ to set the adjoint values to zero within the cushion.


--[[User:Daven|Daven]] 18:59, 8 January 2012 (EST)
--[[User:Daven|Daven]] 15:26, 9 February 2012 (EST)


== Outstanding issues not yet resolved in v32 ==
== 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.
New features to add to the pipeline for the v33 release:
*CH4 adjoint (Wecht, Harvard)
*Reaction rates / deposition / updated chemistry (Paulot, Harvard)
*Strat fluxes and chemistry (Lee, Colorado)
*ETOH simulations (Millet, UMN)
*NH3 updates (Zhu, Colorado)
*NH3 updates (Zhu, Colorado)
*O3 IRK (Henze, Colorado)
*Population weighting (Vogel, Colorado)
*NO2/SO2 obs operators (Bousserez, Dalhousie)
*NO2/SO2 obs operators (Bousserez, Dalhousie)
*ANISO (Capps, Georgia Tech)
*ANISO (Capps, Georgia Tech)
*Update from strat_chem_201106 to strat_chem_201109 in strat_chem_mod.f ?
Some bug fixes that need implementing
* Fix bug in Makefile.ifort.hdf, .lidort and .netcdf: need to list define.h in set of include files for strat_chem_mod.f
* Trap for xFD not in trop in geos_chem_mod.f, line 1622:  add "IF JLOP(IFD,JFD,LFD) > 0" to trap for LFD not in tropospher
* Fix gckpp_adj_Integrator.f90
* Increase magnitude of filter in partition_adj.f
* Fix NH3 seasonality -- check with Steve
* Fix line 4737 of lidort_mod.f (Jamin Koo, 04/20/2012)
* Remove M0 and M0_ADJ from adj_arrays_mod.f (not used) hml


Things to potentially add to the pipeline for the v32 release:
*Reaction rates / deposition / updated chemistry (Paulot, Harvard)?
*ETOH simulations (Millet, UMN)


And some bugs of note with no fix yet
And some bugs of note with no fix yet
*Can't specify IFD, JFD using latitudes with nested domain
*Can't specify IFD, JFD using latitudes with nested domain
*Can't use NEI emissions with tagged-CO

Latest revision as of 23:07, 3 July 2013

Overview

Under development. Expected release date 02/15/2012

What's new in this version

Updates to forward model


Update tagged CO simulation (adj32_017)

Several updates from Zhe Jiang have been implemented for the tagged-CO implementation. See also changes to the adjoint model, GEOS-Chem_Adjoint_v32#Update_tagged_CO_simulation_.28adj32_017.29

Update to tagged_co_mod.f

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.

Also, F_OF_PBL and ED_CO were declared thread PRIVATE in several parallel do loops.


Updates to time_mod.f for tagged-CO

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 )
      ENDIF
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)

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.

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

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

emissions_mod.f

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

file_mod.f

file_mod.f was updated to include IU_BC_05x06.

CMN_SIZE

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

commsoil.h

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

define.h

define.h was updated to include the NESTED_SD option.

tpcore_bc_mod.f

tpcore_bc_mod.f has been updated to include

  • OPEN_BC_FILE_05x0666
  • SAVE_GLOBAL_TPCORE_BC_05x0666
  • DO_WINDOW_TPCORE_BC_05x0666
  • CLEAN_WINDOW_TPCORE_BC_05x0666
  • READ_WINDOW_TPCORE_BC_05x0666
  • ITS_TIME_FOR_BC_05x0666INIT_TPCORE_BC_05x0666

and the BC_05x06 array.

regrid_1x1_mod.f

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

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

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

transport_mod.f

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

co2_mod.f

co2_mod.f has been updated to support nested domains.


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


nei05 and input_mod.f

nei2005_anthro_mod.f and input_mod.f have been updates, see GEOS-Chem_Adjoint_v32#Fixes_for_minor_issues_affecting_nested-grid_simulations_.28v9_01_02e.29

Update CH4 simulation (adj32_023)

GEOS-Chem_Adjoint_v32#Include_adjoint_of_CH4_simulation_.28adj32_023.29

Update stratospheric chemistry and prod/loss rates (adj32_025)

GEOS-Chem_Adjoint_v32#Add_strat_production_and_loss_rates_as_active_control_parameters_.28adj32_025.29

Bug fixes in forward model


Prevent overflow / underflow in rpmares (adj32_001)

Ionic strength below zero

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

        !  Trap for negative CRUTES (dkh, 01/06/12, adj32_001)) 
        IF ( CRUTES(1) < 0d0 ) THEN
           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
           CRUTES(1) = MIN( CRUTES(1), CRUTES(2), CRUTES(3) )
           IF ( CRUTES(1) < 0d0 ) THEN
              print*, ' CRUTES = ', CRUTES         , NR
              CALL ERROR_STOP( ' checking for neg val failed',
    &                       ' rpmares_mod.f ' )
           ENDIF
        ENDIF

In addition, reinstate the following

            ! reinstate (dkh, 02/12/12, adj32_001) 
            IF (HPLUS <= 0d0) GOTO 100

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

RPMARES: not safe to divide...exit

When running with LADJ = F, this error would arise. The reason is not entirely clear, but with the exact same inputs, it sometimes goes away by adding a print statement here or there, so is likely owing to overflow / underflow within RPMARES somewhere. The RPMARES_FORADJ routine is not subject to this instability, so now we move the LADJ flag within RPMARES_FORADJ and always use this routine.

--Daven 12:32, 12 February 2012 (EST)


Compiler level fix

Fabien Paulot wrote

Hi Daven,
I am sorry to bug you again about that issue but this time I may have a solution. It seems the problem with RPMARES is related with the compilation option.
I don't get the error ("not safe to divide") when I use the following option
FFLAGS = -cpp -w -O2 -auto -noalign -convert big_endian -fp-model precise and remove the -mp option from the F90 compiler options.
(Note switching to O2 is not the reason why it's working but Bob suggested not to use O3).


We now implement the -fp-model precise flag and do not use the -mp flag. The latter was originally implemented to allow KPP code to run on 32-bit systems.

This compiler level fix does address the two issues with rpmares noted above on prospero. But we will leave those patches in the code for now, as they don't hurt, and we'll see how this plays out on other platforms.

--Daven 00:34, 1 March 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: http://people.seas.harvard.edu/~ltmurray/LNOx.v9-01-01.Release.Notes.pd

--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)

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,  
      SPECIES_ID_SAVE = (/ IDTNOX,  IDTCO,   IDTSO2,  IDTSO4, IDTNH3,
     $                     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
 
  ENDIF

--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

Fixes for tagged CO EMEP emissions (adj32_019)

To allow for tagged-CO simulation to use the EMEP emissions, use hardwired tracer numbers in emep_mod.f This should eventually be migrated into the standard forward model code.

--Daven 15:23, 9 February 2012 (EST)

Fix for nested boundary condition (adj32_021)

In tpcore_bc_mod.f, subroutine DO_WINDOW_TPCORE_BC_05x0666, the OpenMP statements were causing problems, so remove these for now.

--Daven 15:23, 9 February 2012 (EST)

Fix for nested boundary condition (adj32_022)

Updated totransport_mod.f and diag3.f. This issue has been corrected in the standard forward model with the following note:

     ! Note: the mass flux diagnostic arrays (MASSFLEW, MASSFLNS and MASSFLUP)
     ! are incremented upside-down (level 1 = top of the atmosphere).
     ! The levels order is reversed only when written out in diag3.f
     ! (ccc, 3/8/10),

--Daven 15:23, 9 February 2012 (EST)

--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


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

Do not call GET_BOUNDING_BOX in input_adj_mod.f

When calculating IFD and JFD from a user specified latitude and longitude, the call to GET_BOUNDING_BOX causes a crash, so we need to modify this part of the code to work on the nested domain. For now, just re-order the code so that this is only called if USEINDEX is false, and run the nested code with USEINDEX true.

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 15:01, 18 January 2012 (EST)

Updates for tagged CO simulation (adj32_017)

Several updates from Zhe Jiang have been implemented for the tagged-CO implementation. See also changes to the forward model, GEOS-Chem_Adjoint_v32#Updates_for_tagged_CO_simulation__.28adj32_017.29

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.

Cleanup and merge CALC_APRIORI

Updates 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

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 15:01, 18 January 2012 (EST)

Include adjoint of CH4 simulation (adj32_023)

Updates from Kevin Wecht. Modifications in the following files:

  • checkpt_mod.f
  • adj_arrays_mod.f
  • chemistry_adj_mod.f
  • cleanup_adj.f
  • define_adj.h
  • emissions_adj_mod.f
  • geos_chem_adj_mod.f
  • input_adj_mod.f
  • a3_read_mod.f
  • dao_mod.f
  • geos_chem_mod.f
  • global_ch4_mod.f
  • global_oh_mod.f
  • CMN_DIAG and diag3.f, diag_mod.f, initialize.f, ndxx_setup.f
  • julday_mod.f (fwd model bug fix)

New files

  • netcdf_util_mod.f
  • global_ch4_adj_mod.f
  • tes_ch4_mod.f
  • scia_ch4_mod.f
  • leo_ch4_mod.f
  • mem_ch4_mod.f
  • geocape_ch4_mod.f
Add support for ND19, ND59 and ND60

Modified CMN_DIAG and diag3.f, diag_mod.f, initialize.f, ndxx_setup.f

Add support for OH_LOSS

Modified diag_oh_mod.f


--Daven 13:26, 12 February 2012 (EST)

Add population weighted cost function option (adj32_024)

Update from Steven Vogel adds a new option for cost function evaluation, LPOP_UGM3, and a new file population_mod.f.

--Daven 01:10, 14 February 2012 (EST)

Add strat production and loss rates as active control parameters (adj32_025)

Hyung-Min Lee has developed the adjoint with respect to production and loss control parameters, LADJ_STRAT, based on the new stratospheric fluxes from Lee Murray Stratospheric_chemistry. Also fix debug printing for JLOP.

Updated files include:

  • logical_adj_mod.f
  • adj_arrays_mod.f
  • checkpoint_mod.f
  • checkpt_mod.f
  • chemdr_adj.f
  • geos_chem_adj_mod.f
  • input_adj_mod.f
  • inverse_driver.f
  • inverse_mod.f
  • linoz_adj_mod.f (completely replaced)
  • gamap_mod.f
  • partition_mod.f
  • tacer_mod.f
  • strat_chem_adj_mod.f (added)


fwd:

  • logical_mod.f

modified:

  • chemdr.f
  • geos_chem_mod.f (remove UPBDFLX_NOY and DO_UPBDFLX)
  • input_mod.f
  • linoz_mod.f (replace w/v9)
  • strat_chem_mod.f (add from v9)
  • Makefile.ifort.* (because new strat chem uses netcdf)

--Daven 15:22, 14 February 2012 (EST)

Calculate array sizes dynamically in CINSPECT (adj32_026)

The arrays VAR_FD and RCONST_FD have been moved to adj_arrays_mod.f and are now dynamically allocated. The following error trap (which would frequently cause problems for new users) is no longer relevant:

     IF ( NCHEM > 1000 ) THEN 
        CALL ERROR_STOP( 'Need to boost NCHEM max', 'chemistry_mod.f')
     ENDIF 

--Daven 16:22, 23 February 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.
        ENDIF

     ENDDO
     IF ( .not. FOUND ) THEN

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

--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_ADJ(:,IDNO2) = CSPEC_ADJ(:,IDNO2) + CSPEC_NO2_ADJ(:)
CSPEC_NO2_ADJ(:) = 0d0
#endif
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:

!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( N, NK, JJ )
      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
            ENDIF
         ENDIF
      ENDDO
!$OMP END PARALLEL DO

--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
         CALL MAKE_UPBDFLX_CHKFILE( NYMD, NHMS, TAU )
      ENDIF 

--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.


Bug fix for unit conversion in aircraft and chemical source terms for CO2 adjoint (adj32_018)

The units of the emissions scaling factor adjoints for the aircraft and chemical source terms were incorrect. These sources are in units per volume, not area, so their conversion needed to be

             E_CO2_ADJ      = E_CO2_ADJ * BOXVL(I,J,L)
     &                      * DTSRCE    / XNUMOL_CO2

instead of

           !!E_CO2_ADJ      = E_CO2_ADJ * A_CM2 * DTSRCE / XNUMOL_CO2

Also, the ordering of the application of the CO2 emissions scaling factors was made consistent across all sectors such that the ND04 diagnostic is correct (previously the scaling had been applied after the diagnostic for ocean exchange, balanced biosphere, and CO2 surface correction for CO oxidation).

--Daven 10:32, 9 February 2012 (EST)

Bug fix for parallel loops in tagged_co_mod.f (adj32_017)

F_OF_PBL and ED_CO were declared thread PRIVATE in several parallel do loops.

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

Bug fix for adjoints in nested cushions (adj32_020)

In transport_mod.f, now call DO_WINDOW_TPCORE_BC_ADJ to set the adjoint values to zero within the cushion.

--Daven 15:26, 9 February 2012 (EST)

Outstanding issues not yet resolved in v32

New features to add to the pipeline for the v33 release:

  • Reaction rates / deposition / updated chemistry (Paulot, Harvard)
  • ETOH simulations (Millet, UMN)
  • NH3 updates (Zhu, Colorado)
  • NO2/SO2 obs operators (Bousserez, Dalhousie)
  • ANISO (Capps, Georgia Tech)
  • Update from strat_chem_201106 to strat_chem_201109 in strat_chem_mod.f ?

Some bug fixes that need implementing

  • Fix bug in Makefile.ifort.hdf, .lidort and .netcdf: need to list define.h in set of include files for strat_chem_mod.f
  • Trap for xFD not in trop in geos_chem_mod.f, line 1622: add "IF JLOP(IFD,JFD,LFD) > 0" to trap for LFD not in tropospher
  • Fix gckpp_adj_Integrator.f90
  • Increase magnitude of filter in partition_adj.f
  • Fix NH3 seasonality -- check with Steve
  • Fix line 4737 of lidort_mod.f (Jamin Koo, 04/20/2012)
  • Remove M0 and M0_ADJ from adj_arrays_mod.f (not used) hml


And some bugs of note with no fix yet

  • Can't specify IFD, JFD using latitudes with nested domain
  • Can't use NEI emissions with tagged-CO