Difference between revisions of "Centralized chemistry time step"

Revision as of 17:59, 20 October 2011

NOTE: Page under construction

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.

Overview

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. However, 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:

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

--Bob Y. 16:55, 7 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 #ITS_A_NEW_YEAR.

Problem when calling planeflight diagnostics

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()-->FLAG = ( MOD( ELAPSED_MIN, TS_DYN ) == 0 )--> Do dynamics
ITS_TIME_FOR_CHEM()-->FLAG = ( MOD( ELAPSED_MIN, TS_CHEM ) == 0 )--> 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() -->FLAG = ( MOD( ELAPSED_MIN, TS_DIAG ) == 0 )-->including call Planeflight.

Since we move chemistry time step to the middle of transport, TTCTT, chemdr.f won't call setup_planeflight anymore. So the model will crash when it calls planeflight during diag output (no initialization).

Fix provided by Bob Y.

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
         print*, '### Called SETUP_PLANEFLIGHT'
         CALL SETUP_PLANEFLIGHT
         DATE_PREV = DATE
      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)

Small fix for chemistry time step

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)