GEOS-Chem v8-03-01

From Geos-chem
Revision as of 23:01, 17 February 2011 by Ccarouge (Talk | contribs) (Outstanding issues not yet resolved in v8-03-01)

Jump to: navigation, search

NOTE: GEOS-Chem v8-03-01 refers to the version previously noted as GEOS-Chem v8-02-05. The updates in this version (e.g. new SOA tracers, new isoprene option, switch to Git version control, etc.) necessitated a full rewrite of the GEOS-Chem User's Manual.

Therefore, the GEOS-Chem Support Team has decided to make this version a public release and change the version number accordingly to GEOS-Chem v8-03-01.


PUBLIC RELEASE -- 04 May 2010

Will contain everything in GEOS-Chem v8-02-04, plus:

  1. Option to use Caltech isoprene chemistry (F. Paulot)
  2. ISORROPIA II package for aerosol thermodynamical equilibrium (H. Pye, T. Nenes)
  3. Updated aerosol optical properties for FAST-J photolysis (R. Martin, C. Heald)
  4. Modifications to SOA formation (Aerosols Working Group)
  5. SOA formation from aromatics (D. Henze)
  6. Speciated biogenic emissions from MEGAN v2.1 now used in SOA code (H. O. T. Pye)
  7. Global 1 x 1.25 simulation capability (L. Lamsal)
  8. Modifications to / Creation of the GEOS-5 0.5 x 0.667 N. American and European Nested Grid simulation (A. van Donkelaar)
  9. TOMAS aerosol microphysics option (W. Trivitayanurak, D. Westervelt, J. Pierce)
  10. Extension of annual anthropogenic scale factors to 2006 (A. van Donkelaar)
  11. EMEP emissions extended to 2007 and Seasonality extended to SOx, CO, and NH3 (A. van Donkelaar)
  12. 1st GEOS-Chem version to use the Git version control system

Previous issues now resolved in v8-03-01

Photolysis now computed at midpoint of chemistry timestep

Daniel Jacob wrote:

It was brought to my attention by Lin and Mauricio that GEOS-Chem calculates the Fast-J radiative environment for the beginning of the chemistry time step. This is wrong and produces a small but systematic bias. It should be calculated for the middle of the time step. So for example if you have a chemistry time step from 10 to 11 am the radiative environment should be calculated for 10:30.

Bob Yantosca wrote:

I've modified the code in v8-03-01 so that it computes the cosine of the solar zenith angle (aka SUNCOS) at the midpoint of the chemistry timestep. SUNCOS is then fed to FAST-J, so that it computes the photolysis also at the midpoint of the chemistry timestep.
Here is a plot of the JNO2 diurnal variation at the Harvard forest (70W, 42N) grid box. The black line is the J-value computed at the start of the chem timestep (i.e. the current way) and the red line is the J-value computed at the midpoint (i.e. what we are changing to).
As you can see the J-value curve computed at the midpt of the chemistry timestep starts peaking about a half-hour (which is 1/2 of the 60-min chem timestep) earlier, as you would expect.

--Bob Y. 16:27, 27 April 2010 (EDT)

Bug fix in ISORROPIA for offline aerosol

Eric Sofen wrote:

I've found a small bug in the new version of isoropiaII_mod.f. This bug doesn't appear in a full-chem simulation. In the offline aerosol sim, IDTHNO3 is not defined so that DEN_SAV cannot be calculated (around line 500 of isoropiaII_mod.f) and the run crashes.

Bob Yantosca wrote:

I installed your fix into ISORROPIA, it will be added into GEOS-Chem v8-03-01. You need to declare HNO3_DEN as a REAL*8 variable and put it in the !$OMP+PRIVATE statement:
   !$OMP+PRIVATE( HPLUSTEMP,    NUM_SAV, DEN_SAV, HNO3_DEN              ) 
