Stratospheric chemistry

On this page we list information about the stratospheric chemistry mechanism in GEOS-Chem.

For the influence of stratospheric composition on tropospheric photolysis rates, see our Photolysis mechanism wiki page.

Overview

Though GEOS-Chem is predominantly used for simulating tropospheric chemistry, it is necessary to deal with species that may have net stratospheric sources and transport to the troposphere (e.g., NOx), or net stratospheric sinks of gases transported from the troposphere (e.g., CO). An updated linearized stratospheric chemistry scheme has been prepared for implementation into GEOS-Chem v9-01-03 to complement the linearized ozone chemistry afforded by Linoz.

Authors and collaborators:

• Lee Murray (formerly Harvard, now University of Rochester)

New implementation of linearized stratospheric chemistry

An updated linearized stratospheric chemistry scheme has been implemented into GEOS-Chem v9-01-03o. It currently works with GCAP / GEOS-4 / GEOS-5 / MERRA/ GEOS-FP simulations at all resolutions, for full chemistry, SOA, dicarbonyl, isoprene, H2/HD, tagged Ox, and tagged CO.

One may turn on the linearized stratospheric chemistry scheme in the chemistry menu of the input.geos file

        %%% CHEMISTRY MENU %%%  : 
        Turn on Chemistry?      : T
        Use linear. strat. chem?: T
         => Use Linoz for O3?   : T
        Chemistry Timestep [min]: 60
        Read and save CSPEC_FULL: T
        USE solver coded by KPP : F

At the beginning of a new month, the model reads in the archived 3D monthly mean production rate, P (v/v/s), and loss frequency, k (s-1), for each species. The new concentration at the end of each chemistry time step is then determined by solving

        ! Simple analytic solution to dM/dt = P - kM over [0,t]
        if ( k .gt. 0d0 ) then
           STT(I,J,L,N) = M0 * exp(-k*t) + (P/k)*(1d0-exp(-k*t))
        else
           STT(I,J,L,N) = M0 + P*t
        endif

This is applied for each pertinent species to every grid box for which the tropospheric chemistry solver is not applied in a given chemistry time step.

Also, monthly mean stratospheric OH fields are read in and used to generated decay for certain species (currently, CH3Br, CHBr3, CH2Br2).

The bromine species are hard-wired to be ignored for now. Br_y in the stratosphere is prescribed using monthly mean nighttime and daytime concentrations from the GEOS CCM as in its original (v9-01-03m) implementation. For more information, see the bromine chemistry mechanism wiki page.

Note: Linoz is still the recommended option for stratospheric ozone, and ozone prod/loss rates were only implemented here for completeness. If one turns off Linoz, then Synoz is applied in the code.

--Bob Y. 12:28, 25 April 2014 (EDT)

Prod/loss rates from UCX

This update was included in v11-02e (approved 24 Mar 2018).

The monthly-mean production, rates, loss frequencies and concentrations for 2013 have been archived from 1-year benchmark GEOS-Chem v11-02d-Run2 using the UCX chemistry mechanism. These values will be used for all the linearized stratospheric chemistry of all species in GEOS-Chem v11-02e and later.

--Melissa Sulprizio (talk) 17:08, 20 November 2017 (UTC)

Prod/loss rates from GMI

Preparation of monthly mean 3D prod/loss rates

The production rates (v/v/s), loss frequencies (s^-1) and concentrations (v/v) have been archived from the GMI "Combo" model simulations, which one can think of as having the GEOS-Chem troposphere + the GMI full stratospheric chemistry model. For GEOS4/GEOS-Chem simulations, long-term monthly means were calculated from the "Aura4" GMI monthly concentrations, which used GEOS-4 assimilated met fields for 2004-2006. For GEOS5/ or MERRA/GEOS-Chem simulations, long-term monthly means were determined from the "MERRA" GMI rerun (containing the Synoz bug fix in the GMI model) using MERRA met fields for 2004-2010.

The GMI model outputs as diagnostic the monthly mean rate of each kinetic or photolysis reaction. That output was aggregated to get the total production rate (v/v/s). Photolysis frequencies are also available, but kinetic loss frequencies had to be estimated from the burden/loss rate, to determine the mean loss frequencies (s^-1).

--Bob Y. 12:28, 25 April 2014 (EDT)

Species with linearized stratospheric chemistry

Numbers indicate number of GMI Combo model reactions aggregated, see [reaction details document] for more information

Species Name	Photolysis Prod	Kinetic Prod	Photolysis Loss	Kinetic Loss
NOx	14	39		33
Ox	9	8		37
PAN		1	1	1
CO	16	20
ALK4				2
ISOP				3
HNO3		38	1	2
H2O2		4	1	2
ACET	1	7	1	1
MEK	1	18	1	2
ALD2	4	26	2	2
RCHO	9	33	1	2
MVK	1	7	3	3
MACR	1	7	2	4
PMN		1		3
PPN		1		1
R4N2		4	1	1
PRPE	1	1		3
C3H8				2
CH2O	10	70	2	5
C2H6				2
N2O5		1	2	4
HNO4		1	2	2
MP		1	1	2
GLYX		5	3	2
MGLY	4	23	2	2
GLYC	2	10	1	1
HAC	2	12	1	1
GPAN		1		1
ACTA		20		1
RIP		2	1	1
MAP		1	1	1
CH4	1			5

