Centralized chemistry time step
On this page we provide information about the time at which the chemistry, emissions, photolysis, and drydep operations are performed.
Contents
Overview
The "chemistry" time step, ΔT(chemistry), must be a multiple of the "transport" time step, ΔT(transport). Typically in GEOS-Chem:
- ΔT(transport) is a function of the grid resolution (30 min for 4° x 5°; 15 min for 2° x 2.5°)
- ΔT(chemistry) has historically been 60 minutes. (This may not always be the optimal setting).
By "transport" we mean the following dynamical operations:
- Advection
- PBL mixing
- Convection
By "chemistry", we mean the solving of the chemical mechanism matrix, plus the related operations that are used to compute reaction rates, i.e.
- Emissions
- Dry deposition
- Photolysis
- Chemistry
Timestep Durations
Optimal Configuration
Sajeev Phillip and Randall Martin recommend the following default GEOS-Chem timestep durations for future GC simulations when sufficient CPU resources are available.
Sajeev Philip wrote:
- Philip et al. (submitted) examines the sensitivity of GEOS-Chem simulations to the duration of chemical and transport operators. Based on that study, we recommend chemical timesteps of 20 minutes and transport timesteps of 10 minutes for simulations when sufficient CPU resources are available. Fine horizontal resolution should take priority over fine temporal resolution. In some cases, it may be beneficial to use coarser timesteps for initial simulations, and the recommended resolution (C20T10) for final simulations. We encourage specification in publications of the duration of operators due to their effects on simulation accuracy.
--Lizzie Lundgren (talk) 16:09, 16 November 2015 (UTC)
Default in v9-01-02 and Later Versions
Currently, in the GEOS-Chem main.F driver program, "transport" is always done before "chemistry". This is due to historical usage. This order may be changed in the future; however we need to make sure that by doing so we do not break the current functionality. We also must keep in mind that the state of the atmosphere (i.e. advected tracer concentrations in the STT array) is used to initialize the "chemistry" operations since (1) gas-phase tracer concentrations are fed into the chemistry solver (SMVGEAR or KPP), and (2) aerosol optical depths are computed from the aerosol tracer concentrations and fed into the Photolysis mechanism. In addition, GEOS-Chem currently computes the photolysis with the sun angles at the center of the "chemistry" timestep. The state of the atmosphere that is used to initialize the "chemistry" should therefore be close to the center of the "chemistry" timestep.
In GEOS-Chem v9-01-02, the time at which "chemistry" is done is closer to the midpoint of the "chemistry" time step, ΔT(chemistry). This update was tested in the 1-month benchmark simulation v9-01-02q and approved on 18 Oct 2011, and is illustrated in the diagram below.
--Bob Y. 16:16, 20 October 2011 (EDT)
Default in Versions Prior to v9-01-02
In GEOS-Chem versions prior to GEOS-Chem v9-01-02, "chemistry" was always done at the beginning of the chemistry timestep (typically at the top of the hour), as illustrated by the following diagram.
--Bob Y. 16:16, 20 October 2011 (EDT)
Source code modifications in v9-01-02
Here we list the source code modifications that were done in order to implement the central chemistry time step algorithm. NOTE: perhaps in the near future a more comprehensive rewrite of the main.F will be necessary in order to fix this more rigorously.
SET_TIMESTEPS
In routine SET_TIMESTEPS (in input_mod.f), we now stop the run if the emissions time step is not the same as the "chemistry" time step, ΔT(chemistry).
! NUNIT is time step in minutes for unit conversion TS_UNIT = -1 ! For the centralizing chemistry timestep modification, the emissions ! should be done at the same time as chemistry. Stop the run if the ! emissions timestep differs from the chemistry timestep. (bmy, 10/7/11) IF ( TS_EMIS /= TS_CHEM ) THEN WRITE(6,*) 'The emission time step should equal the chemistry' WRITE(6,*) 'timestep. This is required for the new central' WRITE(6,*) 'chemistry timestep algorithm.' CALL GEOS_CHEM_STOP ENDIF ! For the centralizing chemistry timestep modification, the convection ! should be done at the same time as dynamics. Stop the run if the ! convection timestep differs from the dynamic timestep. (bmy, 10/7/11) IF ( TS_CONV /= TS_DYN ) THEN WRITE(6,*) 'The convection time step should equal the dynamic' WRITE(6,*) 'timestep. This is required for the new central' WRITE(6,*) 'chemistry timestep algorithm.' CALL GEOS_CHEM_STOP ENDIF
In the same routine, we also no longer need to compute TS_SUN:
!------------------------------------------------------------------------------ ! Prior to 10/7/11: ! We no longer need to pass TS_SUN_2 to SET_TIMESTEPS (bmy, 10/7/11) ! ! We need to compute the cosine of the solar zenith angle at the ! ! center of the relevant timestep. Take 1/2 of the proper time ! ! interval and store for future reference. (bmy, 4/2 7/10) ! IF ( LCHEM .or. LEMIS .or. LDRYD ) THEN ! TS_SUN_2 = TS_CHEM / 2 ! 1/2 of chemistry timestep ! ELSE ! TS_SUN_2 = TS_SMALLEST / 2 ! 1/2 of dynamic (smallest) timestep ! ENDIF !------------------------------------------------------------------------------ ! Initialize timesteps in "time_mod.f" CALL SET_TIMESTEPS( CHEMISTRY = TS_CHEM, EMISSION = TS_EMIS, & DYNAMICS = TS_DYN, UNIT_CONV = TS_UNIT, & CONVECTION = TS_CONV, DIAGNOS = TS_DIAG ) !------------------------------------------------------------------------------ ! Prior to 10/7/11: ! We no longer need to pass TS_SUN_2 to SET_TIMESTEPS (bmy, 10/7/11) ! & SUNCOS = TS_SUN_2 ) !------------------------------------------------------------------------------
ITS_TIME_FOR_CHEM
The current version of function ITS_TIME_FOR_CHEM() in time_mod.f was consistent for computing "chemistry" at the beginning of the chemistry timestep. The following modifications were required in order to compute "chemistry" at the midpoint of the chemistry timestep, ΔT(chemistry).
FLAG = ( MOD( ELAPSED_MIN, TS_CHEM ) == 0 )
to
! changes for proper chemistry time (lzh, ccc, 3/20/10) INTEGER :: M !================================================================= ! ITS_TIME_FOR_CHEM begins here! !================================================================= ! Get multiplier between transport and chemistry: M = TS_CHEM/TS_DYN ! Divide by 2 (get middle). KEEP INTEGERS!!!! M = MAX( M/2, 1 ) ! Is it time for chemistry? !----------------------------------------------------------------------- ! Prior to 9/27/11: ! Chemistry time step is now in the center of transport time steps. The ! previous setup was TCTTT, now TTCTT (lzh, mpayer, 9/27/11) ! FLAG = ( MOD( ELAPSED_MIN, TS_CHEM ) == 0 ) !----------------------------------------------------------------------- ! Chemistry time step is now in the center of transport time steps FLAG = ( MOD( ELAPSED_MIN, TS_CHEM ) == (M-1)*TS_DYN )
ITS_TIME_FOR_EMIS
Because emissions are done on the "chemistry" time step, we must make the same modifications in ITS_TIME_FOR_EMIS as we did in ITS_TIME_FOR_CHEM.
! changes for proper chemistry time (lzh, ccc, 3/20/10) INTEGER :: M !================================================================= ! ITS_TIME_FOR_EMIS begins here! !================================================================= ! Get multiplier between transport and chemistry: M = TS_EMIS/TS_DYN ! Divide by 2 (get middle). KEEP INTEGERS!!!! M = MAX( M/2, 1 ) ! Is it time for emissions? !----------------------------------------------------------------------- ! Prior to 9/27/11: ! Emission time step is now in the center of transport time steps. The ! previous setup was TCTTT, now TTCTT (lzh, mpayer, 9/27/11) ! FLAG = ( MOD( ELAPSED_MIN, TS_EMIS ) == 0 ) !----------------------------------------------------------------------- ! Emission time step is now in the center of transport time steps FLAG = ( MOD( ELAPSED_MIN, TS_EMIS ) == (M-1)*TS_DYN )
ITS_A_NEW_YEAR
The routine ITS_A_NEW_YEAR (in time_mod.f) returns .TRUE. if it is the first timestep of a new year where we have to read data from disk. We needed to rewrite this routine so that it returns .TRUE. not at 00:00 GMT on Jan 1st, but at the midpoint of the first chemistry timestep on Jan 1st.
NOTE: Some code is omitted for clarity.
INTEGER :: HH, MM, SS REAL*8 :: GMTb, TS ! Emissions timestep [hours] TS = DBLE( TS_EMIS ) / 60d0 ... !============================================================== ! FOR JANUARY 1st OF THE YEAR !============================================================== IF ( MONTH == 1 .and. DAY == 1 ) THEN ... ! Here, we are using the central chemistry timestep option. ! Therefore, the first data read of the year will occur not at ! 00:00 GMT on the Jan 1st, but offset by a small amount (as ! diagnosed by function ITS_TIME_FOR_EMIS). IS_NEW_YEAR = ( GMT < TS .and. ITS_TIME_FOR_EMIS() ) ELSE IF ( NYMD == NYMDb ) THEN !============================================================== ! FOR THE FIRST DAY OF THE SIMULATION ! (i.e. for simulations that start at other times of the year) !============================================================== ... ! Split starting time into hour, minute, second CALL YMD_EXTRACT( NHMSb, HH, MM, SS ) ! Compute GMT at the start of the simulation GMTb = DBLE( HH ) + ( DBLE( MM ) / 60d0 ) ! Here, we are using the central chemistry timestep option. ! Therefore, the first data read of the year will occur not ! at 00:00 GMT on Jan 1st, but offset by a small amount (as ! diagnosed by function ITS_TIME_FOR_EMIS). IS_NEW_YEAR = ( GMT < GMTb+TS .and. ITS_TIME_FOR_EMIS() ) ... ELSE !============================================================== ! FOR ALL OTHER DAYS !============================================================== ! It isn't time for the first data read of the year; return FALSE IS_NEW_YEAR = .FALSE. ENDIF
ITS_A_NEW_MONTH
Source code modifications are similar to those in function ITS_A_NEW_YEAR.
ITS_A_NEW_DAY
Source code modifications are similar to those in function ITS_A_NEW_YEAR.
ITS_MIDMONTH
This function is used to determine when it is time to read a new month of leaf area index data. The only change we made is to now call function ITS_A_NEW_DAY instead of testing if it is 00:00 GMT on the 16th of the month.
! Test for the 16th of the month at 0 GMT !---------------------------------------------------------------------------- ! Prior to 10/14/11: ! Now use ITS_A_NEW_DAY to tell if it's the first data read time of the day. ! This will account for the central chemistry timestep option. (bmy, 10/14/11) ! IS_MIDMONTH = ( DAY == 16 .and. NHMS == 000000 ) !---------------------------------------------------------------------------- IS_MIDMONTH = ( DAY == 16 .and. ITS_A_NEW_DAY() )
COSSZA
Routine COSSZA in dao_mod.F is used to compute the cosine of the solar zenith angle, which is required for photolysis. We have rewritten this routine to return two quantities:
- cosine( solar zenith angle ) @ the current time.
- This can be used to test if it is the end of a day.
- cosine( solar zenith angle ) @ the midpoint of the chemistry timestep.
- This is passed to the chemistry, emissions, drydep, and photolysis routines.
The subroutine now looks like this (some code omitted for clarity).:
SUBROUTINE COSSZA( JDAY, SUNCOS, SUNCOS_MID ) ... USE GRID_MOD, ONLY : GET_YMID_R USE TIME_MOD, ONLY : GET_HOUR USE TIME_MOD, ONLY : GET_MINUTE USE TIME_MOD, ONLY : GET_LOCALTIME USE TIME_MOD, ONLY : GET_TS_CHEM ... INTEGER, INTENT(IN) :: JDAY ! Day of year (0-365) ... REAL*8, INTENT(OUT) :: SUNCOS(MAXIJ) ! cos(SZA) @ current time REAL*8, INTENT(OUT) :: SUNCOS_MID(MAXIJ) ! cos(SZA) @ midpt of chem step ... INTEGER :: FACTOR, I, IJLOOP, J INTEGER :: HOUR, MINUTE, TS_SUN REAL*8 :: AHR, GMT_MID, R, TIMLOC REAL*8 :: DEC, C_DEC, S_DEC REAL*8 :: YMID_R, C_YMID_R, S_YMID_R ! Coefficients for solar declination angle REAL*8, PARAMETER :: A0 = 0.006918d0 REAL*8, PARAMETER :: A1 = 0.399912d0 REAL*8, PARAMETER :: A2 = 0.006758d0 REAL*8, PARAMETER :: A3 = 0.002697d0 REAL*8, PARAMETER :: B1 = 0.070257d0 REAL*8, PARAMETER :: B2 = 0.000907d0 REAL*8, PARAMETER :: B3 = 0.000148d0 !================================================================= ! Initialization !================================================================= ! Quantities for central chemistry timestep TS_SUN = GET_TS_CHEM() ! Chemistry interval HOUR = GET_HOUR() ! Current hour MINUTE = GET_MINUTE() ! Current minutes FACTOR = MINUTE / TS_SUN ! Multiplying factor ! GMT at the midpoint of the chemistry timestep GMT_MID = ( DBLE( GET_HOUR() ) ) ! Start of hour & + ( DBLE( TS_SUN * FACTOR ) / 60d0 ) ! + steps already done & + ( DBLE( TS_SUN / 2 ) / 60d0 ) ! + 1/2 chem timestep ! Path length of earth's orbit traversed since Jan 1 [radians] R = ( 2d0 * PI / 365d0 ) * DBLE( JDAY - 1 ) ! Solar declination angle (low precision formula) [radians] DEC = A0 - A1*COS( R ) + B1*SIN( R ) & - A2*COS( 2d0*R ) + B2*SIN( 2d0*R ) & - A3*COS( 3d0*R ) + B3*SIN( 3d0*R ) ! Pre-compute sin & cos of DEC outside of DO loops (for efficiency) S_DEC = SIN( DEC ) C_DEC = COS( DEC ) !================================================================= ! Compute cosine of solar zenith angle !================================================================= !$OMP PARALLEL DO !$OMP+DEFAULT( SHARED ) !$OMP+PRIVATE( I, J, YMID_R, S_YMID_R, C_YMID_R, IJLOOP, TIMLOC, AHR ) ! Loop over latitude DO J = 1, JJPAR ! Latitude of grid box [radians] YMID_R = GET_YMID_R( J ) ! Pre-compute sin & cos of DEC outside of I loop (for efficiency) S_YMID_R = SIN( YMID_R ) C_YMID_R = COS( YMID_R ) ! Loop over longitude DO I = 1, IIPAR ! 1-D grid box index IJLOOP = ( (J-1) * IIPAR ) + I !=========================================================== ! Compute cosine of SZA at the current GMT time !=========================================================== ! Local time at box (I,J) [hours] TIMLOC = GET_LOCALTIME( I ) ! Hour angle at box (I,J) [radians] AHR = ABS( TIMLOC - 12d0 ) * 15d0 * PI_180 ! Cosine of solar zenith angle at box (I,J) [unitless] SUNCOS(IJLOOP) = ( S_YMID_R * S_DEC ) & + ( C_YMID_R * C_DEC * COS( AHR ) ) !=========================================================== ! Compute cosine(SZA) at the midpoint of the chem timestep ! Required for photolysis, chemistry, emissions, drydep !=========================================================== ! Local time [hours] at box (I,J) at 30 minutes after the GMT hour TIMLOC = GET_LOCALTIME( I, GMT=GMT_MID ) ! Hour angle at box (I,J) [radians] AHR = ABS( TIMLOC - 12d0 ) * 15d0 * PI_180 ! Corresponding cosine( SZA ) at box (I,J) [unitless] SUNCOS_MID(IJLOOP) = ( S_YMID_R * S_DEC ) & + ( C_YMID_R * C_DEC * COS( AHR ) ) ENDDO ENDDO !$OMP END PARALLEL DO END SUBROUTINE COSSZA
chemdr.F
In routine chemdr.F, we must pass the SUNCOS_MID array to the FAST_J and PHYSPROC routines. This will make sure that photolysis and chemistry are evaluated at the midpoint of the chemistry timestep.
USE DAO_MOD, ONLY : SUNCOS_MID ... CALL FAST_J( SUNCOS_MID, OPTD, UVALBEDO ) ... CALL PHYSPROC( SUNCOS_MID )
DO_DRYDEP
We now pass SUNCOS_MID to the dry deposition routine DEPVEL (in drydep_mod.f)
USE DAO_MOD, ONLY : SUNCOS_MID ... CALL DEPVEL( MAXIJ, RADIAT, TC0, SUNCOS_MID, F0, HSTAR, & XMW, AIROSOL, USTAR, CZ1, OBK, CFRAC, & ZH, LSNOW, DVEL, AZO, RHB )
emissdr.F
We now have to pass SUNCOS_MID to the soil NOx emissions routines:
& CALL SOILNOXEMS( SUNCOS_MID )
As well as to the biogenic emissions routines:
SC = SUNCOS_MID(IJLOOP) ... EMIS = EMISOP( I,J,IJLOOP,SUNCOS_MID,TMMP,XNUMOL_C ) ... EMMB = EMISOP_MB( I, J, IJLOOP, SUNCOS_MID, ... GRASS = EMISOP_GRASS(I,J,IJLOOP,SUNCOS_MID,TMMP,XNUMOL_C)
carbon_mod.F
We have make similar changes in carbon_mod.f as in emissdr.F
USE DAO_MOD, ONLY : SUNCOS_MID ... SC = SUNCOS_MID(IJLOOP) ... SC = SUNCOS_MID(IJLOOP) ... IF ( SUNCOS_MID(IJLOOP) > 0d0 .and. TCOSZ(I,J) > 0d0 ) THEN
main.F
We must pass the SUNCOS_MID array to routine COSSZA in main.F:
USE DAO_MOD, ONLY: SUNCOS_MID ... ! SUNCOS_MID = cos(SZA) at the midpt of the chemistry timestep CALL COSSZA( DAY_OF_YEAR, SUNCOS, SUNCOS_MID )
We have
USE DAO_MOD, ONLY : SUNCOS_MID SC = SUNCOS_MID(IJLOOP) SC = SUNCOS_MID(IJLOOP) IF ( SUNCOS_MID(IJLOOP) > 0d0 .and. TCOSZ(I,J) > 0d0 ) THEN
--Bob Y. 16:34, 20 October 2011 (EDT)
diag_pl_mod.F
Lin Zhang provided the following fix. In GeosCore/diag_pl_mod.f, change line 1021
IF ( MOD( HOUR + CHEM, 24d0 ) == 0 ) THEN
to
IF ( HOUR + CHEM .GE. 24d0 ) THEN
--Melissa Payer 14:36, 27 September 2011 (EDT)
Modifications for planeflight diagnostics
The planeflight following diagnostic (ND40) had to be modified slightly for the central chemistry timestep algorithm:
Description of problem
This code caused the model to crash when you call planeflight diagnostics.
Because originally in main.f, TS_DYN = 15min TS_CHEM = 60min TS_DIAG = 60min TCTTT (T is "transport", C is "chemistry") ITS_TIME_FOR_DYN() --> Do dynamics ITS_TIME_FOR_CHEM() --> Do chemistry (call setup_planeflight based on ITS_A_NEW_DAY in chemdr.f) CALL TIMESTAMP_DIAG CALL SET_ELAPSED_MIN --> ELAPSED_MIN = ELAPSED_MIN + TS_DYN (after this is done, it is not a new day any more.) CALL SET_CURRENT_TIME ITS_TIME_FOR_DIAG() --> call the Planeflight diagnostics
Because we have now centralized the "chemistry" in between a symmetric number of "transport" operations (TTCTT), routine chemdr.f won't call routine SETUP_PLANEFLIGHT anymore. Therefore, the model will crash when it calls planeflight during diag output (no initialization).
Solution (by Bob Y.)
To fix this situation, replace these lines in chemdr.f
IF ( ND40 > 0 .and. ITS_A_NEW_DAY() ) THEN CALL SETUP_PLANEFLIGHT ENDIF
with:
INTEGER :: DATE INTEGER, SAVE :: DATE_PREV = -1 ... !================================================================= ! At the beginning of each new day, call SETUP_PLANEFLIGHT ! to see if there are any plane flight points to be processed !================================================================= ! Get todays' date DATE = GET_NYMD() ! If this is the first chem timestep of a new day, then we need to ! call SETUP_PLANEFLIGHT. If chemistry is turned on, then we need ! to place this call here, so as to make sure that the chemical ! mechanism files (read by READER and READCHEM) have been loaded. IF ( ND40 .and. DATE /= DATE_PREV ) THEN CALL SETUP_PLANEFLIGHT DATE_PREV = DATE ENDIF
NOTE: This probably could also be expressed as:
IF ( ND40 .and. ITS_A_NEW_DAY() ) THEN CALL SETUP_PLANEFLIGHT ENDIF
The idea is that you want to call SETUP_PLANEFLIGHT on the first chemical timestep of a new day. However, this is no longer just when NHMS == 000000. A better way to test that just to see when the date is not the same as the previous date.
--Jmao 14:17, 14 June 2011 (EDT)
--Bob Y. 16:16, 20 October 2011 (EDT)
References
- Santillana, M., L. Zhang, D.J. Jacob, R. Yantosca, Estimating numerical errors due to operator splitting in global atmospheric chemistry models: Transport and chemistry, in preparation, 2011.
--Bob Y. 11:05, 21 October 2011 (EDT)