Difference between revisions of "GEOS-Chem v8-02-04"

From Geos-chem
Jump to: navigation, search
(Run1)
(Run2)
Line 127: Line 127:
 
| Red
 
| Red
 
| v8-02-04 Run0
 
| v8-02-04 Run0
| "new" hybrid MEGAN<br>with AVHRR-derived LAI
+
| "new" hybrid MEGAN
 
| OFF
 
| OFF
 
| "old" restart file
 
| "old" restart file

Revision as of 21:15, 26 January 2010

Overview

BETA RELEASE -- Date TBA

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

Note about Linoz validation

Dylan Jones (dbj@atmosp.physics.utoronto.ca) wrote:

Testing [the Linoz code in v8-02-04] was more difficult that I thought. I began by trying to compare the output of v8-02-04 with our previous runs with v8-02-01. I accounted for the changes in the transport_mod.f and I tried to undo the changes in when the diagnostics are archived in v8-02-04, but I was still getting large differences between v8-02-04 and v8-02-01. I finally gave up on this since I may have made a mistake in reverting to the old way of doing the diagnostics in v8-02-04. In the end I took the new linoz code from v8-02-04 and used it in v8-02-01. I ran two GEOS-5 full chemistry simulations for 2007 and the output were consistent over the full year.
I think that it is safe to release [Linoz in v8-02-04]. However, we should acknowledge that it was [only] tested in v8-02-01, since I was not able to assess the quality of the output in v8-02-04.

--Bob Y. 16:52, 18 November 2009 (EST)

1-year benchmark simulations

Run0

Three GEOS-Chem model versions were compared to each other:

Color Quantity Plotted Anthro Emissions Biogenic Emissions Photolysis and Chemistry
Mechanism
PBL scheme Annual Mean OH
[105 molec/cm3]
Red v8-01-04 Run2 GEIA/Piccot
EDGAR emissions
EMEP European emissions
BRAVO Mexican emissions
David Streets 2006 emissions
CAC Canadian emissions
EPA/NEI99 with ICARTT fix
EDGAR ship emissions
ARCTAS ship SO2 emissions
Anthro scale year 2005
"old" MEGAN "old" jv_spec.dat.
Chemistry mostly unchanged
from v5-07-08
except for
a few revisions
Fully mixed scheme 11.099
Green v8-02-01 Run0 " " " " "updated" jv_spec.dat
"updated" chemical mechanism
(cf. J. Mao, D. Millet,
T-M. Fu, Palmer Group
@ U. Edinburgh)
Fully mixed scheme 11.812
Blue v8-02-04 Run0 GEIA/Piccot
EDGAR emissions
EMEP European emissions
BRAVO Mexican emissions
David Streets 2006 emissions
CAC Canadian emissions
NEI 2005
EDGAR ship emissions
EMEP ship emissions
ARCTAS ship SO2 emissions
Anthro scale year 2005
"new" hybrid MEGAN " " Non-Local scheme 10.450
Black Observations          

Run1

Three GEOS-Chem model versions were compared to each other:

Color Quantity Plotted Anthro Emissions Biogenic Emissions PBL scheme Annual Mean OH
[105 molec/cm3]
Red v8-02-01 Run0 GEIA/Piccot
EDGAR emissions
EMEP European emissions
BRAVO Mexican emissions
David Streets 2006 emissions
CAC Canadian emissions
EPA/NEI99 with ICARTT fix
EDGAR ship emissions
ARCTAS ship SO2 emissions
Anthro scale year 2005
"old" MEGAN Fully mixed scheme 11.812
Green v8-02-04 Run0 GEIA/Piccot
EDGAR emissions
EMEP European emissions
BRAVO Mexican emissions
David Streets 2006 emissions
CAC Canadian emissions
NEI 2005
EDGAR ship emissions
EMEP ship emissions
ARCTAS ship SO2 emissions
Anthro scale year 2005
"new" hybrid MEGAN Non-Local scheme 10.450
Blue v8-02-04 Run1 " " "new" hybrid MEGAN with
MODIS-derived LAI
" " 11.048
Black Observations        

Run2

Three GEOS-Chem model versions were compared to each other:

Color Quantity Plotted Biogenic Emissions Linoz Restart file Annual Mean OH
[105 molec/cm3]
Red v8-02-04 Run0 "new" hybrid MEGAN OFF "old" restart file 10.450
Green v8-02-04 Run1 "new" hybrid MEGAN with
MODIS-derived LAI
OFF " " 11.048
Blue v8-02-04 Run2 " " ON Restart file with O3
concentrations from an
offline Ox simulation
spin-up with Linoz
11.067
Black Observations        

Previous issues now resolved in v8-02-04

Memory issue with KPP

The interface between GEOS-Chem and KPP was modified. The memory load of a GEOS-Chem simulation with KPP is now only slightly higher than with SMVGEAR. The results are not affected at all by the change.

Please refers to the interfacing wiki page if you need to re-build KPP for GEOS-Chem version 8-02-04 and higher.

--Ccarouge 14:30, 10 December 2009 (EST)

Updated reactions in sulfate_mod.f

Jenny Fisher (jafisher@fas.harvard.edu) wrote:

I made two very small changes to sulfate_mod.f that should be included in the standard code. The change updates the reaction rates used for the offline aerosol simulation so that they match the rates Jingqiu implemented in the full chemistry simulation.
The changes are:
  1. In CHEM_DMS, definition of RK2.
  2. In CHEM_SO2, definitions of K0 and Ki. Note that I also removed the definition of Ki earlier in the subroutine (line 1497) and moved it to this section so it can be easily updated for future rate constant revisions.

--Bob Y. 09:50, 15 October 2009 (EDT)

Bug fix for EMEP ship emissions

Monika Kopacz (mkopacz@princeton.edu) wrote:

I noticed that the code crashes when EMEP is set to false but ship emissions from EMEP are true in the input file. Sure, it makes sense, but I would suggest the following addition to the input_mod.f, around say around line 1550 to prevent the code from crashing.
   ! Make sure we're not using EMEP ship emissions, if EMEP is off
   IF (.not. LEMEP ) THEN
      LEMEPSHIP = .FALSE.
   ENDIF 

Bob Yantosca (yantosca@seas.harvard.edu) replied:

Good catch. We've put this code into READ_EMISSIONS_MENU in input_mod.f:
    !%%% Bug fix!  If LEMEPSHIP is turned on but LEMEP is turned %%%
    !%%% off, this will cause an error (because arrays are not   %%%
    !%%% allocated, etc.).  For now, just turn off LEMEPSHIP     %%%
    !%%% and print a warning message.  Whoever wants to fix this %%%
    !%%% in a more robust way is welcome to do so.               %%%
    !%%% (mak, bmy, 10/19/09)                                    %%%
    IF ( LEMEPSHIP .and. ( .not. LEMEP ) ) THEN
       LEMEPSHIP = .FALSE.
       WRITE( 6, '(a)' ) REPEAT( '=', 79 )
       WRITE( 6, '(a)' ) 'WARNING! EMEP emissions are turned off,'
       WRITE( 6, '(a)' ) 'so also turn off EMEP ship emissions'
       WRITE( 6, '(a)' ) 'in order to avoid crashing GEOS-Chem!'
       WRITE( 6, '(a)' ) REPEAT( '=', 79 )
    ENDIF
As you pointed out the simple solution is to just turn off LEMEPSHIP if LEMEP=F. However, it may still be possible to keep LEMEPSHIP turned on...one would have to go thru emep_mod.f to make sure that the ship emissions arrays get allocated even if the other arrays don't. Whoever wants to can feel free to do that...

--Bob Y. 10:57, 19 October 2009 (EDT)

Minor fix for molecular weights in drydep_mod.f

Bob Yantosca (yantosca@seas.harvard.edu) replied:

I found a minor inconsistency in INIT_DRYDEP (in drydep_mod.f) which only affects the SOA tracers. In the big loop where we initialize everything, the following updates were made:
       ! ALPH (Alpha-pinene)
       ELSE IF ( N == IDTALPH ) THEN
          ...

          !----------------------------------------------------------
          ! Prior to 10/19/09:
          ! Molwt should be 136.23 not 136 even (bmy, 10/19/09)
          !XMW(NUMDEP)     = 136d-3
          !----------------------------------------------------------
          XMW(NUMDEP)     = 136.23d-3

          ...            

       ! LIMO (Limonene)
       ELSE IF ( N == IDTLIMO ) THEN
          ...

          !----------------------------------------------------------
          ! Prior to 10/19/09:
          ! Molwt should be 136.23 not 136 even (bmy, 10/19/09)
          !XMW(NUMDEP)     = 136d-3
          !----------------------------------------------------------
          XMW(NUMDEP)     = 136.23d-3
  