Then add HNO3_DEN to the following IF block:
            ! Special handling for HNO3 [kg]
            IF ( IDTHNO3 > 0 ) THEN

               ! COUPLED SIMULATION

               ! HNO3 [mole/m3] is in GAS(2); convert & store in STT [kg]
               STT(I,J,L,IDTHNO3) = MAX( 63.d-3 * VOL * GAS(2), CONMIN )

               ! Save for use in DEN_SAV expression below (sofen, 4/21/10)
               HNO3_DEN           = STT(I,J,L,IDTHNO3 )


               ! OFFLINE SIMULATION:

               ! Convert total inorganic nitrate from [mole/m3] to [ug/m3]
               ! and save for next time
               ! WT(4) is in [mole/m3] -- unit conv is necessary!
               CALL SET_HNO3( I, J, L, 63.d6 * WT(4) )

               ! Save for use in sulfate_mod (SEASALT_CHEM) for offline
               ! aerosol simulations (bec, 4/15/05)
               GAS_HNO3(I,J,L) = GAS(2)

               ! Save for use in DEN_SAV expression below (sofen, 4/21/10)
               HNO3_DEN        = GAS(2) * VOL * 63d-3

And in the expression to DEN_SAV below:
            DEN_SAV            = ( STT(I,J,L,IDTSO4)  / 96d0   * 2d0     +
   &                               STT(I,J,L,IDTNIT)  / 62d0             +
   &                               HNO3_DEN           / 63d0             +
   &                               STT(I,J,L,IDTSALA) * 0.55d0 / 35.45d0 ) 

--Bob Y. 11:11, 21 April 2010 (EDT)

Minor bug fix in SRCSO2

Eric Leibensperger wrote:

I think that there is possibly a small bug in sulfate_mod.f relating to Streets biofuel emissions. Streets' emissions include biofuel emissions and so the default global emissions should be zeroed out. This is done for EDGAR, but not GEIA. EDGAR is the default for SO2 and I don't think anyone uses GEIA emissions anymore, but it's good to be consistent.

The following line needs to be added in routine SRCSO2 of sulfate_mod.f:


        ! Use GEIA SOx emissions


              ! If we are using David Streets' emissions
              ! Remember: those include BF if Year is GE 2006
              IF ( LSTREETS ) THEN

                 ! If we are over the SE Asia region
                 IF ( GET_SE_ASIA_MASK( I, J ) > 0d0 ) THEN

                    ! Overwrite GEIA SO2 w/ David Streets' SO2 [kg SO2/s]
                    ESO2_an(I,J,1) = GET_STREETS_ANTHRO( I, J, IDTSO2, 
    &                                                    KG_S=.TRUE. )

                    ! Zero 2nd level of emissions
                    ESO2_an(I,J,2) = 0d0

                    !%%%%%% ADD THESE LINES! %%%%%%
                    ! Streets 2006 includes biofuels.
                    IF ( SCALEYEAR >= 2006 ) ESO2_bf(I,J) = 0d0


--Bob Y. 16:21, 15 March 2010 (EDT)


The software patch for the bug in GET_EMMONOT_MEGAN from v8-02-04 has now been standardized in this version.

--Bob Y. 16:33, 9 March 2010 (EST)


The software patch for the bug in READ_ANTHRO_NH3 from v8-02-04 has now been standardized in this version.

--Bob Y. 10:47, 5 March 2010 (EST)

Fix for initialization of EMEP ship emissions

Helen Macintyre wrote:

I am trying to do some simulations with 1985 emissions, but I find that I run into an SST error.
Looking through the log file reveals this:
   M O N T H L Y   E M E P   E U R O P E A N     E M I S S I O N S
   ( INCL. SHIP )

   Base Year :1985

   EMEP anthropogenic NOx  for 1985/07:     -Infinity Tg N
   EMEP anthropogenic CO   for 1985/07:           NaN Tg
   EMEP anthropogenic SO2  for 1985/07:      2.217651 Tg S
   EMEP anthropogenic NH3  for 1985/07:      0.771631 Tg

Claire Carouge wrote:

It looks like a bug in initialisation of EMEP ship emissions for CO and NOx.
In INIT_EMEP (in emep_mod.f) there are the lines:
     IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' )
     EMEP_SO2_SHIP = 0d0
     IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' )
     EMEP_SO2_SHIP = 0d0
     IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' )
     EMEP_SO2_SHIP = 0d0
So after allocating EMEP_CO_SHIP or EMEP_NOx_SHIP we initialize EMEP_SO2_SHIP to 0d0. A typical problem of copy/paste.
You should change the lines to initialize the appropriate ship arrays to 0d0:
     IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' )
     EMEP_SO2_SHIP = 0d0
     IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_CO_SHIP' )
     EMEP_CO_SHIP = 0d0
     IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_NOx_SHIP' )
     EMEP_NOx_SHIP = 0d0

--Ccarouge 12:54, 26 January 2010 (EST)

Bug fix in AVGPOLE

Mike Barkley wrote:

Whilst constructing my Amazon grid I noticed that in the sub-routine AVGPOLE there are the following lines of code (lines 326-332):
   #if   defined( GRID1x1 )
   #if   defined( NESTED_CH ) || defined( NESTED_NA )
        ! NOTE: Only do this for 1x1 nested grids (bmy, 12/1/04)
        ! 1x1 window grid does not extend to poles
If someone runs a nested 0.5x0.666 simulation then the adjustment for the 'poles' will still be carried out (which will actually be at the edges of the nested window); it should be modified to:
   #if   defined( GRID1x1 ) || defined( GRID05x0666 )

This has now been fixed in v8-02-05.

--Bob Y. 10:50, 18 December 2009 (EST)

Bug fix in EMITHIGH (carbon_mod.f)

Bob Yantosca wrote:

Upon further inspection in routine EMITHIGH (in carbon_mod.f), it was discovered that the variable EMIS_SAVE had incorrect indices for the LIMO and ALCO tracers. EMIS_SAVE(:,:,...) should be EMIS_SAVE(I,J,...) in both cases.
Therefore, the following fix needs to be applied to the DO loop below.
         DO J = 1, JJPAR
         DO I = 1, IIPAR
            IF ( IS_BCPI ) EMIS_SAVE(I,J,IDTBCPI) = BCSRC(I,J,1)
            IF ( IS_OCPI ) EMIS_SAVE(I,J,IDTOCPI) = OCSRC(I,J,1)
            IF ( IS_BCPO ) EMIS_SAVE(I,J,IDTBCPO) = BCSRC(I,J,2)
            IF ( IS_OCPO ) EMIS_SAVE(I,J,IDTOCPO) = OCSRC(I,J,2)
            IF ( IS_LIMO ) THEN
               ! Prior to 1/11/09:
               ! Bug fix: Should be EMIS_SAVE(I,J,IDTLIMO) since we are
               ! within a parallel loop.
               !EMIS_SAVE(:,:,IDTLIMO) = BIOG_LIMO(I,J)
               EMIS_SAVE(I,J,IDTLIMO) = BIOG_LIMO(I,J)
               ! lead to too much ORVC_TERP in the 1st layer?
               ORVC_TERP(I,J,1)   = ORVC_TERP(I,J,1) + BIOG_TERP(I,J)
            IF ( IS_ALCO ) THEN
               ! Prior to 1/11/09:
               ! Bug fix: Should be EMIS_SAVE(I,J,IDTALCO) since we are
               ! within a parallel loop.
               !EMIS_SAVE(:,:,IDTALCO) = BIOG_ALCO(I,J)
               EMIS_SAVE(I,J,IDTALCO) = BIOG_ALCO(I,J)

               ! lead to too much ORVC_SESQ in the 1st layer?
               ORVC_SESQ(I,J,1)   = ORVC_SESQ(I,J,1) + BIOG_SESQ(I,J)

NOTE: This error only affects the GEOS-Chem 54-tracer simulation if the non-local PBL mixing scheme is turned on. This error would therefore not show up in the 43-tracer benchmarks (both 1-month and 1-year).

--Bob Y. 15:18, 11 January 2010 (EST)