Additionally available species

The following species from the GMI model have also had their monthly mean production rates (v/v/s), loss frequencies (s-1), and concentrations (v/v) archived for GEOS-Chem. The linearized chemistry will automatically have prod/loss applied if a tracer sharing one of the following names is introduced into the model's STT array of transported tracers and DO_STRAT_CHEM is called.

H, H2, HCOOH, H2O, HO2, MO2, MOH, N, N2O, O, O1D, OH, Br, BrCl, BrO, BrONO2, HBr, HOBr, Cl, Cl2, ClO, Cl2O2, ClONO2, HCl, HOCl, OClO, CH3Br, CH3Cl, CH3CCl3, CCl4, CFCl3, CF2Cl2, CFC113, CFC114, CFC115, HCFC22, HCFC141b, HCFC142b, CF2Br2, CF2ClBr, CF3Br, H24O2, A3O2, ATO2, B3O2, EOH, ETO2, ETP, GCO3, GP, IALD, IAO2, IAP, INO2, INPN, ISN1, ISNP, KO2, MAN2, MAO3, MAOP, MCO3, MRO2, MRP, MVN2, PO2, PP, PRN1, PRPN, R4N1, R4O2, R4P, RA3P, RB3P, RCO3, RCOOH, RIO1, RIO2, ROH, RP, VRO2, VRP.

The bromine species are hard-wired to be ignored for now. Br_y in the stratosphere is prescribed using monthly mean nighttime and daytime concentrations from the GEOS CCM, located in the bromine_201205/CCM_stratosphere_Bry data directory.

Also note, the input file may be accessed from anywhere to obtain monthly mean concentrations (e.g., 4D stratospheric OH) for use within the model.

Code structure

All stratospheric-relevant chemistry routines (including Linoz) have been consolidated into a single driving routine DO_STRAT_CHEM in GeosCore/strat_chem_mod.F90, called from DO_CHEMISTRY in GeosCore/chemistry_mod.F for full chemistry, tagged Ox and H2/HD simulations. The tagged CO simulation has been updated in GeosCore/tagged_co_mod.F.

• GeosCore/strat_chem_mod.F: Source code file with stratospheric chemistry. The Linoz routines are called from here.
• GEOS_NATIVE/strat_chem_201206/gmi.clim.*.nc: Species-specific vertical and horizontal resolution dependent files containing mean production rates and loss frequencies from the GMI combo model in the primary data directory

Note, the input file is in the netCDF file format, and therefore also requires GEOS-Chem be compiled with netCDF.

Data I/O via HEMCO in GEOS-Chem v10-01 and newer versions

In GEOS-Chem v10-01 and newer versions, the stratospheric chemistry module (GeosCore/strat_chem_mod.F90) relies upon the HEMCO emissions component to read the netCDF files containing the stratospheric production and loss data archived from the GMI model. Routines GET_RATES and GET_RATES_INTERP have therefore been retired.

We have created new GMI data files (in COARDS-compliant netCDF format) for use with HEMCO. These new data files are contained in the HEMCO data directory tree. For detailed instructions on how to download these data files to your disk server, please see our Downloading the HEMCO data directories wiki post.

--Bob Y. 13:26, 3 March 2015 (EST)

Treatment of the mesosphere in UCX simulations

If using UCX-based simulations, stratospheric chemistry is handled by the chemical mechanism. The simple linearized chemistry in strat_chem_mod.F90 is therefore applied to the mesosphere instead.

Sebastian Eastham wrote:

In the course of the evaluation of the UCX by the GMI group, they actually ended up extending the UCX up into the mesosphere. This removed the need for any mesospheric paramaterizations, and did not involve adding any new chemistry or code. They found that this gave the same results as the GMI model - no big surprise since the GMI and UCX models are very similar, and the GMI model is already run up to the model top. As a result we may want to consider (in v11-03) adding an "extended UCX" option which runs UCX directly up to the model top which will resolve this problem when generating P/L data. It will also allow interested users to avoid reading in all the P/L data (a not-insignificant startup cost) or using mesospheric paramaterizations.

--Melissa Sulprizio (talk) 19:07, 16 March 2018 (UTC)

Original implementation

The original stratospheric chemistry scheme through GEOS-Chem v9-01-02 (v9-01-03n) consisted of:

• A fixed stratosphere-to-tropophere flux of ozone from Synoz (an alternative and recommended linearized ozone chemistry Linoz option introduced in v8-02-04)