So in other words, the mol wt of ALPH and LIMO was listed as 136d-3 but it should be 136.23d-3. It doesn't change things drastically, just in the last few decimal places of the surface resistance array RSURFC.

--Bob Y. 14:46, 19 October 2009 (EDT)

Bug fix in biomass_mod.f

Monika Kopacz (mkopacz@princeton.edu) and Jenny Fisher (jafisher@fas.harvard.edu) wrote:

We have a bug fix for biomass_mod.f. This is based on the GEOS-Chem v8-02-03 code, though we believe the problem has been there since at least GEOS-Chem v8-02-01. The corrected file is located via anonymous FTP at:
   ftp://ftp.gfdl.noaa.gov/home/m1k/biomass_mod.f.corrected
This fix required moving the co-emitted VOC scaling of emissions (in routine COMPUTE_BIOMASS_EMISSIONS) to places where the BIOMASS array is being reset or recalculated to prevent rescaling at every time step and thus blowing up. It might not have been a problem when using GFED2 emissions, but it was definitely a problem for other inventories that do not update the BIOMASS array at every time step.

This fix will be included in GEOS-Chem v8-02-04.

--Bob Y. 16:36, 6 November 2009 (EST)

Minor bug fix in gamap_mod.f

Daven Henze (Daven.Henze@Colorado.EDU) wrote:

The wiki says that MOLC was switched from 3 to 2 for ALD2 in gamap_mod.f, but the v8-02-03 code looks like the switch was accidentally applied to PRPE instead,
                 CASE( 6  )
                    NAME (T,28) = 'ALD2'
                    INDEX(T,28) = 11 + ( SPACING * 45 )
                    MWT  (T,28) = 12e-3
                    MOLC (T,28) = 3
                    UNIT (T,28) = 'atoms C/cm2/s'
                 CASE( 7  )
                    NAME (T,28) = 'PRPE'
                    INDEX(T,28) = 18 + ( SPACING * 45 )
                    MWT  (T,28) = 12e-3
                    !-------------------------------------------------
                    ! Prior to 10/9/09:
                    ! ALD2 has 2 carbons, not 3 (dbm, bmy, 10/9/09)
                    !MOLC (T,28) = 3
                    !-------------------------------------------------
                    MOLC (T,28) = 2
                    UNIT (T,28) = 'atoms C/cm2/s'

Bob Yantosca (yantosca@seas.harvard.edu) replied:

Yes it was a typo. I've since fixed it, it'll go into GEOS-Chem v8-02-04. Indeed, ALD2 should have 2 carbons and PRPE should have 3 carbons.

--Bob Y. 12:36, 19 November 2009 (EST)

Fixes and updates in seasalt_mod.f

Updates to the sea-salt emissions scheme

Lyatt Jaeglé (jaegle@atmos.washington.edu) wrote:

I recommend 3 fixes to the Monahan sea-salt emissions scheme in GEOS-Chem:
1. The source function is for wet aerosol radius (RH=80%, with a radius twice the size of dry aerosols) so in function SRCSALT in seasalt_mod.f, BETHA should be set to 2 instead of 1. Change the line:
     REAL*8, PARAMETER      :: BETHA = 1.d0
to:
     REAL*8, PARAMETER      :: BETHA = 2.d0