Fix to prevent numerical overflow in EXP function

Win Trivitayanurak wrote:

I'm reporting a bug fix, or maybe more appropriately a numerical error preventing fix, to sulfate_mod.f
The symptom i found was a full-chem run (in both v.8-01-04 and v.8-02-01) would crash in nested-grid (china) simulation precisely on 2008/04/24 01:00 with no proper error diagnostic in the log file. After tracking down the error, i found that there is an infinity occuring in sulfate_mod.f in SUBROUTINE CHEM_SO2 where SO2 loss via aqueous chemistry is calculated. Infinity occurred because there is a number in the EXP term that is too large to handle, i.e. >700. My solution was to put a check to this term entering EXP.
I'm not super certain whether there is something else wrong here, i.e. why would the so2_ss term be so large or H2O2 and/or O3 be so tiny. But this numerical check is sure to fix this simulation crash.

Bob Yantosca wrote:

I've implemented Win's exponential-check fix into the v8-02-05. I made the code a little more general. I created routines IS_SAFE_EXP and SAFE_EXP in error_mod.f. These routines allow you to check ahead if the exponential value would be too large before you actually do the exponential:
Thus the code becomes:

        . . .

        REAL*8 :: XX
        . . .

              ! Compute aqueous rxn rates for SO2
              CALL AQCHEM_SO2( LWC, TK,    PATM,    SO2_ss, H2O20,
       &                       O3,  HPLUS, KaqH2O2, KaqO3 )

   ! Prior to 1/4/10:
   !            ! Aqueous phase SO2 loss rate (v/v/timestep):
   !            L2  = EXP( ( SO2_ss - H2O20 ) * KaqH2O2 * DTCHEM ) 
   !            L3  = EXP( ( SO2_ss - O3    ) * KaqO3   * DTCHEM )      !
   !            ! Loss by H2O2
   !            L2S = SO2_ss * H2O20 * (L2 - 1.D0) / ((SO2_ss * L2) - H2O20) !
   !            ! Loss by O3
   !            L3S = SO2_ss * O3    * (L3 - 1.D0) / ((SO2_ss * L3) - O3)   

              ! Compute loss by H2O2.  Prevent floating-point exception
              ! by not allowing the exponential to go to infinity if
              ! the argument is too large.  (win, bmy, 1/4/09)

              ! Argument of the exponential
              XX  = ( SO2_ss - H2O20 ) * KaqH2O2 * DTCHEM

              ! Test if EXP(XX) can be computed w/o numerical exception
              IF ( IS_SAFE_EXP( XX ) ) THEN

                 ! Aqueous phase SO2 loss rate w/ H2O2 [v/v/timestep]
                 L2  = EXP( XX )

                 ! Loss by H2O2
                 L2S = SO2_ss * H2O20 * ( L2 - 1.D0 ) /
   &                   ( (SO2_ss * L2) - H2O20 )            ELSE

                 ! If exponential can't be computed, set to
                 ! the initial H2O2 concentration
                 L2S = H2O20


              ! Compute loss by O3.  Prevent floating-point exception
              ! by not allowing the exponential to go to infinity if
              ! the argument is too large. (win, bmy, 1/4/09)

              ! Argument of the exponential
              XX = ( SO2_ss - O3 ) * KaqO3 * DTCHEM

              ! Test if EXP(XX) can be computed w/o numerical exception
              IF ( IS_SAFE_EXP( XX ) ) THEN

                 ! Aqueous phase SO2 loss rate w/ O3 [v/v/timestep]
                 L3  = EXP( XX )

                 ! Loss by O3
                 L3S = SO2_ss * O3 * (L3 - 1.D0) / ((SO2_ss * L3) - O3)

                 ! If exponential can't be computed, set to
                 ! the initial O3 concentration
                 L3S = O3

              SO2_ss = MAX( SO2_ss - ( L2S + L3S ), MINDAT )
              H2O20  = MAX( H2O20  - L2S,           MINDAT )