And performed a simple linearized chemistry scheme from archived monthly mean production rates, J-values, and OH fields from the 2D stratospheric chemistry model of Hans Schneider and Dylan Jones using GEOS-1 and/or GEOS-STRAT met fields (36 vertical levels) that have since been regridded to subsequent GEOS grid vertical resolutions.

• Fixed monthly zonal mean production rates of NOy in pnoy_200106
• Fixed monthly zonal mean production rates and loss frequencies of CO in pco_lco_200203
• Loss of selected species by reaction with archived monthly zonal mean stratospheric OH fields in stratOH_200203
• Loss of selected species by photolysis using archived monthly zonal mean J-values archived in stratjv_200203

These were variously performed by the subroutines

• upbdflx_mod.f
  • DO_UPBDFLX - Called from main.f before transport, calls Synoz/Linoz ozone.
  • UPBDFLX_NOY - Called first from DO_UPBDFLX for NOy production, and then again after transport in main.f to repartition NOy into NOx and HNO3.
• schem.f - Called just after SMVGEAR/KPP solves tropospheric chemistry in chemdr.f. Performs
  • Photolysis loss
  • CO prod/loss
  • Reaction with OH

Original Species

Species	Photolysis Loss	Loss by OH	Special Treatment
NOx			pnoy_200106 (P&L)
Ox			Synoz (P) or Linoz (P&L)
PAN	X
CO			pco_lco_200203 (P&L)
ALK4		X
ISOP		X
HNO3			pnoy_200106 (P&L)
H2O2	X	X
ACET	X	X
MEK		X
ALD2	X	X
RCHO	X	X
MVK	X	X
MACR	X	X
PMN		X
R4N2	X	X
PRPE		X
C3H8		X
CH2O	X	X
C2H6		X
N2O5	X
HNO4	X	X
MP	X	X

--Bob Y. 12:29, 25 April 2014 (EDT)

STE fluxes

Disable STE calculation printed to log file

This update was included in GEOS-Chem 12.0.0.

Melissa Sulprizio wrote:

There are three methods for computing O3 STE that we’ve used recently in benchmarking:

1. Use O3 STE value found at the end of the GEOS-Chem log file. This value is computed in routine CALC_STE as STE = (P-L) – dStrat/dt where P-L is computed using the species concentration before and after LINOZ.
2. Compute STE as the residual term using GEOS-Chem output: STE = (Initial Mass – Final Mass) – Net POx-LOx – Drydep – Wetdep
3. Compute STE as the total vertical flux of O3 at 100 hPa using GEOS-Chem output.

We compared these three methods back in October when investigating O3 STE in v11-02c for inclusion in our benchmarks (see emails copied below). I also compared the three methods of computing O3 STE with 1-month tropchem simulations at 4x5 and 2x2.5:

                            4x5         2x2.5
     Log file               543 Tg a-1  360 Tg a-1
     Residual               626 Tg a-1  460 Tg a-1
     100 hPa vertical flux  636 Tg a-1  534 Tg a-1

The differences between 4x5 and 2x2.5 are much larger when comparing O3 STE obtained from the log file. I’m wondering if we should remove that output from the log file and use one of the other two methods for reporting O3 STE. I don’t see anything obviously wrong in the calculation of the values printed to the log files, but something strange is clearly going on. It doesn’t seem to be a regridding issue because as far as I can tell.

Daniel Jacob wrote:

I have never liked this "log file" approach to calculating the STE. The reason is that it’s a small difference between huge budget terms in the stratosphere. I think we should get rid of it. The best approach to estimate the STE is with the 100 hPa flux, and there the difference between 4x5 and 2x2.5 makes sense. The residual approach is not great because again it’s a small residual of large numbers but that’s not nearly as bad as for the stratosphere, it has the advantage of closing the tropospheric Ox budget, and the 100 hPa flux is not a perfect measure of STE anyway. So we can keep it.

--Melissa Sulprizio (talk) 11:58, 10 August 2018 (UTC)

Fix for STE calculation

This update was included in v11-02c and approved on 21 Sep 2017.

Barron Henderson wrote:

I think we should make the stratosphere-troposphere-exchange a regular part of our benchmark evaluation. Mat Evans pointed me to the logs and what I found surprised me. From version 9 to present we changed the stratosphere from a net source of ozone (540Tg) to a net sink (-32Tg).

Right now, I suspect artifacts in the calculation STE methodology related to the UCX update (some change could be science). I went through individual versions and found that the reported STE changed from 542 to -213 TgO3 in v10-01c with and without UCX chemistry (see table below).

This raises two questions. First, can someone review the STE calculation methodology (see strat_chem_mod.F90*) and make sure it using the right quantities in the mass balance given the update to UCX. Second, can we get STE added as a quantity for review in the benchmark?

    STE                   TgO3/year
    ----------------------------------------
    v9-02l                 536
    v9-02                  541
    v10-01c-Run0           542
    v10-01c-Run1           -213
    v10-01-public-release  -32.7
    v11-01-public-release  -33996**
    v11-02a                -30
    ----------------------------------------

