Difference between revisions of "GEOS-Chem v8-03-01"
(→Overview) |
(→Fix to prevent numerical overflow in EXP function) |
||
Line 296: | Line 296: | ||
--[[User:Bmy|Bob Y.]] 16:05, 11 January 2010 (EST) | --[[User:Bmy|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 | ||
== Outstanding issues not yet resolved in v8-03-01 == | == Outstanding issues not yet resolved in v8-03-01 == |
Revision as of 14:45, 19 April 2010
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.
Contents
Overview
PUBLIC RELEASE -- March/April 2010
Will contain everything in GEOS-Chem v8-02-04, plus:
- Option to use Caltech isoprene chemistry (F. Paulot)
- ISORROPIA II (H. Pye, T. Nenes)
- Updated aerosol optical properties for FAST-J photolysis (R. Martin, C. Heald)
- Modifications to SOA formation (Aerosols Working Group)
- SOA formation from aromatics (D. Henze)
- Speciated biogenic emissions from MEGAN v2.1 now used in SOA code (H. O. T. Pye)
- Global 1 x 1.25 simulation capability (L. Lamsal)
- Modifications to / Creation of the GEOS-5 0.5 x 0.667 N. American and European Nested Grid simulation (A. van Donkelaar)
- TOMAS aerosol microphysics option (W. Trivitayanurak, D. Westervelt, J. Pierce)
- Extension of annual anthropogenic scale factors to 2006 (A. van Donkelaar)
- EMEP emissions extended to 2007 and Seasonality extended to SOx, CO, and NH3 (A. van Donkelaar)
Previous issues now resolved in v8-03-01
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:
ELSE !============================================================== ! 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 ENDIF ENDIF
--Bob Y. 16:21, 15 March 2010 (EDT)
Bug in GET_EMMONOT_MEGAN
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)
Bug in READ_ANTHRO_NH3
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:
ALLOCATE( EMEP_SO2_SHIP( IIPAR, JJPAR ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' ) EMEP_SO2_SHIP = 0d0 ALLOCATE( EMEP_CO_SHIP( IIPAR, JJPAR ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' ) EMEP_SO2_SHIP = 0d0 ALLOCATE( EMEP_NOx_SHIP( IIPAR, JJPAR ), STAT=AS ) 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:
ALLOCATE( EMEP_SO2_SHIP( IIPAR, JJPAR ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_SO2_SHIP' ) EMEP_SO2_SHIP = 0d0 ALLOCATE( EMEP_CO_SHIP( IIPAR, JJPAR ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'EMEP_CO_SHIP' ) EMEP_CO_SHIP = 0d0 ALLOCATE( EMEP_NOx_SHIP( IIPAR, JJPAR ), STAT=AS ) 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 RETURN #endif #endif
- 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 ) etc.
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.
!$OMP PARALLEL DO DEFAULT( SHARED ) PRIVATE( I, J ) 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_ALPH ) EMIS_SAVE(I,J,IDTALPH) = BIOG_ALPH(I,J) 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) ENDIF 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) ENDIF ENDDO ENDDO !$OMP END PARALLEL DO
- 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:
USE ERROR_MOD, ONLY : IS_SAFE_EXP . . . 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 ENDIF !---------------------------------------------------------- ! 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) ELSE ! If exponential can't be computed, set to ! the initial O3 concentration L3S = O3 ENDIF SO2_ss = MAX( SO2_ss - ( L2S + L3S ), MINDAT ) H2O20 = MAX( H2O20 - L2S, MINDAT ) etc...
--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