--Bob Y. 16:05, 11 January 2010 (EST)

IMPORTANT: You have to declare XX as PRIVATE in the OMP derivative as well, otherwise your simulations will contain significant numerical errors and will not be reliable. This error took me two weeks to figure out -- Jintai Lin

Post-release patches

These patches have been added to the Git repository as of 30 Jun 2010. If you have already downloaded GEOS-Chem v8-03-01, you can get these patches by doing a git pull operation:

git pull git:// master

into your own code directory.

Patch to enable call to SOA_PARA_INIT

Joseph Enberg wrote:

I recently finished a year-long run using GEOS-Chem v8-03-01 with the SOA portions of the code turned on and ran into a problem. The concentrations of all SOA species never increased to any significant values, although the POA concentrations seemed to be reasonable. I explained the problem to Havala Pye and she quickly noticed that the code no longer called the subroutine SOA_PARA_INIT right before the main SOA_CHEMISTRY do loops. Without this, the values of the equilibrium constants and yield parameters remain at their initial values of zero and nothing happens. I implemented the call statement into the code from a previous version of GEOS-Chem and everything seems reasonable now.

--Bob Y. 12:56, 30 June 2010 (EDT)

Patch to fix declaration of IONIC in ISORROPIA II

Havala O. T. Pye wrote:

I noticed that the -r8 flag (ifort) used with isoropiaIIcode.f was causing surface nitrate levels to decrease 25-50% compared to compilation without that flag. I found a problem with a variable in the common blocks.
The IONIC variable was declared REAL in the original Isoropia and declared as REAL*4 in the GEOS-Chem implementation. When the -r8 flag was implemented, IONIC was not converted to REAL*8, and there was an inconsistency between the type of variable in and the use of the variable when it was passed to routines in isoropiaIIcode.f. Since GEOS-Chem is using the -r8 flag, I have declared IONIC to be REAL*8. The surface concentration (and global nitrate burden) are now more similar in magnitude and spatial distribution to RPMARES instead of being significantly lower.

--Bob Y. 12:56, 30 June 2010 (EDT)

Patch to fix SCHEM data input

Lee Murray wrote:

I've discovered a bug in how GEOS-5 does its stratospheric chemistry. This affects all G-5 users running full chemistry, tagged CO or methane simulations.
The 2D vertical prod and loss files in pco_lco_200203/ and other directories have the dimensions of the non-reduced vertical resolution for GEOS3 (LGLOB=48) and GEOS4 (LGLOB=55) and are fine. The GEOS5 files however are using the reduced grid resolution.
The functions in schem.f / tagged_co_mod.f / methane_mod.f are looking for arrays of vertical dimension 1:72, and then collapsing whatever it read to the reduced grid using the transfer_mod.f functions. (I'm not sure how a seg fault wasn't triggered for the December values that should have exceeded the EOF when reading the bpch.) I've attached a plot of the production of CO in January as internally seen in GEOS-Chem to demonstrate the problem.
It appears the stratOH_200203/ and stratjv_200203/ geos5 files have the same issue, though the pnoy_200106/ files are okay.
That being said, I am updating the various strat chem prod / loss files with the 3D fields from the GMI simulations, and will do so with the full 72 levels. But until that's available, you might want to send out a patch for GEOS-5 users.

Bob Yantosca wrote:

I created a software patch for the SCHEM fix. The patch is as follows:
  1. Don't allow people to use GEOS-5 with 72L
  2. Put #if defined(GEOS_5) blocks to prevent calling TRANSFER_ZONAL from the various routines (tagged_co_mod.f, global_ch4_mod.f, schem.f)

Please see this picture to view the before and after effect of this patch. We recommend obtaining this patch as soon as possible.

--Bob Y. 16:43, 2 June 2010 (EDT)

Patch to fix bug in Aromatic SOA

Colette Heald wrote

