Difference between revisions of "Centralized chemistry time step"
(→Overview) |
(→Overview) |
||
Line 18: | Line 18: | ||
* Chemistry | * Chemistry | ||
− | The "chemistry" time step ΔT(chemistry) is done must be a multiple of the "transport" time step ΔT(transport). Typically in GEOS-Chem: | + | 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 | + | * Δ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). | * ΔT(chemistry) has historically been 60 minutes. (This may not always be the optimal setting). | ||
− | Currently, in the GEOS-Chem <tt>main.F</tt> driver program, "transport" is 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. | + | Currently, in the GEOS-Chem <tt>main.F</tt> 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 <tt>STT</tt> array) are used to initialize the "chemistry" operations as follows: | |
− | * Gas-phase tracer concentrations are fed into the | + | * 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]]. | * 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. | + | 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: |
− | + | ||
− | However, in GEOS-Chem versions prior to [[GEOS-Chem v9-01-02]], | + | |
[[Image:Tcttt_configuration.jpg]] | [[Image:Tcttt_configuration.jpg]] | ||
− | + | 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). | |
[[Image:Ttctt_configuration.jpg]] | [[Image:Ttctt_configuration.jpg]] | ||
− | + | The source code modifications that were necessary are described below. NOTE: perhaps in the near future a more comprehensive rewrite of the <tt>main.F</tt> will be necessary in order to | |
− | + | --[[User:Bmy|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 <tt>SET_TIMESTEPS</tt> (in <tt>input_mod.f</tt>), 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 <tt>TS_SUN</tt>: | |
− | In | + | |
− | < | + | |
− | + | ||
− | </ | + | |
− | to | + | !------------------------------------------------------------------------------ |
− | + | ! 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 <tt>ITS_TIME_FOR_CHEM()</tt> in <tt>time_mod.f</tt> 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 === | ||
In this way, chemistry time step is put in the center of transport time steps. | In this way, chemistry time step is put in the center of transport time steps. | ||
Revision as of 17:47, 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.
Contents
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
In this way, chemistry time step is put in the center of transport time steps.
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)