*strat_chem_mod.F90 - I wonder if the prod/loss terms are not being correctly paired with the stratospheric mass change?

**Obviously v11 had issues with the P/L calculations, so we can ignore that.

Seb Eastham wrote:

Your suspicion is correct. It looks like STE is calculated between times t and t+dt as

     STE = (TEND - DSTRAT)/dt

where TEND is the P(O3) - L(O3) in the "stratosphere" and DSTRAT is the total change in ozone mass above the tropopause between t and t+dt.

However, TEND is calculated as the total increase in ozone reported by the linearized chemistry routines which are only used outside the chemistry grid. Since the UCX was implemented, the chemistry grid is extended to the stratopause, so TEND is just the tiny amount of ozone produced in the mesosphere, while DSTRAT is still calculated as the change in total ozone mass above the tropopause - hence the negative "STE".

I don't think it would take long to modify strat_chem_mod. It would just be a question of storing the chemical tendency of ozone in the stratosphere by accumulating [stratospheric ozone mass after chemistry] - [stratospheric ozone mass before chemistry] at each timestep.

--Melissa Sulprizio (talk) 21:53, 12 May 2017 (UTC)

References

1. Murray, L. T., J. A. Logan, and D. J. Jacob, Interannual variability in tropical tropospheric ozone and OH: the role of lightning, J. Geophys. Res., 118, 1-13, doi:10.1002/jgrd.50857, 2013.
2. Murray, L. T., D. J. Jacob, J. A. Logan, R. C. Hudman, and W. J. Koshak, Optimized regional and interannual variability of lightning in a global chemical transport model constrained by LIS/OTD satellite data, J. Geophys. Res., 117, D20307, doi:10.1029/2012JD017934, 2012.

--Bob Y. 12:26, 25 April 2014 (EDT)

Previous issues that are now resolved

The following issues in the stratospheric chemistry module have now been resolved.

Bug fix: make sure stratospheric BrY concentrations are read properly each month

This update was included in GEOS-Chem 12.5.0, which was released on 09 Sep 2019.

Lyatt Jaeglé and Jiayue Huang wrote:

In running two tropchem multi-year simulations with different start dates we noticed very different stratospheric Bry concentrations, especially in polar regions. We think that this is related to a bug in the HEMCO_Config.rc file which specifies the monthly stratospheric Bry files to be read. Currently it is set to:

     #------------------------------------------------------------------------------
     # --- Stratospheric Bry data from the CCM model  ---
     #------------------------------------------------------------------------------
     * GEOSCCM_Br_DAY        ... BR     2007/$MM/1/0 C xyz pptv * - 60 1
     * GEOSCCM_Br2_DAY       ... BRCL   2007/$MM/1/0 C xyz pptv * - 60 1
     ... etc ...

which means that the files are only read once at the beginning of the run and not updated each month. Changing the time specification from:

     2007/$MM/1/0 

to

     2007/1-12/1/0

for each Bry species leads to updating the files every month and fixes the problem. We suggest that this be changed in the default HEMCO_Config.rc files.

--Bob Yantosca (talk) 14:32, 10 June 2019 (UTC)

Do not allocate memory to array MINIT in strat_chem_mod.F90

This update (Git ID: 729fd2a7) was included in GEOS-Chem 12.3.1, which was released on 08 Apr 2019.

The MINIT array (in module GeosCore/strat_chem_mod.F90) has memory allocated to it, but is no longer actively used. MINIT was used to hold initial mass values for the online diagnostic computed by routine CALC_STE. But it was recently found that CALC_STE was giving erroneous values due to numerical instability. In a prior version, we had commented out the CALC_STE routine, but had forgotten to also comment out the code that allocated the MINIT array.

Therefore, the allocation of MINIT is taking up extra memory, which can potentially impede simulations. It was found in at least one case that nested-grid simulations on the AWS cloud computing platform ran out of memory because of the "useless" memory allocated to MINIT.

For now, we have commented out all instances of MINIT in GeosCore/strat_chem_mod.F90, as well as the subroutine SET_MINIT, where MINIT was initialized. At some further point in time we will delete these commented lines.

ALSO NOTE: For now, we have left the module (logical) variable MINIT_IS_SET because this is set from the GCHP routine gigc_chunk_mod.F90, and also from the interface to GEOS-5. For now this will be a "useless" variable (setting it will have no effect). But we have left this as-is for the time being for fear of breaking the GEOS-Chem to NASA-GEOS interface. We will remove it at a later time.

--Bob Yantosca (talk) 19:46, 8 April 2019 (UTC)

Bug in the GMI stratospheric production rates for HCOOH

This update was validated with the 1-month benchmark simulation v11-01j and approved on 03 Dec 2016

Eloise Marais reported very high (> 1ppm) HCOOH in the mid to upper troposphere in GEOS-Chem v10-01.

Eloise Marais wrote:

I ran GEOS-Chem with the GMI HCOOH prod/loss rates and had a large enhancement in HCOOH in February in the midlatitiudes. It's an isolated band that will presumably shift with time.

I eventually found that it's because GMI stratosphere prod rates far exceed the loss rates. This likely isn't an issue in earlier versions of the model that didn't have GMI prod/loss rates for HCOOH. According to my log file from a simulation using the GMI prod/loss rates the stratosphere is a source of 6000 Tg HCOOH per year.

The standard GC simulation (v10) still has HCOOH set to inactive, so this bug should only be apparent in specialized simulations and when FlexChem is implemented.

Xin Chen in Dylan Millet's group has an updated file of HCOOH GMI prod/loss rates that should be used instead of the current prod/loss rates. We have saved this file into the HEMCO data directory path: HEMCO/GMI/v2015-02/gmi.clim.HCOOH.geos5.2x25.20160912.nc.

• NOTE: We have identified an error in the time stamps in the updated file gmi.clim.HCOOH.geos5.2x25.20160912.nc. Please see this addendum for information about how we fixed this issue.

--Melissa Sulprizio (talk) 21:15, 12 September 2016 (UTC)

Fixed incorrect timestamps in the newest HCOOH prod/loss file

This fix was included in GEOS-Chem v11-01 public release

Seb Eastham wrote:

It looks like the values for the time variable in the GMI prod/loss file for HCOOH are not set correctly. Specifically the file gmi.clim.HCOOH.geos5.2x25.nc (as well as in gmi.clim.HCOOH.geos5.2x25.20160912.nc) has the following entries for its time vector:

     time = 0,744,0,0,0,0,0,0,0,0,0,0

It should instead be:

     time = 0,744,1416,2160,2880,3624,4344,5088,5832,6552,7296,8016

This doesn’t cause a problem for HEMCO, but GCHP expects the time entries to be accurate and chokes on this file (by design).

[Also ...] it looks like there is now another problem with the GMI HCOOH loss/prod frequencies, specifically with the dimension information. I noticed that there were significantly more dimension variables than had previously been the case, and this was causing dimension identification issues [in GCHP]....I also generated a simpler file, using the ISOP prod/loss frequency file as a template (as all the other GMI files are much simpler than the new HCOOH file) and writing the HCOOH loss/prod/species variables into it. Would it be possible to upload this as an alternative to the 20170103 file? I’ve successfully tested this one with GCHP, and all the other GMI prod/loss files also work with GCHP.

Therefore, we have created a new file called gmi.clim.HCOOH.geos5.2x25.20170108.nc which has the same data as in gmi.clim.HCOOH.geos5.2x25.20160912.nc, but with the correct time stamps, and the extraneous dimensions removed.

We have also updated the relevant HEMCO_Config.template* files in the GEOS-Chem Unit Tester. We deleted the text in RED:

* GMI_LOSS_HCOOH      $ROOT/GMI/v2015-02/gmi.clim.HCOOH.geos5.2x25.20160912.nc loss  ...
* GMI_PROD_HCOOH      $ROOT/GMI/v2015-02/gmi.clim.HCOOH.geos5.2x25.20160912.nc prod  ...

and added the text in GREEN:

* GMI_LOSS_HCOOH      $ROOT/GMI/v2015-02/gmi.clim.HCOOH.geos5.2x25.20170108.nc loss  ...
* GMI_PROD_HCOOH      $ROOT/GMI/v2015-02/gmi.clim.HCOOH.geos5.2x25.20170108.nc prod  ...

--Bob Yantosca (talk) 17:20, 9 January 2017 (UTC)

Fix for monthly stratospheric P/L rates in HEMCO

This update was validated with the 1-month benchmark simulation v11-01j and approved on 03 Dec 2016

Brian Boys wrote:

It appears that the current public release of GC is not updating monthly stratospheric P/L rates. I've read through strat_chem_mod and it says this is supposed to happen through HEMCO (there's no call to GET_MONTH in strat_chem). I suspect this b/c I've noticed a lack of stratospheric NO2 seasonality in year long runs, with stratospheric NO2 columns reflecting the simulation start month seasonality. Is it possible that the stratospheric P/L rates are set at the simulation start month and not updated?

Christoph Keller responded:

Thanks for reporting this. Does your HEMCO_Config.rc have the following entry:

  * GMI_LOSS_NOx        $ROOT/GMI/v2015-02/gmi.clim.NOx.geos5.2x25.nc         loss    2000/$MM/1/0 C xyz s-1   NOx      - 1  1
  * GMI_PROD_NOx        $ROOT/GMI/v2015-02/gmi.clim.NOx.geos5.2x25.nc         prod    2000/$MM/1/0 C xyz v/v/s NOx      - 1  1

The GMI_PROD/LOSS entries should all be:

  * GMI_LOSS_NOx        $ROOT/GMI/v2015-02/gmi.clim.NOx.geos5.2x25.nc         loss    2000/1-12/1/0 C xyz s-1   NOx      - 1  1
  * GMI_PROD_NOx        $ROOT/GMI/v2015-02/gmi.clim.NOx.geos5.2x25.nc         prod    2000/1-12/1/0 C xyz v/v/s NOx      - 1  1