There's a small bug in the implementation of the aromatic SOA in v8-03-01. In aerosol_mod.f the SOA5 concentrations should be added into OCPI. Therefore ~line 395 the code should read:
              ! Hydrophilic primary OC plus SOA [kg/m3A.  We need
              ! to multiply by OCF to account for the mass of other
              ! components which are attached to the OC aerosol.
              ! (rjp, bmy, 7/15/04)
              ! (clh, 05/19/10) bug corrected SOA5 missing from OCPI
              OCPI(I,J,L) = ( STT(I,J,L,IDTOCPI) * OCF
    &                     +   STT(I,J,L,IDTSOA1)
    &                     +   STT(I,J,L,IDTSOA2)
    &                     +   STT(I,J,L,IDTSOA3)
    &                     +   STT(I,J,L,IDTSOA4)
    &                     +   STT(I,J,L,IDTSOA5) ) / AIRVOL(I,J,L)
You also need to add this line to the top of the routine:

--Bob Y. 09:53, 20 May 2010 (EDT)

Patch to scale AOD diagnostic output to other wavelengths

Bob Yantosca wrote

Dear GEOS-Chem Users,
Colette Heald pointed out to us that we had omitted an update for the aerosol optical depth diagnostics in v8-03-01. Therefore, we have just added this as a patch.
The aerosol optical depth diagnostic output (in ND21, the planeflight output, and in the timeseries diagnostics) can now be scaled from 400nm (i.e. the wavelength band in which they are computed for FAST-J) to any arbitrary wavelength. This will facilitate comparing AOD output from GEOS-Chem to various satellite instruments.
There is a new file jv_spec_aod.dat which contains information about the "target" wavelength to which the AOD output will be scaled. The default setting for "jv_spec_aod.dat" is 550nm. You can modify the file for other wavelenghts if so desired. More information can be found HERE.
We have added the jv_spec_aod.dat to all of the GEOS-Chem run directory repositories. You can download these with Git as well. Please see our Using Git with GEOS-Chem for the download commands.
Once again, this patch does not change any scientific results, but just contains updates for diagnostic options.

--Bob Y. 16:07, 28 May 2010 (EDT)

Outstanding issues not yet resolved in v8-03-01

Modified scale factors for anthropogenic emissions

This is more of a clarification than an outstanding issue. The new anthropogenic scaling factors used in GEOS-Chem v8-03-01 (in the anth_scale_factors_200911 subdirectory) contain minor differences from the data used in GEOS-Chem v8-02-04 (in the anth_scale_factors_200905 subdirectory). Please see this discussion for a full explanation.

--Bob Y. 14:22, 22 April 2010 (EDT)

Code slow down: parallelization problem in setemis.f

The implementation of a scaling of NOx emissions in GEOS-Chem version 8-03-01 resulted in a parallelization problem in setemis.f. In consequence, some simulations would run about 30% slower as expected. To solve the problem, modify the lines:

                         ! Save aircraft NOx [molec NO/cm3/s] in REMIS
                         REMIS(JLOOP,N) = REMIS(JLOOP,N) + EMIS_BL
                !FP 15/12/2009
                ! Add possible scaling of NOx emissions


                         ! Save aircraft NOx [molec NO/cm3/s] in REMIS
                         REMIS(JLOOP,N) = REMIS(JLOOP,N) + EMIS_BL
                      ! NOx scaling moved here for optimisation. (ccc, 11/10/10)
                      ! Add possible scaling of NOx emissions
 !--- Prior to (ccc, 11/10/10). Moved to previous loop over grid cells
 !               !FP 15/12/2009
 !               ! Add possible scaling of NOx emissions
 !               REMIS(:,N)=NOx_SCALING*REMIS(:,N)             

In short, you need to move the NOx scaling into the loop over grid cells. This problem was corrected in version 8-03-02 by issuing an after-release patch.

--Ccarouge 18:21, 11 November 2010 (EST)

Double counting biofuel emissions over Asia

Double counting of biofuel emissions over Asia happened when using Streets emissions for simulations of the years between 2001 and 2005. It will be fixed in the version v9-01-02

--Ccarouge 18:01, 17 February 2011 (EST)