Centralized chemistry time step: Difference between revisions
| Line 385: | Line 385: | ||
USE DAO_MOD, ONLY : SUNCOS_MID | USE DAO_MOD, ONLY : SUNCOS_MID | ||
... | |||
CALL DEPVEL( MAXIJ, RADIAT, TC0, SUNCOS_MID, F0, HSTAR, | CALL DEPVEL( MAXIJ, RADIAT, TC0, SUNCOS_MID, F0, HSTAR, | ||
& XMW, AIROSOL, USTAR, CZ1, OBK, CFRAC, | & XMW, AIROSOL, USTAR, CZ1, OBK, CFRAC, | ||
& ZH, LSNOW, DVEL, AZO, RHB ) | & ZH, LSNOW, DVEL, AZO, RHB ) | ||
=== emissdr.F === | === emissdr.F === | ||
Revision as of 20:29, 20 October 2011
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
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:
- 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<tt 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)
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)