402 lines
16 KiB
Fortran
402 lines
16 KiB
Fortran
! $Id: schem.f,v 1.1 2009/06/09 21:51:52 daven Exp $
|
|
SUBROUTINE SCHEM
|
|
!
|
|
!******************************************************************************
|
|
! Subroutine SCHEM performs simplified stratospheric chemistry, which means
|
|
! only reactions with OH and photolysis are considered. The production and
|
|
! loss of CO and NOy in the stratosphere are taken from Dylan Jones' 2-D
|
|
! model. (qli, bmy, 11/20/1999, 10/25/05)
|
|
!
|
|
! NOTES:
|
|
! (1 ) Now read all inputs (stratospheric OH, monthly mean J-values,
|
|
! P(CO) rates, and L(CO) rates) from binary punch file format.
|
|
! (bmy, 12/10/99)
|
|
! (2 ) Uses READ_BPCH2 to read from binary file format (bmy, 12/10/99)
|
|
! (3 ) Make sure the DO-loops go in the order N-L-J-I to avoid disk
|
|
! swapping problems (bmy, 12/10/99)
|
|
! (4 ) Remove reactions for HNO3 photolysis and HNO3 + OH. The HNO3
|
|
! concentrations that we read in from disk are from Dylan's 2-D
|
|
! model, where chemistry is already taken into account.
|
|
! (qli, bmy, 12/23/99)
|
|
! (5 ) Remove obsolete code from 12/23/99. (bmy, 4/18/00)
|
|
! (6 ) Bug fixes: Cap RDLOSS so that it does not exceed 1.0.
|
|
! Now declare RDLOSS, T1L, RC, K0, K1, K2, K3, M as REAL*8
|
|
! Cosmetic changes & update comments (bmy, 5/4/00)
|
|
! (7 ) Reference F90 module "bpch2_mod" which contains routine "read_bpch2"
|
|
! for reading data from binary punch files (bmy, 6/28/00)
|
|
! (8 ) Now all monthly mean J-values are in the same file (bmy, 6/30/00)
|
|
! (9 ) Now use function GET_TAU0 (from "bpch2_mod.f") to return the TAU0
|
|
! value used to index the binary punch file. (bmy, 7/20/00)
|
|
! (10) Declared arrays for reading data from disk to be both ALLOCATABLE
|
|
! and SAVE. Also cosmetic changes & some cleanup. (bmy, 9/8/00)
|
|
! (11) Activated parallel DO-loops (bmy, 12/12/00)
|
|
! (12) Now use 3 arguments (M/D/Y) in call to GET_TAU0. ARRAY needs to be
|
|
! of size (IGLOB,JGLOB). Use JGLOB,LGLOB in calls to READ_BPCH2.
|
|
! Use TRANSFER_ZONAL (from "transfer_mod.f") to cast from REAL*4 to
|
|
! REAL*8 and resize arrays to (JJPAR,LLPAR). Updated comments,
|
|
! made cosmetic changes. (bmy, 9/27/01)
|
|
! (13) Removed obsolete commented out code from 9/01 (bmy, 10/24/01)
|
|
! (14) Now read COprod and COloss files directly from the
|
|
! DATA_DIR/pco_lco_200203/ subdirectory. Also read stratOH files
|
|
! directly from the DATA_DIR/stratOH_200203/ subdirectory. Also
|
|
! read stratjv files directly from the DATA_DIR/stratjv_200203/
|
|
! subdirectory. (bmy, 4/2/02)
|
|
! (15) Now reference AD and T from "dao_mod.f". Also reference routine
|
|
! ALLOC_ERR from "error_mod.f". Now reference IDTOX, IDTNOX, etc.
|
|
! from "tracerid_mod.f". (bmy, 11/6/02)
|
|
! (16) Now use functions GET_TS_CHEM, GET_MONTH and GET_TAU, and
|
|
! TIMESTAMP_STRING from the new "time_mod.f". Also call READ_BPCH2
|
|
! with QUIET=.TRUE., which prevents info from being printed to the
|
|
! log file. (bmy, 3/14/03)
|
|
! (17) LINUX has a problem putting a function call w/in a WRITE statement.
|
|
! Now save output from TIMESTAMP_STRING to STAMP and print that.
|
|
! (bmy, 9/29/03)
|
|
! (18) Now reference STT and TRACER_MW_KG from "tracer_mod.f". Now reference
|
|
! DATA_DIR from "directory_mod.f". Bug fix: now loop over N_TRACERS
|
|
! and not NNPAR. NNPAR is the max # of tracers but may not be the
|
|
! actual number of tracers. (bmy, 7/20/04)
|
|
! (19) Now references GET_MIN_TPAUSE_LEVEL and ITS_IN_THE_STRAT from
|
|
! "tropopause_mod.f". Now remove reference to CMN, it's obsolete.
|
|
! (bmy, 8/22/05)
|
|
! (20) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
|
|
! (21) Now references XNUMOLAIR from "tracer_mod.f" (bmy, 10/25/05)
|
|
!******************************************************************************
|
|
!
|
|
! References to F90 modules
|
|
USE BPCH2_MOD, ONLY : GET_NAME_EXT, GET_RES_EXT
|
|
USE BPCH2_MOD, ONLY : GET_TAU0, READ_BPCH2
|
|
USE DAO_MOD, ONLY : AD, T
|
|
USE DIRECTORY_MOD, ONLY : DATA_DIR
|
|
USE ERROR_MOD, ONLY : ALLOC_ERR
|
|
USE TIME_MOD, ONLY : GET_MONTH, GET_TAU
|
|
USE TIME_MOD, ONLY : GET_TS_CHEM, TIMESTAMP_STRING
|
|
USE TRACER_MOD, ONLY : N_TRACERS, STT
|
|
USE TRACER_MOD, ONLY : TRACER_MW_KG, XNUMOLAIR
|
|
USE TRACERID_MOD, ONLY : IDTACET, IDTALD2, IDTALK4, IDTC2H6
|
|
USE TRACERID_MOD, ONLY : IDTC3H8, IDTCH2O, IDTH2O2, IDTHNO4
|
|
USE TRACERID_MOD, ONLY : IDTISOP, IDTMACR, IDTMEK, IDTMP
|
|
USE TRACERID_MOD, ONLY : IDTMVK, IDTPMN, IDTPRPE, IDTR4N2
|
|
USE TRACERID_MOD, ONLY : IDTRCHO
|
|
USE TRANSFER_MOD, ONLY : TRANSFER_ZONAL
|
|
USE TROPOPAUSE_MOD, ONLY : GET_MIN_TPAUSE_LEVEL, ITS_IN_THE_STRAT
|
|
|
|
IMPLICIT NONE
|
|
|
|
# include "CMN_SIZE" ! Size parameters
|
|
|
|
! Local variables
|
|
LOGICAL, SAVE :: FIRST = .TRUE.
|
|
|
|
INTEGER :: I, IOS, J, L, N, NN, LMIN
|
|
INTEGER, SAVE :: MONTHSAVE = 0
|
|
|
|
! Number of photolysis species (currently is 13)
|
|
INTEGER, PARAMETER :: NSPHOTO = 13
|
|
|
|
! Tracers that undergo photolysis loss in the stratosphere
|
|
INTEGER :: SPHOTOID(NSPHOTO) = (/
|
|
& 3, 8, 9, 10, 11, 12, 13,
|
|
& 14, 17, 20, 22, 23, 24 /)
|
|
|
|
! Character variables
|
|
CHARACTER(LEN=16 ) :: STAMP
|
|
CHARACTER(LEN=255) :: FILENAME
|
|
|
|
! REAL*4 arrays -- for reading from binary data files
|
|
REAL*4 :: ARRAY(1,JGLOB,LGLOB)
|
|
REAL*4, ALLOCATABLE, SAVE :: STRATOH(:,:)
|
|
REAL*4, ALLOCATABLE, SAVE :: SJVALUE(:,:,:)
|
|
REAL*4, ALLOCATABLE, SAVE :: COPROD(:,:)
|
|
REAL*4, ALLOCATABLE, SAVE :: COLOSS(:,:)
|
|
|
|
! REAL*8 variables
|
|
REAL*8 :: k0, k1, k2, k3, XTAU
|
|
REAL*8 :: DTCHEM, RDLOSS, T1L, M, TK, RC
|
|
|
|
! External functions
|
|
REAL*8, EXTERNAL :: BOXVL
|
|
|
|
!=================================================================
|
|
! SCHEM begins here!
|
|
!=================================================================
|
|
|
|
! Chemistry timestep [s]
|
|
DTCHEM = GET_TS_CHEM() * 60d0
|
|
|
|
! Echo info
|
|
STAMP = TIMESTAMP_STRING()
|
|
WRITE( 6, 100 ) STAMP
|
|
100 FORMAT( ' - SCHEM: Strat chemistry at ', a )
|
|
|
|
!=================================================================
|
|
! If it is the first call to SCHEM, allocate arrays for reading
|
|
! data. These arrays are declared SAVE so they will be preserved
|
|
! between calls.
|
|
!=================================================================
|
|
IF ( FIRST ) THEN
|
|
ALLOCATE( STRATOH( JJPAR, LLPAR ), STAT=IOS )
|
|
IF ( IOS /= 0 ) CALL ALLOC_ERR( 'STRATOH' )
|
|
STRATOH = 0e0
|
|
|
|
ALLOCATE( SJVALUE( JJPAR, LLPAR, NSPHOTO ), STAT=IOS )
|
|
IF ( IOS /= 0 ) CALL ALLOC_ERR( 'SJVALUE' )
|
|
SJVALUE = 0e0
|
|
|
|
ALLOCATE( COPROD( JJPAR, LLPAR ), STAT=IOS )
|
|
IF ( IOS /= 0 ) CALL ALLOC_ERR( 'COPROD' )
|
|
COPROD = 0e0
|
|
|
|
ALLOCATE( COLOSS( JJPAR, LLPAR ), STAT=IOS )
|
|
IF ( IOS /= 0 ) CALL ALLOC_ERR( 'COLOSS' )
|
|
COLOSS = 0e0
|
|
ENDIF
|
|
|
|
!=================================================================
|
|
! If it is a new month (or the first call to SCHEM),
|
|
! do the following:
|
|
!
|
|
! (1) Read archived J-values and store in SJVALUE
|
|
! (2) Read archived CO production rates and store in COPROD
|
|
! (3) Read archived CO loss rates and store in COLOSS
|
|
!
|
|
! NOTES
|
|
! (a) All of the above-mentioned data are stored in binary punch
|
|
! files, for ease of use.
|
|
!
|
|
! (b) STRATOH, SJVALUE, CO_PROD, and CO_LOSS are now declared
|
|
! as both ALLOCATABLE and SAVE. If SCHEM is called, then
|
|
! data will be declared for these arrays, and the values in
|
|
! these arrays will be preserved between calls.
|
|
!
|
|
! (c) If SCHEM is never called (i.e. if you are running another
|
|
! type of chemistry simulation), then memory never gets
|
|
! allocated to STRATOH, SJVALUE, CO_PROD, and CO_LOSS.
|
|
! This saves on computational resources.
|
|
!=================================================================
|
|
IF ( GET_MONTH() /= MONTHSAVE .or. FIRST ) THEN
|
|
MONTHSAVE = GET_MONTH()
|
|
|
|
! TAU value at the beginning of this month
|
|
XTAU = GET_TAU0( GET_MONTH(), 1, 1985 )
|
|
|
|
!==============================================================
|
|
! Read this month's OH
|
|
!==============================================================
|
|
FILENAME = TRIM( DATA_DIR ) // 'stratOH_200203/stratOH.' //
|
|
& GET_NAME_EXT() // '.' //
|
|
& GET_RES_EXT()
|
|
|
|
! Read data
|
|
CALL READ_BPCH2( FILENAME, 'CHEM-L=$', 1,
|
|
& XTAU, 1, JGLOB,
|
|
& LGLOB, ARRAY, QUIET=.TRUE. )
|
|
|
|
! Cast from REAL*4 to REAL*8 and resize to (JJPAR,LLPAR)
|
|
CALL TRANSFER_ZONAL( ARRAY(1,:,:), STRATOH )
|
|
|
|
!==============================================================
|
|
! Read in monthly mean archived J-values
|
|
!==============================================================
|
|
FILENAME = TRIM( DATA_DIR ) // 'stratjv_200203/stratjv.' //
|
|
& GET_NAME_EXT() // '.' //
|
|
& GET_RES_EXT()
|
|
|
|
DO NN = 1, NSPHOTO
|
|
N = SPHOTOID(NN)
|
|
|
|
! Read data
|
|
CALL READ_BPCH2( FILENAME, 'JV-MAP-$', N,
|
|
& XTAU, 1, JGLOB,
|
|
& LGLOB, ARRAY, QUIET=.TRUE. )
|
|
|
|
! Cast from REAL*4 to REAL*8 and resize to (JJPAR,LLPAR)
|
|
CALL TRANSFER_ZONAL( ARRAY(1,:,:), SJVALUE(:,:,NN) )
|
|
ENDDO
|
|
|
|
!==============================================================
|
|
! Read in CO production rates
|
|
!==============================================================
|
|
FILENAME = TRIM( DATA_DIR ) // 'pco_lco_200203/COprod.' //
|
|
& GET_NAME_EXT() // '.' //
|
|
& GET_RES_EXT()
|
|
|
|
! Read data
|
|
CALL READ_BPCH2( FILENAME, 'PORL-L=$', 9,
|
|
& XTAU, 1, JGLOB,
|
|
& LGLOB, ARRAY, QUIET=.TRUE. )
|
|
|
|
! Cast from REAL*4 to REAL*8 and resize to (JJPAR,LLPAR)
|
|
CALL TRANSFER_ZONAL( ARRAY(1,:,:), COPROD )
|
|
|
|
!==============================================================
|
|
! Read in CO loss rates
|
|
!==============================================================
|
|
FILENAME = TRIM( DATA_DIR ) // 'pco_lco_200203/COloss.' //
|
|
& GET_NAME_EXT() // '.' //
|
|
& GET_RES_EXT()
|
|
|
|
! Read data
|
|
CALL READ_BPCH2( FILENAME, 'PORL-L=$', 10,
|
|
& XTAU, 1, JGLOB,
|
|
& LGLOB, ARRAY, QUIET=.TRUE. )
|
|
|
|
! Cast from REAL*4 to REAL*8 and resize to (JJPAR,LLPAR)
|
|
CALL TRANSFER_ZONAL( ARRAY(1,:,:), COLOSS )
|
|
|
|
ENDIF
|
|
|
|
!=================================================================
|
|
! Do photolysis for selected tracers with this
|
|
! month's archived J-values
|
|
!=================================================================
|
|
|
|
! Get the minimum level extent of the ann mean tropopause
|
|
LMIN = GET_MIN_TPAUSE_LEVEL()
|
|
|
|
!$OMP PARALLEL DO
|
|
!$OMP+DEFAULT( SHARED )
|
|
!$OMP+PRIVATE( I, J, L, N, NN )
|
|
!$OMP+SCHEDULE( DYNAMIC )
|
|
DO NN = 1, NSPHOTO
|
|
N = SPHOTOID(NN)
|
|
|
|
DO L = LMIN, LLPAR
|
|
DO J = 1, JJPAR
|
|
DO I = 1, IIPAR
|
|
|
|
! Only proceed for stratospheric boxes
|
|
IF ( ITS_IN_THE_STRAT( I, J, L ) ) THEN
|
|
|
|
! Compute photolysis loss
|
|
STT(I,J,L,N) = STT(I,J,L,N) *
|
|
& EXP( -SJVALUE(J,L,NN) * DTCHEM )
|
|
ENDIF
|
|
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
!print*, 'In schem, done with photolysis'
|
|
|
|
!=================================================================
|
|
! CO is special --
|
|
! use archived P, L rates for CO chemistry in stratosphere
|
|
!=================================================================
|
|
CALL CO_STRAT_PL( COPROD, COLOSS )
|
|
|
|
!=================================================================
|
|
! Reaction with OH -- compute rate constants for each tracer
|
|
!=================================================================
|
|
!print*, 'In schem, before reaction with OH'
|
|
|
|
!$OMP PARALLEL DO
|
|
!$OMP+DEFAULT( SHARED )
|
|
!$OMP+PRIVATE( I, J, L, N, M, TK, RC, k0, k1, RDLOSS, T1L )
|
|
!$OMP+SCHEDULE( DYNAMIC )
|
|
DO N = 1, N_TRACERS
|
|
DO L = LMIN, LLPAR
|
|
DO J = 1, JJPAR
|
|
DO I = 1, IIPAR
|
|
|
|
! Only proceed for stratospheric boxes
|
|
IF ( ITS_IN_THE_STRAT( I, J, L ) ) THEN
|
|
|
|
! Density of air at grid box (I,J,L) in molec/cm3
|
|
M = AD(I,J,L) / BOXVL(I,J,L) * XNUMOLAIR
|
|
|
|
! Temperature at grid box (I,J,L) in K
|
|
TK = T(I,J,L)
|
|
|
|
! Select proper reaction rate w/ OH for the given tracer
|
|
! Some rates are temperature or density dependent
|
|
IF ( N == IDTALK4 ) THEN
|
|
RC = 8.20D-12 * EXP( -300.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTISOP ) THEN
|
|
RC = 2.55D-11 * EXP( 410.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTH2O2 ) THEN
|
|
RC = 2.90D-12 * EXP( -160.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTACET ) THEN
|
|
RC = 1.70D-12 * EXP( -600.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTMEK ) THEN
|
|
RC = 2.92D-13 * EXP( 414.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTALD2 ) THEN
|
|
RC = 1.40D-12 * EXP( -1860.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTRCHO ) THEN
|
|
RC = 2.00D-11
|
|
|
|
ELSE IF ( N == IDTMVK ) THEN
|
|
RC = 4.13D-12 * EXP( 452.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTMACR ) THEN
|
|
RC = 1.86D-11 * EXP( -175.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTPMN ) THEN
|
|
RC = 3.60D-12
|
|
|
|
ELSE IF ( N == IDTR4N2 ) THEN
|
|
RC = 1.30D-12
|
|
|
|
ELSE IF ( N == IDTPRPE ) THEN
|
|
k0 = 8.0D-27 * ( 300.D0 / TK )**3.5
|
|
k1 = 3.0D-11
|
|
|
|
RC = k1 * k0 * M / ( k1 + k0*M )
|
|
RC = RC * 0.5 ** (1 / ( 1 + LOG10( k0*M/k1 )**2 ) )
|
|
|
|
ELSE IF ( N == IDTC3H8 ) THEN
|
|
RC = 8.00D-12 * EXP( -590.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTCH2O ) THEN
|
|
RC = 1.00D-12
|
|
|
|
ELSE IF ( N == IDTC2H6 ) THEN
|
|
RC = 7.9D-12 * EXP( -1030.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTHNO4 ) THEN
|
|
RC = 1.30D-12 * EXP( 380.D0 / TK )
|
|
|
|
ELSE IF ( N == IDTMP ) THEN
|
|
RC = 1.14D-12 * EXP( 200.D0 / TK )
|
|
|
|
ELSE
|
|
RC = 0d0
|
|
|
|
ENDIF
|
|
|
|
! Compute loss with OH based on the rate constants from above
|
|
! Cap RDLOSS so that it does not exceed 1.0 (bmy, 5/4/00)
|
|
RDLOSS = RC * STRATOH(J,L) * DTCHEM
|
|
RDLOSS = MIN( RDLOSS, 1d0 )
|
|
|
|
! T1L is the absolute amount of STT lost to rxn with OH
|
|
! Subtract T1L from STT
|
|
T1L = STT(I,J,L,N) * RDLOSS
|
|
STT(I,J,L,N) = STT(I,J,L,N) - T1L
|
|
|
|
! Oxidation of PRPE as source of ACET with 80% yield
|
|
IF ( N == IDTPRPE ) THEN
|
|
STT(I,J,L,IDTACET) = STT(I,J,L,IDTACET) +
|
|
& 0.8d0 * T1L *
|
|
& TRACER_MW_KG(IDTACET) / TRACER_MW_KG(IDTPRPE)
|
|
ENDIF
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
! Set FIRST = .FALSE. -- we have been thru SCHEM at least once now
|
|
FIRST = .FALSE.
|
|
|
|
! Return to calling program
|
|
END SUBROUTINE SCHEM
|