Centralized chemistry time step

On this page we provide information about the modification made to GEOS-Chem v9-01-02 which moves the time at which the chemistry, emissions, photolysis, and drydep operations are performed.

This update was tested in the 1-month benchmark simulation v9-01-02q and approved on 18 Oct 2011.

Overview

Terminology

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

The "chemistry" time step, ΔT(chemistry), is done 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).

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 must keep in mind that the state of the atmosphere (i.e. advected tracer concentrations in the STT array) are used to initialize the "chemistry" operations as follows:

• Gas-phase tracer concentrations are fed into the chemistry solver (SMVGEAR or KPP)
• Aerosol optical depths are computed from the aerosol tracer concentrations and fed into the Photolysis mechanism.

GEOS-Chem currently computes the photolysis with the sun angles at the center of the "chemistry" timestep. This being the case, we should also make sure that the state of the atmosphere that is used to initialize the "chemistry" should be close to the center of the "chemistry" timestep.

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 (i.e. typically at the top of the hour). The following diagram illustrates this:

v9-01-02 and higher versions

In GEOS-Chem v9-01-02, we have implemented a fix to move the time at which "chemistry" is done closer to the midpoint of the "chemistry" time step, ΔT(chemistry).

The source code modifications that were necessary are described below. NOTE: perhaps in the near future a more comprehensive rewrite of the main.F will be necessary in order to fix this more rigorously.

--Bob Y. 15:52, 20 October 2011 (EDT)

Source code modifications

Here we list the source code modifications that were done in order to implement the central chemistry time step algorithm.

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:

1. cosine( solar zenith angle ) @ the current time.
   • This can be used to test if it is the end of a day.
2. 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 )

--Bob Y. 16:34, 20 October 2011 (EDT) 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

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

1. 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)