2. It turns out that the LOG in the Monahan source function is a log base 10. Therefore these lines:
           ! Sea salt base source [kg/m2]
           SRC(R,N)  = CONST * SS_DEN( N )
    &           * ( 1.d0 + 0.057d0*( BETHA * RMID(R,N) )**1.05d0 )
    &           * 10d0**( 1.19d0*
    &                  EXP(-((0.38d0-LOG(BETHA*RMID(R,N)))/0.65d0)**2))
    &           / BETHA**2

           ! Sea salt base source [#/m2] (bec, bmy, 4/13/05)
           SRC_N(R,N) = CONST_N * (1.d0/RMID(R,N)**3)
    &           * (1.d0+0.057d0*(BETHA*RMID(R,N))**1.05d0)
    &           * 10d0**(1.19d0*EXP(-((0.38d0-LOG(BETHA*RMID(R,N)))
    &           /0.65d0)**2))/ BETHA**2
should be replaced with
           !-----------------------------------------------------------
           ! IMPORTANT NOTE!
           !
           ! In mathematics, "LOG" means "log10". 
           ! In Fortran,     "LOG" means "ln" (and LOG10 is "log10").
           !
           ! The following equations require log to the base 10, so 
           ! we need to use the Fortran function LOG10 instead of LOG. 
           ! (jaegle, bmy, 11/23/09)
           !-----------------------------------------------------------

           ! Sea salt base source [kg/m2]
           SRC(R,N)  = CONST * SS_DEN( N )
    &           * ( 1.d0 + 0.057d0*( BETHA * RMID(R,N) )**1.05d0 )
    &           * 10d0**( 1.19d0*
    &                  EXP(-((0.38d0-LOG10(BETHA*RMID(R,N)))/0.65d0)**2))
    &           / BETHA**2

           ! Sea salt base source [#/m2] (bec, bmy, 4/13/05)
           SRC_N(R,N) = CONST_N * (1.d0/RMID(R,N)**3)
    &           * (1.d0+0.057d0*(BETHA*RMID(R,N))**1.05d0)
    &           * 10d0**(1.19d0*EXP(-((0.38d0-LOG10(BETHA*RMID(R,N)))
    &           /0.65d0)**2))/ BETHA**2

3. Finally I think that we should change the dry size bins for the coarse mode aerosols in input.geos. Please see the discussion on our Sea salt aerosols wiki page for more information.

--Bob Y. 10:44, 23 November 2009 (EST)

Updated hygroscopic growth factors

Becky Alexander (beckya@u.washington.edu) wrote:

The hygroscopic growth factors in routine GET_ALK (seasalt_mod.f) are incorrect. These:
       ! hygroscopic growth factor for sea-salt from Chin et al. (2002)
       IF ( IRH < 100 ) HGF = 2.2d0
       IF ( IRH < 99  ) HGF = 1.9d0
       IF ( IRH < 95  ) HGF = 1.8d0
       IF ( IRH < 90  ) HGF = 1.6d0
       IF ( IRH < 80  ) HGF = 1.5d0
       IF ( IRH < 70  ) HGF = 1.4d0
       IF ( IRH < 50  ) HGF = 1.0d0
should be replaced with:
       ! hygroscopic growth factor for sea-salt from Chin et al. (2002)
       IF ( IRH < 100 ) HGF = 4.8d0
       IF ( IRH < 99  ) HGF = 2.9d0
       IF ( IRH < 95  ) HGF = 2.4d0
       IF ( IRH < 90  ) HGF = 2.0d0
       IF ( IRH < 80  ) HGF = 1.8d0
       IF ( IRH < 70  ) HGF = 1.6d0
       IF ( IRH < 50  ) HGF = 1.0d0

This has now been done in GEOS-Chem v8-02-04.

--Bob Y. 09:48, 23 November 2009 (EST)

Bug fix in emfossil.f for 0.5 x 0.666 nested grid tagged-CO option

Yuxuan Wang (zjiang@atmosp.physics.utoronto.ca) wrote:

Zhe Jiang at University of Toronto found a bug in emfossil.f for the nested_CH tagged-CO option in code version v8-02-01 and higher. In emfossil.f, we need to make the following fix:
        IF ( ITS_A_TAGCO_SIM() ) THEN           
           !----------------------------------------------------------------
           ! Prior to 11/23/09:
           ! The reason for this change is because LNEI99 is set as .FALSE. 
           ! for the NESTED_CH option so GET_USA_MASK doesn't have a 
           ! definition which causes the run to crash.  So we need to
           ! add an ELSE to the IF ( LICARTT ) block to cover for this.
           ! (zie, yxw, bmy, 11/23/09)
           !IF ( LICARTT ) THEN
           !   IF ( GET_USA_MASK(I,J) > 0.d0 ) THEN
           !      IF ( NN == IDTCO ) EMX(1) = EMX(1) * 1.39d0
           !   ELSE
           !      IF ( NN == IDTCO ) EMX(1) = EMX(1) * 1.19d0
           !   ENDIF
           !ENDIF
           !----------------------------------------------------------------
           IF ( LICARTT ) THEN
              IF ( GET_USA_MASK(I,J) > 0.d0 ) THEN
                 IF ( NN == IDTCO ) EMX(1) = EMX(1) * 1.39d0
              ELSE
                 IF ( NN == IDTCO ) EMX(1) = EMX(1) * 1.19d0
              ENDIF
           ELSE
              IF ( NN == IDTCO ) EMX(1) = EMX(1) * 1.19d0
           ENDIF
        ELSE
           IF ( NN == IDTCO ) EMX(1) = EMX(1) * 1.02d0
        ENDIF

This bug is now fixed in GEOS-Chem v8-02-04.

--Bob Y. 16:40, 23 November 2009 (EST)

Bug fix in DIAG20 (diag_pl_mod.f)

Both Dylan Jones and Yuxuan Wang found a bug in the ND20 diagnostic (production and loss rate of Ox for tagged Ox simulation) in code version 8-02-01 and higher. In routine DIAG20 of diag_pl_mod.f, the COUNT variable for ND20 is not zeroed correctly at the beginning of each new day.

The following lines of code:

        ! Zero counter
        COUNT(I,J,L) = 0    <=== ERROR: (I,J,L) indices used outside loop!

        ! Zero PL24H array
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, N )
        DO N = 1, 2
        DO L = 1, LD65
        DO J = 1, JJPAR
        DO I = 1, IIPAR
           PL24H(I,J,L,N) = 0d0
        ENDDO
        ENDDO
        ENDDO
        ENDDO
!$OMP END PARALLEL DO 

can now be replaced by the F90 array assignment statements:

        ! Zero arrays
        COUNT = 0
        PL24H = 0d0

and then there will be no chance of an out-of-bounds error.

--Bob Y. 10:05, 1 December 2009 (EST)

Div-by-zero error encountered in arsl1k.f

Bob Yantosca (yantosca@seas.harvard.edu) wrote:

While doing some debugging with the SunStudio 12 compiler I came across an error in the chemistry that is a classic div-by-zero error. The Totalview debugger shows that the code is dying at this line in arsl1k.f:
       ! Compute ARSL1K according to the formula listed above
       ARSL1K = AREA / ( RADIUS/DFKG + 2.749064E-4*SQM/(STKCF*STK) )
because the variable STKCF = 0. The value of XSTKCF is computed in calcrate.f for the HO2 uptake in this section:
                ! Test if HO2 het uptake reaction
                ELSE IF ( NK == NKHO2 ) THEN
                   ...
                   XSTKCF = HO2( XRADIUS, T3K(KLOOP), XDENA, XSQM,
   &                             HO2_MOLEC_CM3, N , CONTINENTAL_PBL)
                   ...
                ELSE

                   ! Get GAMMA for species other than N2O5
                   XSTKCF = BRR(NK,NCS)

                ENDIF
and is passed down to arsl1k.f. I noticed in the function HO2 within calcrate.f there are a couple of places where GAMMA can equal zero:
          ! Error check that B^2-(4d0*A*C) is not negative
          TEST= B**2-(4d0*A*C)
          IF ( TEST < 0d0 ) THEN
              GAMMA = 0d0
          ELSE
              ! Calculate the concentration of O2- in the aerosol
              o2_ss= ( -B  -sqrt(B**2-(4d0*A*C)) )/(2d0*A)

              ! Calculate the reactive flux
              fluxrxn = kmt*HO2DENS - o2_ss*Na*kmt/H_eff/Raq/TEMP

              IF ( fluxrxn <= 0d0 ) THEN
                 GAMMA = 0d0
              ELSE
                 ! Gamma for HO2 at TEMP, ho2, and RADIUS given
                 GAMMA = 1./( ( ( HO2DENS/fluxrxn ) -
   &                            ( RADIUS/DFKG ) ) * w / 4.d0 )
              ENDIF
          ENDIF


here the GAMMA is the return value of the function, and this is what gets assigned to XSTCKF in calcrate.f.
It seems the best way to fix this error is to just extend the IF statement within the arsl1k.f function, such that if any of the arguments (which also happen to be in denominators in the equations below) are zero (or at least less than an epsilon value like 1e-30), then we just return the default value of ARSL1K=1d-3. The IF statement within function arsl1k.f should be modified accordingly:
    !----------------------------------------------------------------------
    ! Prior to 12/3/09:
    ! Also check other values to avoid div-by-zero errors (bmy, 12/3/09)
    !IF ( AREA < 0d0 .or. RADIUS < 1d-30 ) THEN
    !----------------------------------------------------------------------
    IF ( AREA < 0d0   .or. DENAIR < 1d-30 .or. RADIUS < 1d-30  .or.
   &     SQM  < 1d-30 .or. STK    < 1d-30 .or. STKCF  < 1d-30 ) THEN

       ! Use default value if any of the above values are zero
       ! This will prevent div-by-zero errors in the eqns below
       ARSL1K = 1.D-3

    ELSE

       ! DFKG = Gas phase diffusion coeff [cm2/s] (order of 0.1)
       DFKG  = 9.45D17/DENAIR * STK * SQRT(3.472D-2 + 1.D0/(SQM*SQM))

       ! Compute ARSL1K according to the formula listed above
       ARSL1K = AREA / ( RADIUS/DFKG + 2.749064E-4*SQM/(STKCF*STK) )

    ENDIF 
The SunStudio compiler checks for NaN, div-by-zero, overflow, underflow, etc. errors when you use the -fast optimization option. The IFORT compiler will not check for these automatically (you can set a compiler switch to invoke this checking).
Long story short: The best defense against bugs like these is to always be careful when writing divisions, and to make sure that the denominator can never be zero (or if it can be zero, to avoid doing the division in the first place).

--Bob Y. 17:01, 3 December 2009 (EST)

Fix for diagnostic arrays in TPCORE

The following previously-reported error:

We found that passing unallocated ND24/25/26 (mass flux) arrays into TPCORE can cause a run to crash. It happened with a modified version of tagged CO simulation in GC v8-01-04 that uses Flambe emissions and was running at 2x2.5 and on multiple processors. Since the mass flux arrays are not optional but intent(inout) for TPCORE, they should be allocated even if they are not used, i.e. even if the diagnostic are turned off.
A general fix that allocates these arrays but lets you turn off the output diagnostic will be in the next GC version (scheduled to be v8-02-01). Meanwhile we strongly advise users to **always** switch on the ND24/25/26 mass flux in their input.geos.
--phs 10:40, 2 April 2009 (EDT)

has now been resolved in GEOS-Chem v8-02-04.

If the ND24, ND25, or ND26 diagnostics are turned off, we simply allocate the MASSFLEW, MASSFLNS, and MASSFLUP arrays using dimensions of 1. This allocates a minimum amount of storage for the arrays, which is enough to prevent the above-described error.

In subroutine ndxx_setup.f, we now have the following IF statement for the ND24 diagnostic:

     !=================================================================
     ! ND24: Eastward mass flux from transport [kg/s] 
     !       --> uses MASSFLEW array (allocatable)
     !=================================================================
     IF ( ND24 > 0 ) THEN
        LD24 = MIN( ND24, LLPAR )
        NMAX = MIN( N_TRACERS, NNPAR )
     
        ALLOCATE( MASSFLEW( IIPAR, JJPAR, LLPAR, NMAX ), STAT=AS) 
        IF ( AS /= 0 ) CALL ALLOC_ERR( 'MASSFLEW' )
     ELSE
        ALLOCATE( MASSFLEW( 1, 1, 1, 1 ), STAT=AS) 
        IF ( AS /= 0 ) CALL ALLOC_ERR( 'MASSFLEW' )
     ENDIF

along with similar IF statements for ND25 and ND26 diagnostics.

--Bob Y. 10:08, 12 January 2010 (EST)

Outstanding issues not yet resolved in v8-02-04