The STRAT_OH entry is (also) an error and should be fixed in the HEMCO configuration file.

So the entry should be

  * STRAT_OH            $ROOT/GMI/v2015-02/gmi.clim.OH.geos5.2x25.nc          species 2000/1-12/1/0 C xyz v/v   *        - 1  1

Instead of

  * STRAT_OH            $ROOT/GMI/v2015-02/gmi.clim.OH.geos5.2x25.nc          species 2000/$MM/1/0 C xyz v/v   *        - 1  1

The difference between these two entries is that the current entry (2000/$MM/1/0) reads the data only once (at the beginning) and then never updates it, whereas the correct entry (2000/1-12/1/0) updates the data every month – which is what we want.

--Melissa Sulprizio (talk) 13:57, 18 July 2016 (UTC)

Fixed parallelization errors in DO_STRAT_CHEM

This fix was included in GEOS-Chem v11-01g (approved 28 Sep 2016).

We discovered a parallelization error in routine DO_STRAT_CHEM (in module GeosCore/strat_chem_mod.F90), in the section where the Br_y concentrations are set:

       !===============================
       ! Prescribe Br_y concentrations
       !===============================
 
       !$OMP PARALLEL DO &
       !$OMP DEFAULT( SHARED ) &
       !$OMP PRIVATE( NN, BEFORE, I, J, L, BryTmp ) &
       !$OMP PRIVATE( LCYCLE )
       DO NN=1,6

          IF ( GC_Bry_TrID(NN) > 0 ) THEN

             ! Make note of inital state for determining tendency later
             ! NOTE: BEFORE has to be made PRIVATE to the DO loop since
             ! it only has IJL scope, but the loop is over IJLN!
             ! (bmy, 8/7/12)
             BEFORE = STT(:,:,:,GC_Bry_TrID(NN))
 
             ! Is this Br2?
             ISBR2 = ( TRIM(Input_Opt%TRACER_NAME(Strat_TrID_GC(NN))) == 'Br2' )

Actually there are two errors:

1. The ISBR2 variable (highlighted in RED) is a scalar that gets assigned a value within a parallel loop. Therefore, it should be declared as PRIVATE to the loop.
2. The test to determine if the species is Br2 used an incorrect algorithm. It returned FALSE even when the current species was indeed Br2.

The correct code should be:

      !--------------------------------------------------------------------
      ! Prescribe Br_y concentrations
      !--------------------------------------------------------------------

      !$OMP PARALLEL DO                                           &
      !$OMP DEFAULT( SHARED                                     ) &
      !$OMP PRIVATE( NN, BEFORE, ISBR2, L, J, I, LCYCLE, BryTmp )
      DO NN = 1,6

         IF ( GC_Bry_TrID(NN) > 0 ) THEN

            ! Make note of inital state for determining tendency later
            ! NOTE: BEFORE has to be made PRIVATE to the DO loop since
            ! it only has IJL scope, but the loop is over IJLN!
            ! (bmy, 8/7/12)
            BEFORE = Spc(:,:,:,GC_Bry_TrID(NN))

            ! Is this Br2?
            ISBR2  = ( Gc_Bry_TrId(NN) == id_Br2 )

We have now placed the ISBR2 in the !$OMP PRIVATE statement. We also now test the Gc_Bry_TrId array against the id_Br2 variable, which is computed as id_Br2 = Ind_('Br2') (see our Species indexing in GEOS-Chem wiki page for more information).

--Bob Yantosca (talk) 18:42, 3 January 2017 (UTC)

Bug fix in routine GET_RATES_INTERP

This update was validated in the 1-month benchmark simulation v10-01c and approved on 29 May 2014.

Routine GET_RATES_INTERP has now been removed from GEOS-Chem v10-01 and newer versions. GEOS-Chem's stratospheric chemistry module now relies upon the HEMCO emissions component to read production and loss rates that have been archived from the GMI model.

We have noted a bug in routine GET_RATES_INTERP (in module GeosCore/strat_chem_mod.F90). We fixed it as follows (by commenting out the offending line):

       ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR if necessary
       ptr_3d => LOSS(:,:,:,N)
!-----------------------------------------------
! Prior to 4/3/14:
! Array2 should be ptr_3d (bmy, 4/3/14)
!       call transfer_3D( array, array2 )
!-----------------------------------------------
       call transfer_3D( array, ptr_3D )
       NULLIFY( ptr_3D )

Routine GET_RATES_INTERP is called for nested-grid full-chemistry simulations, but not for global simulations.

--Bob Y. 17:10, 30 May 2014 (EDT)

Reduce memory footprint of the stratospheric chemistry module

This update was validated in the 1-month benchmark simulation v10-01c and approved on 29 May 2014.

Routine GET_RATES has now been removed from GEOS-Chem v10-01 and newer versions. GEOS-Chem's stratospheric chemistry module now relies upon the HEMCO emissions component to read production and loss rates that have been archived from the GMI model.

We have made a modification that should reduce the amount of memory that the stratospheric chemistry module requires, especially for simulations at 2° x 2.5° or higher horizontal resolution.

Bob Yantosca wrote:

When I ran the 2° x 2.5° standard code (compiled with MET=geosfp GRID=2x25) in our regular 8-CPU computational queue, the run died with a segmentation fault error. This error usually indicates that the code is using more than the available amount of memory. The seg fault occurred in the transport routine TPCORE_FVDAS, so for a while I thought the issue was there. But when I turned off the stratospheric chemistry, the run proceeded to completion. So this meant that the problem was in the strat_chem_mod.F90 module.

Long story short: what was happening is that the code was allocating BIG arrays in strat_chem_mod.F90 at the start of the simulation. This meant that there wasn't enough memory leftover by the time you arrived at the TPCORE transport routines. To fix this problem, I declared the following arrays as REAL*4 instead of REAL*8:

   !------------------------------------------------------------------------------
   ! Prior to 4/3/14:
   ! Change arrays from REAL*8 to REAL*4 to reduce memory footprint (bmy, 4/3/14)
   !  REAL*8,  ALLOCATABLE, TARGET :: PROD(:,:,:,:)   ! Production rate [v/v/s]
   !  REAL*8,  ALLOCATABLE, TARGET :: LOSS(:,:,:,:)   ! Loss frequency [s-1]
   !  REAL*8,  ALLOCATABLE, TARGET :: STRAT_OH(:,:,:) ! Monthly mean OH [v/v]
   !------------------------------------------------------------------------------
     REAL*4,  ALLOCATABLE, TARGET :: PROD(:,:,:,:)   ! Production rate [v/v/s]
     REAL*4,  ALLOCATABLE, TARGET :: LOSS(:,:,:,:)   ! Loss frequency [s-1]
     REAL*4,  ALLOCATABLE, TARGET :: STRAT_OH(:,:,:) ! Monthly mean OH [v/v]

     . . .

   !----------------------------------------------------------------------------
   ! Prior to 4/3/14:
   ! Change arrays from REAL*8 to REAL*4 to reduce memory footprint (bmy, 4/3/14)
   !  REAL*8, ALLOCATABLE  :: MInit(:,:,:,:)      ! Init. atm. state for STE period
   !  REAL*8, ALLOCATABLE  :: SChem_Tend(:,:,:,:) ! Stratospheric chemical tendency
   !----------------------------------------------------------------------------
     REAL*4, ALLOCATABLE  :: MInit(:,:,:,:)      ! Init. atm. state for STE period
     REAL*4, ALLOCATABLE  :: SChem_Tend(:,:,:,:) ! Stratospheric chemical tendency
                                                 !   (total P - L) [kg period-1]

Then I had to make the correesponding changes in module routines GET_RATES:

   !------------------------------------------------------------------------------
   ! Prior to 4/3/14:
   ! We no longer need this (bmy, 4/3/14)
   !    REAL*8             :: ARRAY2( IIPAR, JJPAR, LLPAR )  ! Actual vertical res
   !------------------------------------------------------------------------------

       . . .
  
       ! Pointers
       REAL*4,  POINTER   :: ptr_3D(:,:,:)

       . . .
 
       ! Intialize arrays
   !------------------------------------------------------------------------------
   ! Prior to 4/3/14:
   ! Initialize w/ 0e0 for REAL*4 (bmy, 4/3/14)
   !    LOSS = 0d0
   !    PROD = 0d0
   !------------------------------------------------------------------------------
       LOSS = 0e0
       PROD = 0e0

       . . .

   !-----------------------------------------------------------------------------
   ! Prior to 4/3/14:
   ! PROD is now REAL*4 (bmy, 4/13/14)
   !       ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR
   !       call transfer_3D( array, array2 )
   !
   !       PROD(:,:,:,N) = ARRAY2
   !-----------------------------------------------------------------------------
          ! Resize the data to LLPAR levels (if necessary)
          ptr_3D => PROD(:,:,:,N)
          CALL TRANSFER_3D( array, ptr_3D )
          NULLIFY( ptr_3D )

       . . .

   !-----------------------------------------------------------------------------
   ! Prior to 4/3/14:
   ! LOSS is now REAL*4 (bmy, 4/13/14)
   !       ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR
   !       call transfer_3D( array, array2 )
   !
   !       LOSS(:,:,:,N) = ARRAY2
   !-----------------------------------------------------------------------------
          ! Resize the data to LLPAR levels (if necessary)
          ptr_3D => LOSS(:,:,:,N)
          CALL TRANSFER_3D( array, ptr_3D )
          NULLIFY( ptr_3D )

Then I had to make the modifications listed above also to moduleroutine GET_RATES_INTERP.

Furthermore, I had to modify routine GeosUtil/transfer_mod.F in the following way:

(1) Rename routine TRANSFER_3D to TRANSFER_3D_R8

(2) Copy TRANSFER_3D_R8 to a new routine TRANSFER_3D_R4, anc declare the OUT argument as REAL*4 instead of REAL*8:

       !
       ! !OUTPUT PARAMETERS:
       !
             REAL*4,  INTENT(OUT) :: OUT(IIPAR,JJPAR,LLPAR)   ! Output data

(3) Add the following code at the top of the GeosUtil/transfer_mod.F:

         INTERFACE TRANSFER_3D
            MODULE PROCEDURE TRANSFER_3D_R4
            MODULE PROCEDURE TRANSFER_3D_R8
         END INTERFACE

           . . .

         PRIVATE :: TRANSFER_3D_R4
         PRIVATE :: TRANSFER_3D_R8

The modifications in GeosUtil/transfer_mod.F will allow the 72-level REAL*4 data to be resized and stored in the REAL*4 arrays in GeosCore/strat_chem_mod.F90.

All of these modifications will reduce the memory footprint of 2° x 2.5° (and higher-resolution) full-chemistry simulations. This was sufficient for me to run the 2° x 2.5° standard simulation in our normal 8-CPU computational queue. But if your simulations have a lot of diagnostics turned on, your code may still hit a hard memory limit.

--Bob Y. 17:10, 30 May 2014 (EDT)

Eliminate floating-point invalid in INIT_STRAT_CHEM

This fix was included GEOS-Chem v9-02o (Approved 03 Sep 2013).

In routine INIT_STRAT_CHEM (in module GeosCore/strat_chem_mod.F90), we need to initialize the REAL*8 variable TPAUSEL_CNT as follows:

    Old code:
    TpauseL_Cnt              = 0.

    New code:
    TpauseL_Cnt              = 0.d0

The old code causes a floating point invalid situation when compiling with FPE=yes.

--Bob Y. 12:11, 13 September 2013 (EDT)

Minor bug in parallel DO loop for Bry concentrations

A minor parallelization bug existed in the benchmark v9-01-03o. We have now fixed this in benchmark v9-01-03p.

The BEFORE variable was declared as a 3-D allocatable array at the top of module strat_chem_mod.F90:

 REAL*8, ALLOCATABLE :: Before(:,:,:)      ! Init atm. state each chem. dt

and was allocated at the start of the GEOS-Chem simulation by routine INIT_STRAT_CHEM. BEFORE was also used in this parallel loop in routine DO_STRAT_CHEM:

       !===============================
       ! Prescribe Br_y concentrations
       !===============================

       !$OMP PARALLEL DO &
       !$OMP DEFAULT( SHARED ) &
       !$OMP PRIVATE( NN, I, J, L, BryDay, BryNight, IJWINDOW )
       DO NN=1,6

          IF ( GC_Bry_TrID(NN) > 0 ) THEN

             ! Make note of inital state for determining tendency later
             BEFORE = STT(:,:,:,GC_Bry_TrID(NN))

             ... etc ...

As you can see, the BEFORE array has (I,J,L) scope, but is located within an OpenMP parallel DO loop with (I,J,L,N) scope. Therefore, BEFORE should have been declared in the !$OMP PRIVATE listing, but was not.

To correct this situation we changed BEFORE from an allocatable array to a local array within routine DO_STRAT_CHEM:

   ! Arrays
   REAL*8                    :: STT0(IIPAR,JJPAR,LLPAR,N_TRACERS)
   REAL*8                    :: BEFORE(IIPAR,JJPAR,LLPAR)

and also added it to the !OMP PRIVATE declaration in this parallel DO loop:

       !===============================
       ! Prescribe Br_y concentrations
       !===============================

       !$OMP PARALLEL DO &
       !$OMP DEFAULT( SHARED ) &
       !$OMP PRIVATE( NN, BEFORE, I, J, L, BryDay, BryNight, IJWINDOW )
       DO NN=1,6

          IF ( GC_Bry_TrID(NN) > 0 ) THEN

             ! Make note of inital state for determining tendency later
             ! NOTE: BEFORE has to be made PRIVATE to the DO loop since
             ! it only has IJL scope, but the loop is over IJLN!
             ! (bmy, 8/7/12)
             BEFORE = STT(:,:,:,GC_Bry_TrID(NN))

             ... etc ...

This ensures that the strat-trop exchange diagnostic will be computed properly.

--Bob Y. 17:29, 8 August 2012 (EDT)

Bug fix for Intel Fortran Compiler 12

Prasad Kasibhatla wrote:

I am experiencing the IFORT 12 compile failure of GEOS-Chem v9-01-03 with nested grid options turned on in define.h. The compile seems to be failing in the Stratospheric chemistry module GeosCore/strat_chem_mod.F90.

Please see this wiki post on our Intel Fortran Compiler wiki page for a full description of the issue.

--Bob Y. 15:27, 20 December 2012 (EST)