! $Id: sulfate_mod.f,v 1.36 2007/03/29 20:31:23 bmy Exp $
MODULE SULFATE_MOD 3
!
!******************************************************************************
! Module SULFATE_MOD contains arrays and routines for performing either a
! coupled chemistry/aerosol run or an offline sulfate aerosol simulation.
! Original code taken from Mian Chin's GOCART model and modified accordingly.
! (rjp, bdf, bmy, 6/22/00, 12/8/06)
!
! Module Variables:
! ============================================================================
! (1 ) XNUMOL_OH (REAL*8 ) : Molecules OH per kg OH [molec/kg]
! (2 ) XNUMOL_O3 (REAL*8 ) : Molecules O3 per kg O3 [molec/kg]
! (3 ) XNUMOL_NO3 (REAL*8 ) : Molecules NO3 per kg NO3 [molec/kg]
! (4 ) TCVV_S (REAL*8 ) : Ratio: Molwt air / Molwt S [unitless]
! (5 ) DMSo (REAL*8 ) : DMS oceanic emissions [v/v/timestep]
! (6 ) DRYH2O2 (INTEGER) : Pointer to H2O2 in DEPVEL array [unitless]
! (7 ) DRYSO2 (INTEGER) : Pointer to SO2 in DEPVEL array [unitless]
! (8 ) DRYSO4 (INTEGER) : Pointer to SO4 in DEPVEL array [unitless]
! (9 ) DRYSO4s (INTEGER) : Pointer to SO4s in DEPVEL array [unitless]
! (10) DRYMSA (INTEGER) : Pointer to MSA in DEPVEL array [unitless]
! (11) DRYNH3 (INTEGER) : Pointer to NH3 in DEPVEL array [unitless]
! (12) DRYNH4 (INTEGER) : Pointer to NH4 in DEPVEL array [unitless]
! (13) DRYNIT (INTEGER) : Pointer to NIT in DEPVEL array [unitless]
! (14) DRYNITs (INTEGER) : Pointer to NITs in DEPVEL array [unitless]
! (15) DRYSO4aq (INTEGER) : Pointer to SO4aq in DEPVEL array [unitless]
! (16) DRYAS (INTEGER) : Pointer to AS in DEPVEL array [unitless]
! (17) DRYAHS (INTEGER) : Pointer to AHS in DEPVEL array [unitless]
! (18) DRYLET (INTEGER) : Pointer to LET in DEPVEL array [unitless]
! (19) DRYNH4aq (INTEGER) : Pointer to NH4aq in DEPVEL array [unitless]
! (20) ENH3_an (REAL*8 ) : NH3 anthropogenic emissions [kg NH3/box/s]
! (21) ENH3_bb (REAL*8 ) : NH3 biomass emissions [kg NH3/box/s]
! (22) ENH3_bf (REAL*8 ) : NH3 biofuel emissions [kg NH3/box/s]
! (23) ENH3_na (REAL*8 ) : NH73 natural source emissions [kg NH3/box/s]
! (24) ESO2_ac (REAL*8 ) : SO2 aircraft emissions [kg SO2/box/s]
! (25) ESO2_an (REAL*8 ) : SO2 anthropogenic emissions [kg SO2/box/s]
! (26) ESO2_ev (REAL*8 ) : SO2 eruptive volcanic em. [kg SO2/box/s]
! (27) ESO2_nv (REAL*8 ) : SO2 non-eruptive volcanic em. [kg SO2/box/s]
! (28) ESO2_bb (REAL*8 ) : SO2 biomass burning emissions [kg SO2/box/s]
! (29) ESO2_bf (REAL*8 ) : SO2 biofuel burning emissions [kg SO2/box/s]
! (30) ESO2_sh (REAL*8 ) : SO2 ship emissions [kg SO2/box/s]
! (31) ESO4_an (REAL*8 ) : SO4 anthropogenic emissions [kg SO4/box/s]
! (32) JH2O2 (REAL*8 ) : Monthly mean J(H2O2) values [s-1]
! (33) O3m (REAL*8 ) : Monthly mean O3 concentration [v/v]
! (34) PH2O2m (REAL*8 ) : Monthly mean P(H2O2) [molec/cm3/s]
! (35) PMSA_DMS (REAL*8 ) : P(MSA) from DMS [v/v/timestep]
! (36) PSO2_DMS (REAL*8 ) : P(SO2) from DMS [v/v/timestep]
! (37) PSO4_SO2 (REAL*8 ) : P(SO4) from SO2 [v/v/timestep]
! (38) SSTEMP (REAL*8 ) : Sea surface temperatures [K]
! (39) VCLDF (REAL*8 ) : Volume cloud frac. for SO2 aq. [unitless]
! (40) NEV (INTEGER) : Max # of eruptive volcanoes [unitless]
! (41) IEV (INTEGER) : Longitudes of eruptive volcanoes [degrees]
! (42) JEV (INTEGER) : Latitudes of eruptive volcanoes [degrees ]
! (43) IHGHT (INTEGER) : Height of eruptive volcano plume [m]
! (44) IELVe (INTEGER) : Elevation of eruptive volcanoes [m]
! (45) Eev (REAL*8 ) : SO2 em. from eruptive volcanoes [kg SO2/box/s]
! (46) NNV (INTEGER) : Max # of non-eruptive volcanoes [unitless]
! (47) NNVOL (INTEGER) : Number of non-eruptive volcanoes [unitless]
! (48) INV (INTEGER) : Longitude of non-erup volcanoes [degrees]
! (49) JNV (INTEGER) : Latitude of non-erup volcanoes [degrees]
! (50) IELVn (INTEGER) : Elevation of non-erup volcanoes [m]
! (51) Env (INTEGER) : SO2 em. from non-erup volcanoes [kg SO2/box/s]
! (52) TCOSZ (REAL*8 ) : Sum of cos(SZA) for offline run [unitless]
! (53) TTDAY (REAL*8 ) : Total daylight length at (I,J) [minutes]
! (54) SMALLNUM (REAL*8 ) : Small number - prevent underflow [unitless]
! (55) COSZM (REAL*8 ) : Array for MAX(cos(SZA)) at (I,J) [unitless]
!
! Module Routines:
! ===========================================================================
! (1 ) GET_VCLDF : Computes volume cloud fraction for SO2 chemistry
! (2 ) GET_LWC : Computes liquid water content as a function of T
! (3 ) CHEMSULFATE : Driver routine for sulfate/aerosol chemistry
! (4 ) GRAV_SETTLING : Routine to compute settling of SO4s and NITs
! (5 ) CHEM_DMS : Chemistry routine for DMS tracer
! (6 ) CHEM_H2O2 : Chemistry routine for H2O2 tracer
! (7 ) CHEM_SO2 : Chemistry routine for SO2 tracer
! (8 ) SEASALT_CHEM : Computes SO2->SO4 and HNO3->nitrate w/in seasalt
! (9 ) AQCHEM_SO2 : Computes reaction rates for aqueous SO2 chemistry
! (10) CHEM_SO4 : Chemistry routine for SO4 tracer
! (11) PHASE_SO4 : Computes phase transition for crystalline tracers
! (12) PHASE_RADIATIVE : Computes radiative forcing for crystalline tracers
! (13) CHEM_MSA : Chemistry routine for MSA tracer
! (14) CHEM_NH3 : Chemistry routine for ammonia tracer
! (15) CHEM_NH4 : Chemistry routine for ammonium tracer
! (16) CHEM_NIT : Chemistry routine for nitrates tracer
! (17) EMISSSULFATE : Driver routine for sulfate/aerosol emissions
! (18) SRCDMS : Emission routine for DMS tracer
! (19) SRCSO2 : Emission routine for SO2 tracer
! (20) SRCSO4 : Emission routine for SO4 tracer
! (21) SRCNH3 : Emission routine for NH3 tracer
! (22) GET_OH : Returns OH for coupled or offline simulations
! (23) SET_OH : Resets modified OH in SMVGEAR's CSPEC array
! (24) GET_NO3 : Returns NO3 for coupled or offline simulations
! (25) SET_NO3 : Resets modified OH in SMVGEAR's CSPEC array
! (26) GET_O3 : Returns O3 for coupled or offline simulations
! (27) READ_NONERUP_VOLC : Reads SO2 emissions from non-eruptive volcanoes
! (28) READ_ERUP_VOLC : Reads SO2 emissions from eruptive volcanoes
! (29) READ_ANTHRO_SOx : Reads anthropogenic SO2 and SO4 emissions
! (30) READ_OCEAN_DMS : Reads biogenic DMS emissions from oceans
! (31) READ_SST : Reads monthly mean sea-surface temperatures
! (32) READ_BIOMASS_SO2 : Reads SO2 emissions from biomass burning
! (33) READ_AIRCRAFT_SO2 : Reads SO2 emissions from aircraft exhaust
! (34) READ_SHIP_SO2 : Reads SO2 emissions from ship exhaust
! (35) READ_ANTHRO_NH3 : Reads NH3 emissions from anthropogenic sources
! (36) READ_NATURAL_NH3 : Reads NH3 emissions from natural sources
! (37) READ_BIOMASS_NH3 : Reads NH3 biomass burning emissions
! (38) READ_OXIDANT : Reads monthly mean O3 and H2O2 for offline run
! (39) OHNO3TIME : Computes time arrays for scaling offline OH, NO3
! (40) INIT_SULFATE : Allocates & zeroes module arrays
! (41) CLEANUP_SULFATE : Deallocates module arrays
!
! GEOS-CHEM modules referenced by sulfate_mod.f
! ============================================================================
! (1 ) biomass_mod.f : Module w/ routines for biomass burning
! (2 ) bpch2_mod.f : Module w/ routines for binary pch file I/O
! (3 ) bravo_mod.f : Module w/ routines to read BRAVO emissions
! (4 ) comode_mod.f : Module w/ SMVGEAR allocatable arrays
! (5 ) dao_mod.f : Module w/ DAO met field arrays
! (6 ) diag_mod.f : Module w/ GEOS-Chem diagnostic arrays
! (7 ) directory_mod.f : Module w/ GEOS-Chem data & met field dirs
! (8 ) drydep_mod.f : Module w/ GEOS-Chem dry deposition routines
! (9 ) epa_nei_mod.f : Module w/ routines to read EPA/NEI99 data
! (10) error_mod.f : Module w/ NaN, other error check routines
! (11) file_mod.f : Module w/ file unit numbers & error checks
! (12) future_emissions_mod.f : Module w/ routines for IPCC future emissions
! (13) grid_mod.f : Module w/ horizontal grid information
! (14) global_no3_mod.f : Module w/ routines to read 3-D NO3 field
! (15) global_oh_mod.f : Module w/ routines to read 3-D OH field
! (16) isoropia_mod.f : Module w/ ISORROPIA routines for aer thermo eq
! (17) logical_mod.f : Module w/ GEOS-Chem logical switches
! (18) pbl_mix_mod.f : Module w/ routines for PBL height & mixing
! (19) pressure_mod.f : Module w/ routines to compute P(I,J,L)
! (20) seasalt_mod.f : Module w/ routines for seasalt chemistry
! (21) streets_anthro_mod.f : Module w/ routines for David Streets' emiss
! (22) time_mod.f : Module w/ routines to compute time & date
! (23) tracer_mod.f : Module w/ GEOS-Chem tracer array STT etc.
! (24) tracerid_mod.f : Module w/ pointers to tracers & emissions
! (25) transfer_mod.f : Module w/ routines to cast & resize arrays
! (26) tropopause_mod.f : Module w/ routines to read ann mean tropopause
! (27) uvalbedo_mod.f : Module w/ UV albedo array and reader
! (28) wetscav_mod.f : Module w/ routines for wetdep & scavenging
!
! References:
! ============================================================================
! (1 ) Andreae, M.O. & P. Merlet, "Emission of trace gases and aerosols from
! biomass burning", Global Biogeochem. Cycles, 15, 955-966, 2001.
! (2 ) Nightingale et al [2000a], J. Geophys. Res, 14, 373-387
! (3 ) Nightingale et al [2000b], Geophys. Res. Lett, 27, 2117-2120
! (4 ) Wanninkhof, R., "Relation between wind speed and gas exchange over
! the ocean", J. Geophys. Res, 97, 7373-7382, 1992.
!
! NOTES:
! (1 ) All module variables are declared PRIVATE (i.e., they can only
! be seen from within this module (bmy, 6/2/00)
! (2 ) The routines in "sulfate_mod.f" assume that we are doing chemistry
! over the global region (e.g. IIPAR=IGLOB, JJPAR=JGLOB). (bmy, 6/8/00)
! (3 ) Removed obsolete code from DRYDEP_SULFATE (bmy, 12/21/00)
! (4 ) Removed obsolete commented-out code from module routines (bmy, 4/23/01)
! (5 ) Now read data files from DATA_DIR/sulfate_sim_200106/ (bmy, 6/19/01)
! (6 ) Updated comments (bmy, 9/4/01)
! (7 ) XTRA2(IREF,JREF,5) is now XTRA2(I,J). Now reference COSSZA from
! "dao_mod.f". (bmy, 9/27/01)
! (8 ) Removed obsolete commented out code from 9/01 (bmy, 10/24/01)
! (9 ) Minor fixes to facilitate compilation on ALPHA (bmy, 11/15/01)
! (11) Now divide module header into MODULE PRIVATE, MODULE VARIABLES, and
! MODULE ROUTINES sections. Updated comments (bmy, 5/28/02)
! (12) Replaced all instances of IM with IIPAR and JM with JJPAR, in order
! to prevent namespace confusion for the new TPCORE (bmy, 6/25/02)
! (13) Now reference "file_mod.f" (bmy, 6/27/02)
! (14) Now references GET_PEDGE from "pressure_mod.f", which computes P at
! the bottom edge of grid box (I,J,L). Also deleted obsolete,
! commented-out code. (dsa, bdf, bmy, 8/21/02)
! (15) Added updated code from Rokjin Park and Brendan Field, in order to
! perform coupled chemistry-aerosol simulations. Also added parallel
! DO-loops in several subroutines. Updated comments, cosmetic
! changes. Now reference "error_mod.f" and "wetscav_mod.f".
! Now only do chemistry below the tropopause. (rjp, bdf, bmy, 12/6/02)
! (16) Added ENH3_na array to hold natural source NH3 emissions. Also now
! facilitate passing DMS, SO2, SO4, NH3 to SMVGEAR for fullchem
! simulations. Added subroutine READ_NATURAL_NH3. (rjp, bmy, 3/23/03)
! (17) Now references "grid_mod.f" and "time_mod.f". Also made other minor
! cosmetic changes. (bmy, 3/27/03)
! (18) Updated chemistry routines to apply drydep losses throughout the
! entire PBL. (rjp, bmy, 8/1/03)
! (19) Now accounts for GEOS-4 PBL being in meters (bmy, 1/15/04)
! (20) Fix ND44 diag so that we get same results for sp or mp (bmy, 3/24/04)
! (21) Added COSZM array. Now use diurnal varying JH2O2 in CHEM_H2O2.
! (rjp, bmy, 3/39/04)
! (22) Added more parallel DO-loops (bmy, 4/14/04)
! (23) Now add SO2 from ships (bec, bmy, 5/20/04)
! (24) Now references "directory_mod.f", "logical_mod.f" and "tracer_mod.f".
! Now removed IJSURF. (bmy, 7/20/04)
! (25) Can overwrite USA with EPA/NEI99 emissions (rjp, rch, bmy, 11/16/04)
! (26) Modified for AS, AHS, LET, SO4aq, NH4aq (cas, bmy, 1/11/05)
! (27) Now also references "pbl_mix_mod.f". NOTE: Comment out phase
! transition code for now since it is still under development and
! will take a while to be rewritten. (bmy, 3/15/05)
! (28) Modified for SO4s, NITs chemistry (bec, 4/13/05)
! (29) Now reads updated files for SST and offline chemistry. Now read data
! for both GCAP and GEOS grids. Now references "tropopause_mod.f".
! (bmy, 8/22/05)
! (30) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
! (31) Now references XNUMOL & XNUMOLAIR from "tracer_mod.f" (bmy, 10/25/05)
! (32) Now read int'annual SST data on GEOS 1x1 grid (bmy, 11/17/05)
! (33) Bug fix for offline aerosol sim in SEASALT_CHEM (bec, bmy, 3/29/06)
! (34) Bug fix in INIT_DRYDEP (bmy, 5/23/06)
! (35) Now references "bravo_mod.f" (rjp, kfb, bmy, 6/26/06)
! (36) Now references "streets_anthro_mod.f" (yxw, bmy, 8/17/06)
! (37) Now references "biomass_mod.f" (bmy, 9/27/06)
! (38) Now prevent seg fault error in READ_BIOMASS_SO2 (bmy, 11/3/06)
! (39) Bug fix in SEASALT_CHEM (havala, bec, bmy, 12/8/06)
!******************************************************************************
!
IMPLICIT NONE
!=================================================================
! MODULE PRIVATE DECLARATIONS -- keep certain internal variables
! and routines from being seen outside "sulfate_mod.f"
!=================================================================
! Make everything PRIVATE ...
PRIVATE
! ... except these routines
PUBLIC :: CHEMSULFATE
PUBLIC :: EMISSSULFATE
PUBLIC :: CLEANUP_SULFATE
!=================================================================
! MODULE VARIABLES (see descriptions listed above)
!=================================================================
! Time variable
INTEGER :: ELAPSED_SEC
! Logical Flags
LOGICAL, PARAMETER :: LENV = .TRUE.
LOGICAL, PARAMETER :: LEEV = .TRUE.
! Parameters
REAL*8, PARAMETER :: XNUMOL_OH = 6.022d23 / 17d-3
REAL*8, PARAMETER :: XNUMOL_O3 = 6.022d23 / 48d-3
REAL*8, PARAMETER :: XNUMOL_NO3 = 6.022d23 / 62d-3
REAL*8, PARAMETER :: TCVV_S = 28.97d0 / 32d0
REAL*8, PARAMETER :: SMALLNUM = 1d-20
! Allocatable arrays
REAL*8, ALLOCATABLE :: DMSo(:,:)
REAL*8, ALLOCATABLE :: ENH3_an(:,:)
REAL*8, ALLOCATABLE :: ENH3_bb(:,:)
REAL*8, ALLOCATABLE :: ENH3_bf(:,:)
REAL*8, ALLOCATABLE :: ENH3_na(:,:)
REAL*8, ALLOCATABLE :: ESO2_ac(:,:,:)
REAL*8, ALLOCATABLE :: ESO2_an(:,:,:)
REAL*8, ALLOCATABLE :: ESO2_bb(:,:)
REAL*8, ALLOCATABLE :: ESO2_bf(:,:)
REAL*8, ALLOCATABLE :: ESO2_ev(:,:,:)
REAL*8, ALLOCATABLE :: ESO2_nv(:,:,:)
REAL*8, ALLOCATABLE :: ESO2_sh(:,:)
REAL*8, ALLOCATABLE :: ESO4_an(:,:,:)
REAL*8, ALLOCATABLE :: JH2O2(:,:,:)
REAL*8, ALLOCATABLE :: LSO2_AQ(:,:,:)
REAL*8, ALLOCATABLE :: O3m(:,:,:)
REAL*8, ALLOCATABLE :: PH2O2m(:,:,:)
REAL*8, ALLOCATABLE :: PMSA_DMS(:,:,:)
REAL*8, ALLOCATABLE :: PSO2_DMS(:,:,:)
REAL*8, ALLOCATABLE :: PSO4_SO2(:,:,:)
REAL*8, ALLOCATABLE :: PSO4_SS(:,:,:)
REAL*8, ALLOCATABLE :: PNITs(:,:,:)
REAL*8, ALLOCATABLE :: SOx_SCALE(:,:)
REAL*8, ALLOCATABLE :: SSTEMP(:,:)
REAL*8, ALLOCATABLE :: TCOSZ(:,:)
REAL*8, ALLOCATABLE :: TTDAY(:,:)
REAL*8, ALLOCATABLE :: VCLDF(:,:,:)
REAL*8, ALLOCATABLE :: COSZM(:,:)
! Eruptive volcanoes
INTEGER, PARAMETER :: NEV=50
INTEGER :: NEVOL
INTEGER, ALLOCATABLE :: IEV(:), JEV(:)
INTEGER, ALLOCATABLE :: IDAYs(:), IDAYe(:)
INTEGER, ALLOCATABLE :: IHGHT(:), IELVe(:)
REAL*8, ALLOCATABLE :: EEV(:)
! Non-eruptive volcanoes
INTEGER, PARAMETER :: NNV=50
INTEGER :: NNVOL
INTEGER, ALLOCATABLE :: INV(:), JNV(:), IELVn(:)
REAL*8, ALLOCATABLE :: ENV(:)
! Pointers to drydep species w/in DEPSAV
INTEGER :: DRYSO2, DRYSO4, DRYMSA, DRYNH3
INTEGER :: DRYNH4, DRYNIT, DRYSO4s, DRYNITs
INTEGER :: DRYH2O2, DRYSO4aq, DRYAS, DRYAHS
INTEGER :: DRYLET, DRYNH4aq
!=================================================================
! MODULE ROUTINES -- follow below the "CONTAINS" statement
!=================================================================
CONTAINS
!------------------------------------------------------------------------------
SUBROUTINE GET_VCLDF 2,8
!
!******************************************************************************
! Subroutine GET_VCLDF computes the volume cloud fraction for SO2 chemistry.
! (rjp, bdf, bmy, 9/23/02)
!
! References:
! ============================================================================
! (1) Sundqvist et al. [1989]
!
! NOTES:
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD
, ONLY : RH
USE PRESSURE_MOD, ONLY : GET_PCENTER, GET_PEDGE
# include "CMN_SIZE" ! Size parameters
! Local variables
INTEGER :: I, J, L
REAL*8 :: PRES, PSFC, RH2, R0, B0
! Parameters
REAL*8, PARAMETER :: ZRT = 0.60d0, ZRS = 0.99d0
!=================================================================
! GET_VCLDF begins here!
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, PSFC, PRES, RH2, R0, B0 )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
! Surface pressure
PSFC = GET_PEDGE(I,J,1)
! Pressure at the center of the grid box
PRES = GET_PCENTER(I,J,L)
! RH (from "dao_mod.f") is relative humidity [%]
! Convert to fraction and store in RH2
RH2 = RH(I,J,L) * 1.0d-2
! Terms from Sundqvist ???
R0 = ZRT + ( ZRS - ZRT ) * EXP( 1d0 - ( PSFC / PRES )**2.5 )
B0 = ( RH2 - R0 ) / ( 1d0 - R0 )
! Force B0 into the range 0-1
IF ( RH2 < R0 ) B0 = 0d0
IF ( B0 > 1d0 ) B0 = 1d0
! Volume cloud fraction
VCLDF(I,J,L) = 1d0 - SQRT( 1d0 - B0 )
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE GET_VCLDF
!------------------------------------------------------------------------------
FUNCTION GET_LWC( T ) RESULT( LWC ) 2
!
!******************************************************************************
! Function GET_LWC returns the cloud liquid water content at a GEOS-CHEM
! grid box as a function of temperature. (rjp, bmy, 10/31/02, 1/14/03)
!
! Arguments as Input:
! ============================================================================
! (1 ) T (REAL*8) : Temperature value at a GEOS-CHEM grid box [K]
!
! NOTES:
!******************************************************************************
!
! Arguments
REAL*8, INTENT(IN) :: T
! Function value
REAL*8 :: LWC
!=================================================================
! GET_LWC begins here!
!=================================================================
! Compute Liquid water content in [g/m3]
IF ( T > 293d0 ) THEN
LWC = 0.2d0
ELSE IF ( T >= 280.d0 .AND. T <= 293.d0 ) THEN
LWC = 0.32d0 - 0.0060d0 * ( T - 273.D0 )
ELSE IF ( T >= 248.d0 .AND. T < 280.d0 ) THEN
LWC = 0.23d0 + 0.0065d0 * ( T - 273.D0 )
ELSE IF ( T < 248.d0 ) THEN
LWC = 0.07d0
ENDIF
! Convert from [g/m3] to [m3/m3]
LWC = LWC * 1.D-6
! Return to calling program
END FUNCTION GET_LWC
!------------------------------------------------------------------------------
SUBROUTINE CHEMSULFATE 2,47
!
!******************************************************************************
! Subroutine CHEMSULFATE is the interface between the GEOS-CHEM main program
! and the sulfate chemistry routines. The user has the option of running
! a coupled chemistry-aerosols simulation or an offline aerosol simulation.
! (rjp, bdf, bmy, 5/31/00, 3/16/06)
!
! NOTES:
! (1 ) Now reference all arguments except FIRSTCHEM and RH from either F90
! modules or from common block header files. Updated comments,
! cosmetic changes. Added NH3, NH4, NITRATE chemistry routines.
! Also call MAKE_RH and CONVERT_UNITS from "dao_mod.f". Now references
! IDTDMS, IDTSO2 etc. from "tracerid_mod.f". Now make FIRSTCHEM a
! local SAVEd variable. Now reference DEPSAV from "drydep_mod.f".
! Also get rid of extraneous dimensions of DEPSAV. Added NTIME,
! NHMSb arrays for OHNO3TIME. (rjp, bdf, bmy, 12/16/02)
! (2 ) CHEM_DMS is now only called for offline sulfate simulations.
! (rjp, bmy, 3/23/03)
! (3 ) Now remove NTIME, NHMSb from the arg list and call to OHNO3TIME.
! Now references functions GET_MONTH, GET_TS_CHEM, and GET_ELAPSED_SEC
! from the new "time_mod.f". (bmy, 3/27/03)
! (4 ) Now reference STT, TCVV, N_TRACERS, ITS_AN_AEROSOL_SIM from
! "tracer_mod.f". Now reference ITS_A_NEW_MONTH from "time_mod.f".
! Now references LPRT from "logical_mod.f". (bmy, 7/20/04)
! (5 ) Updated for AS, AHS, LET, SO4aq, NH4aq. Now references LCRYST from
! logical_mod.f. Now locate species in the DEPSAV array w/in
! INIT_SULFATE. (bmy, 12/21/04)
! (6 ) Now handle gravitational settling of SO4s, NITs (bec, bmy, 4/13/05)
! (7 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
! (8 ) Remove reference to MAKE_RH, it's not needed here (bmy, 3/16/06)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD, AIRDEN, CLDF
USE DAO_MOD, ONLY : SUNCOS, CONVERT_UNITS
USE DRYDEP_MOD, ONLY : DEPSAV
USE ERROR_MOD, ONLY : DEBUG_MSG
USE GLOBAL_OH_MOD, ONLY : GET_GLOBAL_OH
USE GLOBAL_NO3_MOD, ONLY : GET_GLOBAL_NO3
USE LOGICAL_MOD, ONLY : LCRYST, LPRT
USE TIME_MOD, ONLY : GET_MONTH, GET_TS_CHEM
USE TIME_MOD, ONLY : GET_ELAPSED_SEC, ITS_A_NEW_MONTH
USE TRACER_MOD, ONLY : STT, TCVV
USE TRACER_MOD, ONLY : N_TRACERS, ITS_AN_AEROSOL_SIM
USE TRACERID_MOD, ONLY : IDTNITs, IDTSO4s
# include "CMN_SIZE" ! Size parameters
! Local variables
LOGICAL, SAVE :: FIRSTCHEM = .TRUE.
INTEGER, SAVE :: LASTMONTH = -99
INTEGER :: I, J, L, N, MONTH
REAL*8 :: DTCHEM
! External functions
REAL*8, EXTERNAL :: BOXVL
!=================================================================
! CHEMSULFATE begins here!
!=================================================================
! Get current month
MONTH = GET_MONTH()
! Establish indices w/in DEPSAV array
IF ( FIRSTCHEM ) THEN
! Initialize arrays (if not already done before)
CALL INIT_SULFATE
! Reset first-time flag
FIRSTCHEM = .FALSE.
ENDIF
! If it's an offline simulation ...
IF ( ITS_AN_AEROSOL_SIM() ) THEN
! Then read monthly data files ...
IF ( ITS_A_NEW_MONTH() ) THEN
CALL GET_GLOBAL_OH( MONTH )
CALL GET_GLOBAL_NO3( MONTH )
ENDIF
! And compute time scaling arrays for offline OH, NO3
CALL OHNO3TIME
ENDIF
! Store NTIME in a shadow variable
ELAPSED_SEC = GET_ELAPSED_SEC()
! DTCHEM is the chemistry timestep in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Initialize module arrays
PSO2_DMS = 0d0
PMSA_DMS = 0d0
PSO4_SO2 = 0d0
PSO4_SS = 0d0
PNITs = 0d0
!=================================================================
! Call individual chemistry routines for sulfate/aerosol tracers
!=================================================================
! SO4s [kg] gravitational settling
CALL GRAV_SETTLING( STT(:,:,:,IDTSO4s), 1 )
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: GRAV_SET, SO4S' )
! NITs [kg] gravitational settling
CALL GRAV_SETTLING( STT(:,:,:,IDTNITs), 2 )
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: GRAV_SET, NITS' )
! Convert all tracers in STT from [kg] -> [v/v]
CALL CONVERT_UNITS( 1, N_TRACERS, TCVV, AD, STT )
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CONVERT UNITS' )
! For offline runs only ...
IF ( ITS_AN_AEROSOL_SIM() ) THEN
! DMS (offline only)
CALL CHEM_DMS
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_DMS' )
! H2O2 (offline only)
CALL CHEM_H2O2
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_H2O2' )
ENDIF
! SO2
CALL GET_VCLDF
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a get VCLDF' )
CALL CHEM_SO2
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_SO2' )
! SO4
CALL CHEM_SO4
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_SO4' )
!-----------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% ! Only do the following if the crystalline sulfate & aqueous
!%%% ! tracers (AS, AHS, LET, SO4aq, NH4aq) are defined
!%%% IF ( LCRYST ) THEN
!%%%
!%%% ! Phase change
!%%% CALL PHASE_SO4
!%%% IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a PHASE_SO4' )
!%%%
!%%% ! Radiative forcing
!%%% CALL PHASE_RADIATIVE
!%%% IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a PHASE_RAD' )
!%%%
!%%% ENDIF
!-----------------------------------------------------------------------------
! MSA
CALL CHEM_MSA
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_MSA' )
! NH3
CALL CHEM_NH3
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_NH3' )
! NH4 (gas-phase)
CALL CHEM_NH4
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_NH4' )
!-----------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% ! NH4 (aqueous phase)
!%%% IF ( LCRYST ) THEN
!%%% CALL CHEM_NH4aq
!%%% IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_NH4aq' )
!%%% ENDIF
!-----------------------------------------------------------------------------
! Sulfur Nitrate
CALL CHEM_NIT
IF ( LPRT ) CALL DEBUG_MSG( '### CHEMSULFATE: a CHEM_NIT' )
! Convert STT from [v/v] -> [kg]
CALL CONVERT_UNITS( 2, N_TRACERS, TCVV, AD, STT )
! We have already gone thru one chemistry iteration
FIRSTCHEM = .FALSE.
! Return to calling program
END SUBROUTINE CHEMSULFATE
!------------------------------------------------------------------------------
SUBROUTINE GRAV_SETTLING( TC, N ) 2,12
!
!******************************************************************************
! Subroutine GRAV_SETTLING performs gravitational settling of sulfate
! and nitrate in coarse sea salt (SO4S and NITS).
! (bec, rjp, bmy, 4/20/04, 7/20/04, 10/25/05)
!
! Arguments as Input:
! ============================================================================
! (1 ) TC (REAL*8 ) : Tracer [kg]
! (2 ) N (INTEGER) : N=1 is SO4S; N=2 is NITS
!
! Arguments as Output:
! ============================================================================
! (1 ) TC (REAL*8 ) : Contains modified tracer
!
! NOTES:
! (1 ) Now references SALA_REDGE_um and SALC_REDGE_um from "tracer_mod.f"
! (bmy, 7/20/04)
! (2 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : T, BXHEIGHT, RH
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE PRESSURE_MOD, ONLY : GET_PCENTER
USE TRACER_MOD, ONLY : SALA_REDGE_um, SALC_REDGE_um, XNUMOL
USE TRACERID_MOD, ONLY : IDTSO4s, IDTNITs
USE TIME_MOD, ONLY : GET_ELAPSED_SEC, GET_TS_CHEM
USE GRID_MOD, ONLY : GET_AREA_CM2
# include "CMN_SIZE" ! Size parameters
# include "CMN_GCTM" ! g0
# include "CMN_DIAG" ! ND44
! Arguments
INTEGER, INTENT(IN) :: N
REAL*8, INTENT(INOUT) :: TC(IIPAR,JJPAR,LLPAR)
! Local variables
INTEGER :: I, J, L, DTCHEM
REAL*8 :: DELZ, DELZ1, REFF
REAL*8 :: P, DP, PDP, TEMP
REAL*8 :: CONST, SLIP, VISC, FAC1
REAL*8 :: FAC2, FLUX, AREA_CM2, RHB
REAL*8 :: RCM, RWET, RATIO_R, RHO
REAL*8 :: TOT1, TOT2
REAL*8 :: VTS(LLPAR)
REAL*8 :: TC0(LLPAR)
! Parameters
REAL*8, PARAMETER :: C1 = 0.7674d0
REAL*8, PARAMETER :: C2 = 3.079d0
REAL*8, PARAMETER :: C3 = 2.573d-11
REAL*8, PARAMETER :: C4 = -1.424d0
REAL*8, PARAMETER :: DEN = 2200.0d0 ! [kg/m3] sea-salt density
! Arrays
INTEGER :: IDDEP(2)
INTEGER :: IDTRC(2)
!=================================================================
! GRAV_SETTLING begins here!
!=================================================================
! Return if tracers are undefined
IF ( IDTSO4s == 0 .and. IDTNITs == 0 ) RETURN
! Return if it's the start of the run
IF ( GET_ELAPSED_SEC() == 0 ) RETURN
! Chemistry timestep [s]
DTCHEM = GET_TS_CHEM() * 60d0
! Store in IDDEP array
IDDEP(1) = DRYSO4s
IDDEP(2) = DRYNITs
! Tracer array
IDTRC(1) = IDTSO4s
IDTRC(2) = IDTNITs
! Coarse mode
REFF = 0.5d-6 * ( SALC_REDGE_um(1) + SALC_REDGE_um(2) )
! Sea salt radius [cm]
RCM = REFF * 100d0
! Exponential factors
FAC1 = C1 * ( RCM**C2 )
FAC2 = C3 * ( RCM**C4 )
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, VTS, P, TEMP, RHB, RWET )
!$OMP+PRIVATE( RATIO_R, RHO, DP, PDP, CONST, SLIP, VISC, TC0 )
!$OMP+PRIVATE( DELZ, DELZ1, TOT1, TOT2, AREA_CM2, FLUX )
!$OMP+SCHEDULE( DYNAMIC )
DO J = 1, JJPAR
DO I = 1, IIPAR
! Initialize
DO L = 1, LLPAR
VTS(L) = 0d0
ENDDO
! Loop over levels
DO L = 1, LLPAR
! Pressure at center of the level [kPa]
P = GET_PCENTER(I,J,L) * 0.1d0
! Temperature [K]
TEMP = T(I,J,L)
! Cap RH at 0.99
RHB = MIN( 0.99d0, RH(I,J,L) * 1d-2 )
! Aerosol growth with relative humidity in radius [m]
! (Gerber, 1985)
RWET = 0.01d0*(FAC1/(FAC2-DLOG(RHB))+RCM**3.d0)**0.33d0
! Ratio dry over wet radii at the cubic power
RATIO_R = ( REFF / RWET )**3.d0
! Density of the wet aerosol (kg/m3)
RHO = RATIO_R * DEN + ( 1.d0 - RATIO_R ) * 1000.d0
! Dp = particle diameter [um]
DP = 2.d0 * RWET * 1.d6
! PdP = P * dP [hPa * um]
PDp = P * Dp
! Constant
CONST = 2.d0 * RHO * RWET**2 * g0 / 9.d0
!===========================================================
! NOTE: Slip correction factor calculations following
! Seinfeld, pp464 which is thought to be more accurate
! but more computation required. (rjp, 1/24/02)
!
! # air molecule number density
! num = P * 1d3 * 6.023d23 / (8.314 * Temp)
!
! # gas mean free path
! lamda = 1.d6/( 1.41421 * num * 3.141592 * (3.7d-10)**2 )
!
! # Slip correction
! Slip = 1. + 2. * lamda * (1.257 + 0.4 * exp( -1.1 * Dp
! & / (2. * lamda))) / Dp
!
! NOTE: Eq) 3.22 pp 50 in Hinds (Aerosol Technology)
! which produces slip correction factore with small error
! compared to the above with less computation.
!===========================================================
! Slip correction factor (as function of P*dp)
Slip = 1.d0+(15.60d0 + 7.0d0 * EXP(-0.059d0 * PDp)) / PDp
! Viscosity [Pa*s] of air as a function of temperature
VISC = 1.458d-6 * (Temp)**(1.5d0) / ( Temp + 110.4d0 )
! Settling velocity [m/s]
VTS(L) = CONST * Slip / VISC
ENDDO
! Method is to solve bidiagonal matrix which is
! implicit and first order accurate in z (rjp, 1/24/02)
! Save initial tracer concentration in column
DO L = 1, LLPAR
TC0(L) = TC(I,J,L)
ENDDO
! We know the boundary condition at the model top
L = LLTROP
DELZ = BXHEIGHT(I,J,L)
TC(I,J,L) = TC(I,J,L) / ( 1.d0 + DTCHEM * VTS(L) / DELZ )
DO L = LLTROP-1, 1, -1
DELZ = BXHEIGHT(I,J,L)
DELZ1 = BXHEIGHT(I,J,L+1)
TC(I,J,L) = 1.d0 / ( 1.d0 + DTCHEM * VTS(L) / DELZ )
& * ( TC(I,J,L) + DTCHEM * VTS(L+1) / DELZ1
& * TC(I,J,L+1) )
ENDDO
!==============================================================
! ND44 diagnostic: sea salt loss [molec/cm2/s]
!==============================================================
IF ( ND44 > 0 ) THEN
! Initialize
TOT1 = 0d0
TOT2 = 0d0
! Compute column totals of TCO(:) and TC(I,J,:,N)
DO L = 1, LLPAR
TOT1 = TOT1 + TC0(L)
TOT2 = TOT2 + TC(I,J,L)
ENDDO
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Convert sea salt flux from [kg/s] to [molec/cm2/s]
FLUX = ( TOT1 - TOT2 ) / DTCHEM
FLUX = FLUX * XNUMOL(IDTRC(N)) / AREA_CM2
! Store in AD44 array
AD44(I,J,IDDEP(N),1) = AD44(I,J,IDDEP(N),1) + FLUX
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE GRAV_SETTLING
!------------------------------------------------------------------------------
SUBROUTINE CHEM_DMS 1,16
!
!******************************************************************************
! Subroutine CHEM_DMS is the DMS chemistry subroutine from Mian Chin's
! GOCART model, modified for use with the GEOS-CHEM model.
! (rjp, bdf, bmy, 5/31/00, 10/25/05)
!
! Module Variables used:
! ============================================================================
! (1 ) PSO2_DMS (REAL*8 ) : Array for P(SO2) from DMS [v/v]
! (2 ) PMSA_DMS (REAL*8 ) : Array for P(MSA) from DMS [v/v]
!
! Reaction List (by Mian Chin, chin@rondo.gsfc.nasa.gov)
! ============================================================================
!
! R1: DMS + OH -> a*SO2 + b*MSA OH addition channel
! k1 = { 1.7e-42*exp(7810/T)*[O2] / (1+5.5e-31*exp(7460/T)*[O2] }
! a = 0.75, b = 0.25
!
! R2: DMS + OH -> SO2 + ... OH abstraction channel
! k2 = 1.2e-11*exp(-260/T)
!
! DMS_OH = DMS0 * exp(-(r1+r2)* NDT1)
! where DMS0 is the DMS concentration at the beginning,
! r1 = k1*[OH], r2 = k2*[OH].
!
! R3: DMS + NO3 -> SO2 + ...
! k3 = 1.9e-13*exp(500/T)
!
! DMS = DMS_OH * exp(-r3*NDT1)
! where r3 = k3*[NO3].
!
! R4: DMS + X -> SO2 + ...
! assume to be at the rate of DMS+OH and DMS+NO3 combined.
!
! The production of SO2 and MSA here, PSO2_DMS and PMSA_DMS, are saved
! for use in CHEM_SO2 and CHEM_MSA subroutines as a source term. They
! are in unit of [v/v/timestep].
!
! NOTES:
! (1 ) Now reference AD, AIRDEN, and SUNCOS from "dao_mod.f". Added
! parallel DO-loops. Also now extract OH and NO3 from SMVGEAR
! for coupled chemistry-aerosol runs. (rjp, bdf, bmy, 9/16/02)
! (2 ) Bug fix: remove duplicate definition of RK3 (bmy, 3/23/03)
! (3 ) Now use function GET_TS_CHEM from "time_mod.f". (bmy, 3/27/03)
! (4 ) Now reference STT and ITS_A_FULLCHEM_SIM from "tracer_mod.f"
! Now replace IJSURF w/ an analytic function. (bmy, 7/20/04)
! (5 ) Shift rows 8,9 in AD05 to 9,10 in to make room for P(SO4) from O3
! oxidation in sea-salt aerosols (bec, bmy, 4/13/05)
! (6 ) Now remove reference to CMN, it's obsolete. Now reference
! ITS_IN_THE_STRAT from "tropopause_mod.f". (bmy, 8/22/05)
! (7 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! Reference to F90 modules
USE DAO_MOD, ONLY : AD, AIRDEN, SUNCOS, T
USE DIAG_MOD, ONLY : AD05
USE DRYDEP_MOD, ONLY : DEPSAV
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, ITS_A_FULLCHEM_SIM, XNUMOL
USE TRACERID_MOD, ONLY : IDTDMS
USE TROPOPAUSE_MOD, ONLY : ITS_IN_THE_STRAT
# include "CMN_SIZE" ! Size parameters
# include "CMN_GCTM" ! AIRMW
# include "CMN_DIAG" ! ND05, LD05
! Local variables
LOGICAL :: IS_FULLCHEM
INTEGER :: I, J, L, IJLOOP
REAL*8 :: TK, O2, RK1, RK2, RK3, F
REAL*8 :: DMS, DMS0, DMS_OH, DTCHEM, XOH, XN3
REAL*8 :: XX, OH, OH0, XNO3, XNO30, LOH
REAL*8 :: LNO3
! Parameters
REAL*8, PARAMETER :: FX = 1.0d0
REAL*8, PARAMETER :: A = 0.75d0
REAL*8, PARAMETER :: B = 0.25d0
! From D4: only 0.8 efficiency, also some goes to DMSO and lost.
! So we assume 0.75 efficiency for DMS addtion channel to form
! products.
REAL*8, PARAMETER :: EFF = 1d0
! External functions
REAL*8, EXTERNAL :: BOXVL
!=================================================================
! CHEM_DMS begins here!
!=================================================================
IF ( IDTDMS == 0 ) RETURN
! Flag for fullchem simulation
IS_FULLCHEM = ITS_A_FULLCHEM_SIM()
! DTCHEM is the chemistry timestep in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Factor to convert AIRDEN from kgair/m3 to molecules/cm3:
f = 1000.d0 / AIRMW * 6.022d23 * 1.d-6
!=================================================================
! Do the chemistry over all tropospheric grid boxes!
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, IJLOOP, J, L, TK, O2, DMS0, OH, XNO3, RK1, RK2 )
!$OMP+PRIVATE( RK3, DMS_OH, DMS, OH0, XNO30, XOH, XN3, XX, LOH, LNO3 )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
! Skip stratospheric boxes
IF ( ITS_IN_THE_STRAT( I, J, L ) ) CYCLE
! IJLOOP is the 1-D grid box index for SUNCOS
IJLOOP = ( (J-1) * IIPAR ) + I
! Temperature [K]
TK = T(I,J,L)
! Get O2 [molec/cm3], DMS [v/v], OH [molec/cm3], NO3 [molec/cm3]
O2 = AIRDEN(L,I,J) * f * 0.21d0
DMS0 = STT(I,J,L,IDTDMS)
OH = GET_OH( I, J, L )
XNO3 = GET_NO3( I, J, L )
!==============================================================
! (1) DMS + OH: RK1 - addition channel
! RK2 - abstraction channel
!==============================================================
RK1 = 0.d0
RK2 = 0.d0
RK3 = 0.d0
IF ( OH > 0.d0 ) THEN
RK1 = ( 1.7d-42 * EXP( 7810.d0 / TK ) * O2 ) /
& ( 1.d0 + 5.5d-31 * EXP( 7460.d0 / TK ) * O2 ) * OH
RK2 = 1.2d-11 * EXP( -260.d0 / TK ) * OH
ENDIF
!==============================================================
! (2) DMS + NO3 (only happens at night):
!==============================================================
IF ( SUNCOS(IJLOOP) <= 0d0 ) THEN
RK3 = 1.9d-13 * EXP( 500.d0 / TK ) * XNO3
ENDIF
!==============================================================
! Update DMS concentrations after reaction with OH and NO3,
! and also account for DMS + X assuming at a rate as
! (DMS+OH)*Fx in the day and (DMS+NO3)*Fx at night:
!
! DMS_OH : DMS concentration after reaction with OH
! DMS : DMS concentration after reaction with NO3
! (min(DMS) = 1.0E-32)
!
! NOTE: If we are doing a coupled fullchem/aerosol run, then
! also modify OH and NO3 concentrations after rxn w/ DMS.
!==============================================================
DMS_OH = DMS0 * EXP( -( RK1 + RK2 ) * Fx * DTCHEM )
DMS = DMS_OH * EXP( -( RK3 ) * Fx * DTCHEM )
IF ( DMS < SMALLNUM ) DMS = 0d0
! Archive initial OH and NO3 for diagnostics
OH0 = OH
XNO30 = XNO3
IF ( IS_FULLCHEM ) THEN
! Update OH after rxn w/ DMS (coupled runs only)
OH = OH0 - ( ( DMS0 - DMS_OH ) * AIRDEN(L,I,J) * f )
IF ( OH < SMALLNUM ) OH = 0d0
! Update NO3 after rxn w/ DMS (coupled runs only)
XNO3 = XNO30 - ( ( DMS_OH - DMS ) * AIRDEN(L,I,J) * f )
IF ( XNO3 < SMALLNUM ) XNO3 = 0d0
ENDIF
! Save DMS back to the tracer array
STT(I,J,L,IDTDMS) = DMS
!==============================================================
! Save SO2 and MSA production from DMS oxidation
! in [mixing ratio/timestep]:
!
! SO2 is formed in DMS+OH addition (0.85) and abstraction
! (1.0) channels as well as DMS + NO3 reaction. We also
! assume that SO2 yield from DMS + X is 1.0.
!
! MSA is formed in DMS + OH addition (0.15) channel.
!==============================================================
IF ( ( RK1 + RK2 ) == 0.d0 ) THEN
PMSA_DMS(I,J,L) = 0.d0
ELSE
PMSA_DMS(I,J,L) = ( DMS0 - DMS_OH ) *
& B*RK1 / ( ( RK1 + RK2 ) * Fx ) * EFF
ENDIF
PSO2_DMS(I,J,L) = DMS0 - DMS - PMSA_DMS(I,J,L)
!==============================================================
! ND05 diagnostic: production and loss
!
! For the offline run, we are reading in monthly mean OH, NO3
! from disk. We don't modify these, so LOH = 0 and LNO3 = 0.
!==============================================================
IF ( ND05 > 0 .and. L <= LD05 ) THEN
! P(SO2) from DMS+OH, DMS+NO3, and DMS+X
XOH = ( DMS0 - DMS_OH ) / Fx * AD(I,J,L) / TCVV_S
XN3 = ( DMS_OH - DMS ) / Fx * AD(I,J,L) / TCVV_S
XX = ( ( DMS0 - DMS ) * AD(I,J,L) / TCVV_S ) - XOH - XN3
! Convert L(OH) and L(NO3) from [molec/cm3] to [kg/timestep]
LOH = ( OH0 - OH ) * BOXVL(I,J,L) / XNUMOL_OH
LNO3 = ( XNO30 - XNO3 ) * BOXVL(I,J,L) / XNUMOL_NO3
! Store P(SO2) from DMS + OH [kg S/timestep]
AD05(I,J,L,1) = AD05(I,J,L,1) + XOH
! Store P(SO2) from DMS + NO3 [kg S/timestep]
AD05(I,J,L,2) = AD05(I,J,L,2) + XN3
! Store total P(SO2) from DMS [kg S/timestep]
AD05(I,J,L,3) = AD05(I,J,L,3) +
& ( PSO2_DMS(I,J,L) * AD(I,J,L) / TCVV_S )
! Store P(MSA) from DMS [kg S/timestep]
AD05(I,J,L,4) = AD05(I,J,L,4) +
& ( PMSA_DMS(I,J,L) * AD(I,J,L) / TCVV_S )
! Store L(OH) by DMS [kg OH/timestep]
AD05(I,J,L,9) = AD05(I,J,L,9) + LOH
! Store L(NO3) by DMS [kg NO3/timestep]
AD05(I,J,L,10) = AD05(I,J,L,10) + LNO3
ENDIF
!==============================================================
! For a coupled fullchem/aerosol run, save OH [molec/cm3]
! and NO3 [molec/cm3] back into the CSPEC array of SMVGEAR
!==============================================================
IF ( IS_FULLCHEM ) THEN
CALL SET_OH( I, J, L, OH )
CALL SET_NO3( I, J, L, XNO3 )
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE CHEM_DMS
!------------------------------------------------------------------------------
SUBROUTINE CHEM_H2O2 1,25
!
!******************************************************************************
! Subroutine CHEM_H2O2 is the H2O2 chemistry subroutine for offline sulfate
! simulations. For coupled runs, H2O2 chemistry is already computed by
! the SMVGEAR module. (rjp, bmy, 11/26/02, 10/25/05)
!
! NOTES:
! (1 ) Bug fix: need to multiply DXYP by 1d4 for cm2 (bmy, 3/23/03)
! (2 ) Now replace DXYP(JREF)*1d4 with routine GET_AREA_CM2 of "grid_mod.f"
! Now use functions GET_MONTH and GET_TS_CHEM from "time_mod.f".
! (bmy, 3/27/03)
! (3 ) Now references PBLFRAC from "drydep_mod.f". Now apply dry deposition
! throughout the entire PBL. Added FREQ variable. (bmy, 8/1/03)
! (4 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (5 ) Now use diurnally-varying JO1D. Now use new unit conversion for
! the ND44 diagnostic. (rjp, bmy, 3/30/04)
! (6 ) Now use parallel DO-loop to zero ND44_TMP. Now uses ITS_A_NEW_MONTH
! from time_mod.f. (bmy, 4/14/04)
! (7 ) Now reference STT & TCVV from "tracer_mod.f". Also replace IJSURF
! with an analytic function. Now references DATA_DIR from
! "directory_mod.f". (bmy, 7/20/04)
! (8 ) Now suppress output from READ_BPCH with QUIET keyword (bmy, 1/25/05)
! (9 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP
! from "pbl_mix_mod.f" (bmy, 2/22/05)
! (10) Now read offline files from "sulfate_sim_200508/offline". Now remove
! reference to CMN, it's obsolete. Now reference ITS_IN_THE_STRAT from
! "tropopause_mod.f". (bmy, 8/22/05)
! (11) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
! (12) Now references XNUMOL 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, AIRDEN, OPTD, SUNCOS, T
USE DIAG_MOD, ONLY : AD44
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP
USE TIME_MOD, ONLY : GET_MONTH, GET_TS_CHEM, ITS_A_NEW_MONTH
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTH2O2
USE TRANSFER_MOD, ONLY : TRANSFER_3D_TROP
USE UVALBEDO_MOD, ONLY : UVALBEDO
USE TROPOPAUSE_MOD, ONLY : ITS_IN_THE_STRAT
# include "cmn_fj.h" ! IPAR, JPAR, LPAR, CMN_SIZE
# include "CMN_GCTM" ! AIRMW
# include "CMN_DIAG" ! ND44
! Local variables
LOGICAL :: FIRST = .TRUE.
INTEGER, SAVE :: LASTMONTH = -99
INTEGER :: I, J, L, JLOOP
REAL*4 :: ARRAY(IGLOB,JGLOB,LLTROP)
REAL*8 :: ND44_TMP(IIPAR,JJPAR,LLTROP)
REAL*8 :: DT, Koh, DH2O2, M, F , XTAU
REAL*8 :: H2O20, H2O2, ALPHA, FLUX, FREQ, PHOTJ
REAL*8, PARAMETER :: A = 2.9d-12
CHARACTER(LEN=255) :: FILENAME
!=================================================================
! CHEM_H2O2 begins here!
!=================================================================
IF ( IDTH2O2 == 0 .or. DRYH2O2 == 0 ) RETURN
! Chemistry timestep [s]
DT = GET_TS_CHEM() * 60d0
! Factor to convert AIRDEN from kgair/m3 to molecules/cm3:
F = 1000.d0 / AIRMW * 6.022d23 * 1.d-6
! Zero ND44_TMP array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
ND44_TMP(I,J,L) = 0d0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!=================================================================
! For offline run: read J(H2O2) from disk below
!=================================================================
IF ( ITS_A_NEW_MONTH() ) THEN
! File name to read data
FILENAME = TRIM( DATA_DIR ) //
& 'sulfate_sim_200508/offline/JH2O2.' //
& GET_NAME_EXT() // '.' // GET_RES_EXT()
! Print filename
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - CHEM_H2O2: Reading ', a )
! Get TAU0 value for this month in "generic" year 1985
XTAU = GET_TAU0( GET_MONTH(), 1, 1985 )
! Read J(H2O2) [s-1] from disk (only up to tropopause)
! limit array 3d dimension to LLTROP_FIX, i.e, case of annual mean
! tropopause. This is backward compatibility with
! offline data set.
CALL READ_BPCH2( FILENAME, 'JV-MAP-$', 3,
& XTAU, IGLOB, JGLOB,
& LLTROP_FIX, ARRAY(:,:,1:LLTROP_FIX), QUIET=.TRUE. )
! & XTAU, IGLOB, JGLOB,
! & LLTROP, ARRAY, QUIET=.TRUE. )
! Cast to REAL*8 and resize if necessary
CALL TRANSFER_3D_TROP( ARRAY, JH2O2 )
! Reset LASTMONTH
!LASTMONTH = GET_MONTH()
ENDIF
!=================================================================
! Loop over tropopsheric grid boxes and do chemistry
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, M, H2O20, KOH, FREQ, ALPHA, DH2O2, H2O2, FLUX )
!$OMP+PRIVATE( JLOOP, PHOTJ )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
! Initialize for safety's sake
FLUX = 0d0
FREQ = 0d0
! Skip stratospheric boxes
IF ( ITS_IN_THE_STRAT( I, J, L ) ) CYCLE
! Density of air [molec/cm3]
M = AIRDEN(L,I,J) * f
! Initial H2O2 [v/v]
H2O20 = STT(I,J,L,IDTH2O2)
! Loss frequenty due to OH oxidation [s-1]
KOH = A * EXP( -160.d0 / T(I,J,L) ) * GET_OH(I,J,L)
! H2O2 drydep frequency [1/s]. Account for the fraction
! of grid box (I,J,L) that is located beneath the PBL top.
FREQ = DEPSAV(I,J,DRYH2O2) * GET_FRAC_UNDER_PBLTOP( I, J, L )
! 1-D grid box index for SUNCOS
JLOOP = ( (J-1) * IIPAR ) + I
! Impose a diurnal variation of jH2O2 by multiplying COS of
! solar zenith angle normalized by maximum solar zenith angle
! because the archived JH2O2 is for local noon time
IF ( COSZM(I,J) > 0.d0 ) THEN
PHOTJ = JH2O2(I,J,L) * SUNCOS(JLOOP) / COSZM(I,J)
PHOTJ = MAX( PHOTJ, 0d0 )
ELSE
PHOTJ = 0d0
ENDIF
! Compute loss fraction from OH, photolysis, drydep [unitless].
ALPHA = 1.D0 + ( KOH + PHOTJ + FREQ ) * DT
! Delta H2O2 [v/v]
DH2O2 = PH2O2m(I,J,L) * DT / ( ALPHA * M )
! Final H2O2 [v/v]
H2O2 = ( H2O20 / ALPHA + DH2O2 )
IF ( H2O2 < SMALLNUM ) H2O2 = 0d0
! Store final H2O2 in STT
STT(I,J,L,IDTH2O2) = H2O2
!==============================================================
! ND44 diagnostics: H2O2 drydep loss [molec/cm2/s]
!==============================================================
IF ( ND44 > 0 .AND. FREQ > 0d0 ) THEN
! Convert H2O2 from [v/v] to H2O2 [molec/cm2/s]
FLUX = H2O20 * FREQ * DT / ( 1.D0 + FREQ * DT )
FLUX = FLUX * AD(I,J,L) / TCVV(IDTH2O2)
FLUX = FLUX * XNUMOL(IDTH2O2) / ( GET_AREA_CM2( J ) * DT )
! Save dryd flx in ND44_TMP as a placeholder
ND44_TMP(I,J,L) = ND44_TMP(I,J,L) + FLUX
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, LLTROP
AD44(I,J,DRYH2O2,1) = AD44(I,J,DRYH2O2,1) + ND44_TMP(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_H2O2
!------------------------------------------------------------------------------
SUBROUTINE CHEM_SO2 1,31
!
!******************************************************************************
! Subroutine CHEM_SO2 is the SO2 chemistry subroutine
! (rjp, bmy, 11/26/02, 10/25/05)
!
! Module variables used:
! ============================================================================
! (1 ) PSO2_DMS (REAL*8 ) : Array for P(SO2) from DMS [v/v/timestep]
! (2 ) PSO4_SO2 (REAL*8 ) : Array for P(SO4) from SO2 [v/v/timestep]
! (3 ) LSO2_AQ (REAL*8 ) : Array for L(SO2) from Aqueuos chem [v/v/timestep]
!
! Reaction List (by Rokjin Park, rjp@io.harvard.edu)
! ============================================================================
! (1 ) SO2 production:
! DMS + OH, DMS + NO3 (saved in CHEM_DMS)
!
! (2 ) SO2 loss:
! (a) SO2 + OH -> SO4
! (b) SO2 -> drydep
! (c) SO2 + H2O2 or O3 (aq) -> SO4
!
! (3 ) SO2 = SO2_0 * exp(-bt) + PSO2_DMS/bt * [1-exp(-bt)]
!
! where b is the sum of the reaction rate of SO2 + OH and the dry
! deposition rate of SO2, PSO2_DMS is SO2 production from DMS in
! MixingRatio/timestep.
!
! If there is cloud in the gridbox (fraction = fc), then the aqueous
! phase chemistry also takes place in cloud. The amount of SO2 oxidized
! by H2O2 in cloud is limited by the available H2O2; the rest may be
! oxidized due to additional chemistry, e.g, reaction with O3 or O2
! (catalyzed by trace metal).
!
! NOTES:
! (1 ) Removed duplicate definition of Ki (bmy, 11/15/01)
! (2 ) Eliminate duplicate HPLUS definition. Make adjustments to facilitate
! SMVGEAR chemistry for fullchem runs (rjp, bmy, 3/23/03)
! (3 ) Now replace DXYP(J+J0)*1d4 with routine GET_AREA_CM2 of "grid_mod.f"
! Now use function GET_TS_CHEM from "time_mod.f".
! (4 ) Now apply dry deposition to entire PBL. Now references PBLFRAC array
! from "drydep_mod.f". (bmy, 8/1/03)
! (5 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (6 ) Now use parallel DO-loop to zero ND44_TMP (bmy, 4/14/04)
! (7 ) Now reference STT, TCVV, & ITS_AN_AEROSOL_SIM from "tracer_mod.f".
! Now reference DATA_DIR from "directory_mod.f" (bmy, 7/20/04)
! (8 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f" (bmy, 2/22/05)
! (9 ) Modified for SO4s, NITs. Also modified for alkalinity w/in the
! seasalt chemistry. (bec, bmy, 4/13/05)
! (10) Now remove reference to CMN, it's obsolete. Now reference
! ITS_IN_THE_STRAT from "tropopause_mod.f" (bmy, 8/22/05)
! (11) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! Reference to diagnostic arrays
USE DAO_MOD, ONLY : AD, AIRDEN, T
USE DIAG_MOD, ONLY : AD05, AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE GLOBAL_HNO3_MOD, ONLY : GET_GLOBAL_HNO3
USE GRID_MOD, ONLY : GET_AREA_CM2
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP
USE PRESSURE_MOD, ONLY : GET_PCENTER
USE TIME_MOD, ONLY : GET_TS_CHEM, GET_MONTH
USE TIME_MOD, ONLY : ITS_A_NEW_MONTH
USE TRACER_MOD, ONLY : STT, TCVV, ITS_AN_AEROSOL_SIM
USE TRACER_MOD, ONLY : XNUMOL
USE TRACERID_MOD, ONLY : IDTH2O2, IDTSO2
USE SEASALT_MOD, ONLY : GET_ALK
USE WETSCAV_MOD, ONLY : H2O2s, SO2s
USE TROPOPAUSE_MOD, ONLY : ITS_IN_THE_STRAT
# include "CMN_SIZE" ! Size parameters
# include "CMN_GCTM" ! AIRMW
# include "CMN_DIAG" ! LD05, ND05, ND44
! Local variables
LOGICAL :: IS_OFFLINE
INTEGER :: I, J, L, I1, I2
INTEGER :: II, NSTEP
REAL*8 :: K0, Ki, KK, M, L1
REAL*8 :: L2, L3, Ld, F, Fc
REAL*8 :: RK, RKT, DTCHEM, DT_T, TK
REAL*8 :: F1, RK1, RK2, RK3, SO20
REAL*8 :: SO2_cd, H2O20, O3, L2S, L3S
REAL*8 :: LWC, KaqH2O2, KaqO3, PATM, FLUX
REAL*8 :: ALK, ALK1, ALK2, SO2_ss
REAL*8 :: Kt1, Kt2, AREASS1, AREASS2
REAL*8 :: PSO4E, PSO4F, Kt1N, Kt2N
REAL*8 :: AREA_CM2
REAL*8 :: ND44_TMP(IIPAR,JJPAR,LLTROP)
! Parameters
REAL*8, PARAMETER :: HPLUS = 3.16227766016837953d-5 !pH = 4.5
REAL*8, PARAMETER :: MINDAT = 1.d-20
!=================================================================
! CHEM_SO2 begins here!
!=================================================================
IF ( IDTH2O2 == 0 .or. IDTSO2 == 0 .or. DRYSO2 == 0 ) RETURN
! Is it an offline simulation?
IS_OFFLINE = ITS_AN_AEROSOL_SIM()
! Read HNO3 for offline simulation
IF ( IS_OFFLINE ) THEN
IF ( ITS_A_NEW_MONTH() ) THEN
CALL GET_GLOBAL_HNO3( GET_MONTH() )
ENDIF
ENDIF
! DTCHEM is the chemistry timestep in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Factor to convert AIRDEN from [kg air/m3] to [molec air/cm3]
F = 1000.d0 / AIRMW * 6.022d23 * 1.d-6
Ki = 1.5d-12
! Zero ND44_TMP array
ND44_TMP = 0d0
! Loop over tropospheric grid boxes
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, SO20, H2O20, O3, PATM, TK, K0, M, KK, F1, RK1 )
!$OMP+PRIVATE( RK2, RK, RKT, SO2_cd, L1, Ld, L2, L2S, L3, L3S, FC, LWC )
!$OMP+PRIVATE( KaqH2O2, KaqO3, AREA_CM2, FLUX, ALK, ALK1, ALK2 )
!$OMP+PRIVATE( Kt1, Kt2, AREASS1, AREASS2, SO2_ss, Kt1N, Kt2N )
!$OMP+PRIVATE( PSO4E, PSO4F )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
! Initialize for safety's sake
AREA_CM2 = 0d0
FLUX = 0d0
Ld = 0d0
! Skip stratospheric boxes
IF ( ITS_IN_THE_STRAT( I, J, L ) ) CYCLE
! Initial SO2, H2O2 and O3 [v/v]
SO20 = STT(I,J,L,IDTSO2)
H2O20 = STT(I,J,L,IDTH2O2)
O3 = GET_O3(I,J,L)
! PATM : Atmospheric pressure in atm
PATM = GET_PCENTER( I, J, L ) / 1013.25d0
! TK : Temperature [K]
TK = T(I,J,L)
IF ( IS_OFFLINE ) THEN
! Gas phase SO4 production is done here in offline run only
! RK1: SO2 + OH(g) [s-1] (rjp, bmy, 3/23/03)
K0 = 3.0d-31 * ( 300.d0 / TK )**3.3d0
M = AIRDEN(L,I,J) * F
KK = K0 * M / Ki
F1 = ( 1.d0 + ( LOG10( KK ) )**2 )**( -1 )
RK1 = ( K0 * M / ( 1.d0 + KK ) ) * 0.6d0**F1 * GET_OH(I,J,L)
ELSE
! For online runs, SMVGEAR deals w/ this computation,
! so we can simply set RK1 = 0 (rjp, bmy, 3/23/03)
K0 = 0.d0
M = 0.d0
KK = 0.d0
F1 = 0.d0
RK1 = 0.d0
ENDIF
! SO2 drydep frequency [1/s]. Also accounts for the fraction
! of grid box (I,J,L) that is located beneath the PBL top.
RK2 = DEPSAV(I,J,DRYSO2) * GET_FRAC_UNDER_PBLTOP( I, J, L )
! RK: total reaction rate [1/s]
RK = ( RK1 + RK2 )
! RKT: RK * DTCHEM [unitless] (bmy, 6/1/00)
RKT = RK * DTCHEM
!==============================================================
! Update SO2 conc. after gas phase chemistry and deposition
!==============================================================
IF ( RK > 0.d0 ) THEN
SO2_cd = ( SO20 * EXP( -RKT ) ) +
& ( PSO2_DMS(I,J,L) * ( 1.d0 - EXP( -RKT ) ) / RKT )
L1 = ( SO20 - SO2_cd + PSO2_DMS(I,J,L) ) * RK1/RK
Ld = ( SO20 - SO2_cd + PSO2_DMS(I,J,L) ) * RK2/RK
ELSE
SO2_cd = SO20
L1 = 0.d0
ENDIF
!==============================================================
! Update SO2 conc. after seasalt chemistry (bec, 12/7/04)
!==============================================================
! Get alkalinity of accum (ALK1) and coarse (ALK2) [kg]
CALL GET_ALK( I, J, L, ALK1, ALK2, Kt1, Kt2, Kt1N, Kt2N )
! Total alkalinity [kg]
ALK = ALK1 + ALK2
! If (1) there is alkalinity, (2) there is SO2 present, and
! (3) O3 is in excess, then compute seasalt SO2 chemistry
IF ( ( ALK > MINDAT ) .AND.
& ( SO2_cd > MINDAT ) .AND.
& ( SO2_cd < O3 ) ) THEN
! Compute oxidation of SO2 -> SO4 and condensation of
! HNO3 -> nitrate within the seasalt aerosol
CALL SEASALT_CHEM( I, J, L, ALK1, ALK2,
& SO2_cd, Kt1, Kt2, Kt1N, Kt2N,
& SO2_ss, PSO4E, PSO4F )
ELSE
! Otherwise set equal to zero
SO2_ss = SO2_cd
PSO4E = 0.d0
PSO4F = 0.d0
PNITS(I,J,L) = 0.d0
ENDIF
!==============================================================
! Update SO2 concentration after cloud chemistry
! SO2 chemical loss rate = SO4 production rate [v/v/timestep]
!==============================================================
! Volume cloud fraction (Sundqvist et al 1989) [unitless]
FC = VCLDF(I,J,L)
! Liquid water content in cloudy area of grid box [m3/m3]
LWC = GET_LWC( TK ) * FC
! Zero variables
KaqH2O2 = 0.d0
KaqO3 = 0.d0
L2 = 0.d0
L3 = 0.d0
L2S = 0.d0
L3S = 0.d0
! If (1) there is cloud, (2) there is SO2 present, and
! (3) the T > -15 C, then compute aqueous SO2 chemistry
IF ( ( FC > 0.d0 ) .AND.
& ( SO2_ss > MINDAT ) .AND.
& ( TK > 258.0 ) ) THEN
!===========================================================
! NOTE...Sulfate production from aquatic reactions of SO2
! with H2O2 & O3 is computed here and followings are
! approximations or method used for analytical (integral)
! solution of these computations. Please email us
! (rjp@io.harvard.edu or bmy@io.harvard.edu) if you find
! anything wrong or questionable.
!
! 1) with H2O2(aq)
! [HSO3-] + [H+] + [H2O2(aq)] => [SO4=] (rxn)
! d[SO4=]/dt = k[H+][HSO3-][H2O2(aq)] (M/s) (rate)
!
! we can rewrite k[H+][HSO3-] as K1 pSO2 hSO2,
! where pSO2 is equilibrium vapor pressure of SO2(g)
! in atm, and hSO2 is henry's law constant for SO2
!
! Therefore, rate can be written as
!
! k * K1 * pSO2 * hSO2 * pH2O2 * hH2O2,
!
! where pH2O2 is equilibrium vapor pressure of H2O2(g),
! and hH2O2 is henry's law constant for H2O2. Detailed
! values are given in AQCHEM_SO2 routine.
!
! Let us define a fraction of gas phase of A species
! in equilibrium with aqueous phase as
!
! xA = 1/(1+f),
!
! where f = hA * R * T * LWC,
! hA = Henry's constant,
! R = gas constant,
! T = temperature in kelvin,
! LWC = liquid water content [m3/m3]
!
! As a result, the rate would become:
!
! d[SO4=]
! ------- = k K1 hSO2 hH2O2 xSO2 xH2O2 P P [SO2][H2O2]
! dt
! ^ ^ ^ ^ ^
! | |____________________________| | |
!
! mole/l/s mole/l/s v/v v/v
!
!
! And we multiply rate by (LWC * R * T / P) in order to
! convert unit from mole/l/s to v/v/s
!
! Finally we come to
!
! d[SO4=]
! ------- = KaqH2O2 [SO2][H2O2],
! dt
!
! where
!
! KaqH2O2 = k K1 hSO2 hH2O2 xSO2 xH2O2 P LWC R T,
!
! this new rate corresponds to a typical second order
! reaction of which analytical (integral) solution is
!
! X = A0 B0 ( exp[(A0-B0) Ka t] - 1 )
! / ( A0 exp[(A0-B0) Ka t] - B0 )
!
! inserting variables into solution then we get
! [SO4=] = [SO2][H2O2](exp[([SO2]-[H2O2]) KaqH2O2 t] - 1 )
! / ( [SO2] exp[([SO2]-[H2O2]) KaqH2O2 t] - [H2O2] )
!
! Note...Exactly same method can be applied to O3 reaction
! in aqueous phase with different rate constants.
!===========================================================
! Compute aqueous rxn rates for SO2
CALL AQCHEM_SO2( LWC, TK, PATM, SO2_ss, H2O20,
& O3, HPLUS, KaqH2O2, KaqO3 )
! Aqueous phase SO2 loss rate (v/v/timestep):
L2 = EXP( ( SO2_ss - H2O20 ) * KaqH2O2 * DTCHEM )
L3 = EXP( ( SO2_ss - O3 ) * KaqO3 * DTCHEM )
! Loss by H2O2
L2S = SO2_ss * H2O20 * (L2 - 1.D0) / ((SO2_ss * L2) - H2O20)
! Loss by O3
L3S = SO2_ss * O3 * (L3 - 1.D0) / ((SO2_ss * L3) - O3)
SO2_ss = MAX( SO2_ss - ( L2S + L3S ), MINDAT )
H2O20 = MAX( H2O20 - L2S, MINDAT )
! Update SO2 level, save SO2[ppv], H2O2[ppv] for WETDEP
SO2s( I,J,L) = SO2_ss
H2O2s(I,J,L) = H2O20
ELSE
! Otherwise, don't do aqueous chemistry, and
! save the original concentrations into SO2 and H2O2
H2O2s(I,J,L) = MAX( H2O20, 1.0d-32 )
SO2s(I,J,L ) = MAX( SO2_ss, 1.0d-32 )
L2S = 0.d0
L3S = 0.d0
ENDIF
! Store updated SO2, H2O2 back to the tracer arrays
STT(I,J,L,IDTSO2) = SO2s( I,J,L)
STT(I,J,L,IDTH2O2) = H2O2s(I,J,L)
! SO2 chemical loss rate = SO4 production rate [v/v/timestep]
PSO4_SO2(I,J,L) = L1 + L2S + L3S + PSO4E
PSO4_ss (I,J,L) = PSO4F
!=================================================================
! ND05 Diagnostics [kg S/timestep]
!=================================================================
IF ( ND05 > 0 .and. L <= LD05 ) THEN
! P(SO4) from gas-phase oxidation [kg S/timestep]
AD05(I,J,L,5) = AD05(I,J,L,5) +
& ( L1 * AD(I,J,L) / TCVV_S )
! P(SO4) from aqueous-phase oxidation with H2O2 [kg S/timestep]
AD05(I,J,L,6) = AD05(I,J,L,6) +
& ( L2S * AD(I,J,L) / TCVV_S )
! P(SO4) from aqueous-phase oxidation with O3 [kg S/timestep]
AD05(I,J,L,7) = AD05(I,J,L,7) +
& ( L3S * AD(I,J,L) / TCVV_S )
! P(SO4) from O3 oxidation in sea-salt aerosols [kg S/timestep]
AD05(I,J,L,8) = AD05(I,J,L,8) +
& ( (PSO4E + PSO4F) * AD(I,J,L) / TCVV_S )
ENDIF
!=================================================================
! ND44 Diagnostic: Drydep flux of SO2 [molec/cm2/s]
!=================================================================
IF ( ND44 > 0 .AND. Ld > 0d0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Convert [v/v/timestep] to [molec/cm2/s]
FLUX = Ld * AD(I,J,L) / TCVV(IDTSO2)
FLUX = FLUX * XNUMOL(IDTSO2) / AREA_CM2 / DTCHEM
! Store dryd flx in ND44_TMP as a placeholder
ND44_TMP(I,J,L) = ND44_TMP(I,J,L) + FLUX
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, LLTROP
AD44(I,J,DRYSO2,1) = AD44(I,J,DRYSO2,1) + ND44_TMP(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_SO2
!------------------------------------------------------------------------------
SUBROUTINE SEASALT_CHEM ( I, J, L, ALK1, ALK2, 1,15
& SO2_cd, Kt1, Kt2, Kt1N, Kt2N,
& SO2_ss, PSO4E, PSO4F )
!
!******************************************************************************
! Function SEASALT_CHEM computes SO4 formed from S(IV) + O3 on seasalt
! aerosols as a function of seasalt alkalinity. (bec, bmy, 4/13/05, 12/8/06)
!
! Arguments as Input:
! ============================================================================
! (1 ) I (INTEGER) :
! (2 ) J (INTEGER) :
! (3 ) L (INTEGER) :
! (4 ) O30 (REAL*8) : Initial O3 mixing ratio (v/v]
! (5 ) ALK (REAL*8) : Alkalinity [kg] from seasalt from seasalt_mod
! (6 ) SO2_cd (REAL*8) : SO2 mixing ratio [v/v] after gas phase chemistry
! and dry deposition
! (7 ) Kt1 (REAL*8) : Rate constant [s-1] for sulfate formation on
! fine sea-salt aerosols from GET_ALK
! (8 ) Kt2 (REAL*8) : Rate constant [s-1] for sulfate formation on
! coarse sea-salt aerosols from GET_ALK
!
! Arguments as Output:
! ============================================================================
! (9 ) SO2_ss (REAL*8) : SO2 mixing ratio [v/v] updated after SS chemistry
! (10) SO4E (REAL*8) : SO4E (sulfate produced by S(IV)+O3 on fine seasalt)
! mixing ratio [v/v]
! (11) SO4F (REAL*8) : SO4F(sulfate produced by S(IV)+O3 on coarse seasalt)
! mixing ratio [v/v]
! (12) O3 (REAL*8) : Updated O3 mixing ratio [v/v] for fullchem runs
! only. Otherwise O30=O3.
!
! Chemical reactions:
! ============================================================================
! (R1) SO2 + O3 + ALK => SO4 + O2
! Modeled after Chamedies and Stelson, 1992?
!
! NOTES:
! (1 ) Now references XNUMOLAIR from "tracer_mod.f" (bmy, 10/25/05)
! (2 ) Bug fix: now avoid seg fault error if IDTHNO3 is zero, as it would
! be for an offline aerosol simulation. (bmy, 3/29/06)
! (3 ) Fixed typo in FALK_A_SO2 equation: C_FLUX_C should be C_FLUX_A.
! (havala, bec, bmy, 12/8/06)
!******************************************************************************
!
! References to F90 modules
USE COMODE_MOD, ONLY : CSPEC, JLOP, VOLUME
USE DAO_MOD, ONLY : AD, AIRDEN, AIRVOL
USE TRACERID_MOD
!----------------------------------------------------------------
! DIAGNOSTIC -- leave commented out for now (bec, bmy, 4/13/05)
!USE DIAG_MOD, ONLY : AD09
!----------------------------------------------------------------
USE ERROR_MOD, ONLY : GEOS_CHEM_STOP
USE TIME_MOD, ONLY : GET_TS_CHEM, GET_ELAPSED_SEC
USE ERROR_MOD, ONLY : IT_IS_NAN
USE TRACER_MOD, ONLY : ITS_A_FULLCHEM_SIM, STT
USE TRACER_MOD, ONLY : TCVV, XNUMOLAIR
USE GLOBAL_HNO3_MOD, ONLY : GET_HNO3_UGM3
USE TIME_MOD, ONLY : GET_ELAPSED_SEC, GET_MONTH
USE TIME_MOD, ONLY : ITS_A_NEW_MONTH
USE ISOROPIA_MOD, ONLY : GET_GNO3
# include "CMN_SIZE" ! Size parameters
!---------------------------------------------------------------
! DIAGNOSTIC -- leave commented out for now (bec, bmy, 4/13/05)
!# include "CMN_DIAG" ! ND19
!---------------------------------------------------------------
# include "CMN_GCTM" ! AIRMW
! Arguments
INTEGER, INTENT(IN) :: I, J, L
REAL*8, INTENT(IN) :: SO2_cd, Kt1, Kt2, Kt1N, Kt2N
REAL*8, INTENT(IN) :: ALK1, ALK2
REAL*8, INTENT(OUT) :: SO2_ss, PSO4E, PSO4F
! Local variables
INTEGER :: JLOOP
REAL*8 :: SO2_chem, DTCHEM
REAL*8 :: O3_cspec, O3_lost
REAL*8 :: EQ_1_C, EQ_2_C
REAL*8 :: SO4E, SO2_new, SO4F
REAL*8 :: SO2_eq, N_FLUX_A, N_FLUX_C
REAL*8 :: END_ALK, L5A, L5C
REAL*8 :: EQ1, EQ2, TITR_SO2
REAL*8 :: TITR_HNO3, NIT_vv, NITs_vv
REAL*8 :: NIT0, NITS0
REAL*8 :: F_SO2, FALK_A_SO2, FALK_C_SO2
REAL*8 :: EQ_BEG, F_SO2_A, F_SO2_C
REAL*8 :: ALKA, ALKC, TOTAL_ACID_FLUX
REAL*8 :: HNO3_EQ, TOT_FLUX_A, TOT_FLUX_C
REAL*8 :: FALK_A_HNO3, HNO3_vv
REAL*8 :: FALK_C_HNO3, F_HNO3_A, F_HNO3_C
REAL*8 :: EQ_1_N, EQ_2_N, F_HNO3
REAL*8 :: HNO3_SSA, HNO3_SSC, N_FLUX
REAL*8 :: HNO3_EQ_C, L6A, L6C
REAL*8 :: C_FLUX_A, C_FLUX_C, C_FLUX
REAL*8 :: HNO3_ss, HNO3_kg
REAL*8, PARAMETER :: MINDAT = 1.0d-20
REAL*8, PARAMETER :: TCVV_HNO3 = 28.97d0 / 63.0d0
!=================================================================
! SEASALT_CHEM begins here!
!=================================================================
! DTCHEM is the chemistry timestep in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Convert SO2 [v/v] to [eq/gridbox]
SO2_eq = ( 2.d0 * SO2_cd * AD(I,J,L) ) / ( TCVV(IDTSO2) * 0.064d0)
SO2_eq = MAX( SO2_eq, MINDAT )
IF ( ITS_A_FULLCHEM_SIM() ) THEN
! Convert HNO3 [v/v] to [equivalents]
HNO3_vv = STT(I,J,L,IDTHNO3)
HNO3_eq = HNO3_vv * AD(I,J,L) / ( 28.97d0 / 63d0 ) / 63.d-3
ELSE
! Get gas-phase HNO3 from ISORROPIA code
CALL GET_GNO3( I, J, L, HNO3_kg )
! Convert HNO3 [kg] first to [v/v]
HNO3_vv = HNO3_kg * ( 28.97d0 / 63d0 ) / AD(I,J,L)
! Then convert HNO3 [kg] to [equivalents]
HNO3_eq = HNO3_kg / 63d-3
ENDIF
!-----------
! SO2
!-----------
! Available flux of SO2 to accum sea salt aerosols [v/v/timestep]
L5A = EXP( -Kt1 * DTCHEM )
F_SO2_A = SO2_cd * ( 1.d0 - L5A )
F_SO2_A = MAX( F_SO2_A, 1.d-32 )
! Convert to [eq/timestep]
C_FLUX_A = 2.d0 * F_SO2_A * AD(I,J,L) / TCVV(IDTSO2) / 0.064d0
! Available flux of SO2 to coarse sea salt aerosols [v/v/timestep]
L5C = EXP( - Kt2 * DTCHEM )
F_SO2_C = SO2_cd * ( 1.d0 - L5C )
F_SO2_C = MAX( F_SO2_C, 1.0d-32 )
! Convert to [eq/timestep]
C_FLUX_C = 2.d0 * F_SO2_C * AD(I,J,L) / TCVV(IDTSO2) / 0.064d0
! Total flux of SO2 [v/v/timestep]
F_SO2 = F_SO2_A + F_SO2_C
! Total flux of SO2 [eq/timestep]
C_FLUX = C_FLUX_A + C_FLUX_C
!-----------
! HNO3
!-----------
! Available flux of HNO3 to accum sea salt aerosols [v/v/timestep]
L6A = EXP( - Kt1N * DTCHEM )
F_HNO3_A = HNO3_vv * ( 1.D0 - L6A )
F_HNO3_A = MAX( F_HNO3_A, 1.0D-32 )
! Convert to [eq/timestep]
N_FLUX_A = F_HNO3_A * AD(I,J,L)/( 28.97d0 / 63d0 )/0.063d0
! Available flux of HNO3 to coarse sea salt aerosols [v/v/timestep]
L6C = EXP( - Kt2N * DTCHEM )
F_HNO3_C = HNO3_vv * ( 1.D0 - L6C )
F_HNO3_C = MAX( F_HNO3_C, 1.0D-32 )
! convert to [eq/timestep]
N_FLUX_C = F_HNO3_C * AD(I,J,L)/( 28.97d0 / 63d0 )/0.063d0
! Total flux of HNO3
F_HNO3 = F_HNO3_A + F_HNO3_C ![v/v/timestep]
N_FLUX = N_FLUX_A + N_FLUX_C ![eq/timestep]
!-----------
! Acid
!-----------
! Total acid flux to accum sea-salt aerosols [eq/box/timestep]
TOT_FLUX_A = C_FLUX_A + N_FLUX_A
TOT_FLUX_A = MAX( TOT_FLUX_A, MINDAT )
! Total acid flux to coarse sea-salt aerosols [eq/box/timestep]
TOT_FLUX_C = C_FLUX_C + N_FLUX_C
TOT_FLUX_C = MAX( TOT_FLUX_C, MINDAT )
! Total acid flux to sea salt aerosols
TOTAL_ACID_FLUX = TOT_FLUX_A + TOT_FLUX_C
! Total available alkalinity [eq]
EQ1 = ALK1 * 0.07d0
EQ2 = ALK2 * 0.07d0
!----------------------------------------------------------------------
! NOTE: This was a sensitivity simulation, keep for future reference
! cf Alexander et al 2005 (bec, bmy, 4/13/05)
!! Total available alkalinity [eq] doubled for Sievering run
!EQ1 = ALK1 * 0.14d0
!EQ2 = ALK2 * 0.14d0
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! DIAGNOSTIC -- leave uncommented for now (bec, bmy, 4/13/05)
!! Write out beginning alkalinity [eq SO4]
!EQ_BEG = EQ1 + EQ2
!IF ( ND09 > 0 ) AD09(I,J,L,1) = AD09(I,J,L,1) + EQ_BEG
!----------------------------------------------------------------------
IF ( TOT_FLUX_A > EQ1 ) THEN
! Fraction of alkalinity available for each acid
FALK_A_SO2 = C_FLUX_A / TOT_FLUX_A
FALK_A_HNO3 = N_FLUX_A / TOT_FLUX_A
FALK_A_SO2 = MAX( FALK_A_SO2, MINDAT )
FALK_A_HNO3 = MAX( FALK_A_HNO3, MINDAT )
ELSE
FALK_A_SO2 = 1.0d0
FALK_A_HNO3 = 1.0d0
ENDIF
IF ( TOT_FLUX_C > EQ2 ) THEN
! Fraction of flkalinity available for each acid
FALK_C_SO2 = C_FLUX_C/TOT_FLUX_C
FALK_C_HNO3 = N_FLUX_C/TOT_FLUX_C
FALK_C_SO2 = MAX( FALK_C_SO2, MINDAT )
FALK_C_HNO3 = MAX( FALK_C_HNO3, MINDAT )
ELSE
FALK_C_SO2 = 1.0d0
FALK_C_HNO3 = 1.0d0
ENDIF
! Alkalinity available for S(IV) --> S(VI)
EQ_1_C = EQ1 * FALK_A_SO2
EQ_1_C = MAX( EQ_1_C, MINDAT )
EQ_1_N = EQ1 * FALK_A_HNO3
EQ_1_N = MAX( EQ_1_N, MINDAT )
EQ_2_C = EQ2 * FALK_C_SO2
EQ_2_C = MAX( EQ_2_C, MINDAT )
EQ_2_N = EQ2 * FALK_C_HNO3
EQ_2_N = MAX( EQ_2_N, MINDAT )
!-----------------
! Fine Seasalt
!-----------------
! don't produce more SO4 than available ALK or SO2
SO4E = MIN( C_FLUX_A, EQ_1_C, SO2_eq )
SO4E = MAX( SO4E, MINDAT )
! Update SO2 concentration [eq/box]
SO2_new = SO2_eq - SO4E
SO2_new = MAX( SO2_new, MINDAT )
!-----------------
! Coarse Seasalt
!-----------------
IF ( SO2_new > MINDAT ) THEN
! don't produce more SO4 than available ALK or SO2
SO4F = MIN( C_FLUX_C, SO2_new, EQ_2_C )
SO4F = MAX( SO4F, MINDAT )
!Update SO2 concentration [eq]
SO2_chem = SO2_new - SO4F
SO2_chem = MAX( SO2_chem, MINDAT )
ELSE
SO4F = MINDAT
SO2_chem = MINDAT
ENDIF
! Alkalinity titrated by S(IV) --> S(VI) [eq]
TITR_SO2 = SO4E + SO4F
!-------------------------------------------------------------------
! DIAGNOSTIC -- leave uncommented for now
!! write out in diagnostic
!IF ( ND09 > 0 ) AD09(I,J,L,2) = AD09(I,J,L,2) + TITR_SO2
!-------------------------------------------------------------------
!Modified SO2 [eq] converted back to [v/v]
SO2_ss = SO2_chem * 0.064 * TCVV(IDTSO2)/AD(I,J,L)/2.0d0
SO2_ss = MAX( SO2_ss, MINDAT )
!SO4E produced converted from [eq/timestep] to [v/v/timestep]
PSO4E = SO4E * 0.096 * TCVV(IDTSO4)/AD(I,J,L)/2.0d0
!SO4F produced converted from [eq/timestep] to [v/v/timestep]
PSO4F = SO4F * 0.096 * TCVV(IDTSO4S)/AD(I,J,L)/2.0d0
! Alkalinity titrated by HNO3
HNO3_SSA = MIN(N_FLUX_A, HNO3_EQ, EQ_1_N)
HNO3_SSA = MAX(HNO3_SSA, MINDAT)
HNO3_EQ_C = HNO3_EQ - HNO3_SSA
HNO3_EQ_C = MAX(HNO3_EQ_C, MINDAT)
HNO3_SSC = MIN(N_FLUX_C, HNO3_EQ_C, EQ_2_N)
HNO3_SSC = MAX(HNO3_SSC, MINDAT)
TITR_HNO3 = HNO3_SSA + HNO3_SSC
!----------------------------------------------------------------------
! DIAGNOSTIC -- leave commented out for now
! !write out alkalinity titrated by HNO3(g)
!IF ( ND09 > 0 ) AD09(I,J,L,3) = AD09(I,J,L,3) + TITR_HNO3
!----------------------------------------------------------------------
! HNO3 lost [eq/timestep] converted back to [v/v/timestep]
IF ( IDTHNO3 > 0 ) THEN
! Coupled sim: IDTHNO3 is defined, so use it
HNO3_ss = TITR_HNO3 * 0.063 * TCVV(IDTHNO3) / AD(I,J,L)
STT(I,J,L,IDTHNO3) = MAX( HNO3_vv - HNO3_ss, MINDAT )
ELSE
! Offline aerosol sim: IDTHNO3 isn't defined, use TCVV_HNO3
HNO3_ss = TITR_HNO3 * 0.063 * TCVV_HNO3 / AD(I,J,L)
ENDIF
! NITS produced converted from [eq/timestep] to [v/v/timestep]
PNITs(I,J,L) = HNO3_SSC * 0.063 * TCVV(IDTNITS)/AD(I,J,L)
! Modified accum alkalinity
ALKA = EQ1 - (SO4E + HNO3_SSA)
ALKA = MAX( ALKA, MINDAT )
!------------------------------------------------------------------------
! Uncomment this if you want to transport alkalinity (bec, bmy, 4/13/05)
![eq] --> [kg] --> [v/v] use this only if transporting alkalinity
!ALKAvv = (ALKA * TCVV(IDTSAL1))/(7.0d-2 * AD(I,J,L))
!ALKAvv = MAX( ALKAvv, MINDAT )
!------------------------------------------------------------------------
! Modified accum alkalinity
ALKC = EQ2 - (SO4F + HNO3_SSC)
ALKC = MAX( ALKC, MINDAT )
!------------------------------------------------------------------------
! Uncomment this if you want to transport alkalinity (bec, bmy, 4/13/05)
!! [eq] --> [kg] --> [v/v] use this only if transporting alkalinity
!ALKCvv = (ALKC * TCVV(IDTSAL2))/(7.0d-2 * AD(I,J,L))
!ALKCvv = MAX(ALKCvv, MINDAT)
!------------------------------------------------------------------------
!------------------------------------------------------------------------
! DIAGNOSTIC -- leave commented out for now (bec, bmy, 4/13/05)
!! write out ending alkalinity
!END_ALK = ALKA + ALKC
!IF ( ND09 > 0 ) AD09(I,J,L,4) = AD09(I,J,L,4) + END_ALK
!------------------------------------------------------------------------
! Return to calling program
END SUBROUTINE SEASALT_CHEM
!------------------------------------------------------------------------------
SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, 1
& O3, Hplus, KaqH2O2, KaqO3 )
!
!******************************************************************************
! Function AQCHEM_SO2 computes the reaction rates for aqueous SO2 chemistry.
! (rjp, bmy, 10/31/02, 12/12/02)
!
! Arguments as Input:
! ============================================================================
! (1 ) LWC (REAL*8) : Liquid water content [m3/m3] = 1.E-6*L [g/m3]
! (2 ) T (REAL*8) : Temperature [K]
! (3 ) P (REAL*8) : Pressure [atm]
! (4 ) SO2 (REAL*8) : SO2 mixing ratio [v/v]
! (5 ) H2O2 (REAL*8) : H2O2 mixing ratio [v/v]
! (6 ) O3 (REAL*8) : O3 mixing ratio [v/v]
! (7 ) HPLUS (REAL*8) : Concentration of H+ ion (i.e. the pH) [v/v]
!
! Arguments as Output:
! ============================================================================
! (8 ) KaqH2O2 (REAL*8) : Reaction rate for H2O2
! (9 ) KaqO3 (REAL*8) : Reaction rate for O3
!
! Chemical Reactions:
! ============================================================================
! (R1) HSO3- + H2O2(aq) + H+ => SO4-- + 2H+ + H2O [Jacob, 1986]
!
! d[S(VI)]/dt = k[H+][H2O2(aq)][HSO3-]/(1 + K[H+])
! [Seinfeld and Pandis, 1998, page 366]
!
! (R2) SO2(aq) + O3(aq) =>
! HSO3- + O3(aq) =>
! SO3-- + O3(aq) =>
! [Jacob, 1986; Jacobson, 1999]
!
! d[S(VI)]/dt = (k0[SO2(aq)] + k1[HSO3-] + K2[SO3--])[O3(aq)]
! [Seinfeld and Pandis, 1998, page 363]
!
! Reaction rates can be given as
! Ra = k [H2O2(ag)] [S(IV)] [mole/liter*s] OR
! Krate = Ra LWC R T / P [1/s]
!
! Where:
! LWC = Liquid water content(g/m3)*10-6 [m3(water)/m3(gas)]
! R = 0.08205 (atm L / mol-K), Universal gas const.
! T = Temperature (K)
! P = Pressure (atm)
!
! Procedure:
! ============================================================================
! (a ) Given [SO2] which is assumed to be total SO2 (gas+liquid) in
! equilibrium between gas and liquid phase.
!
! (b ) We can compute SO2(g) using Henry's law
! P(so2(g)) = Xg * [SO2]
! Xg = 1/(1 + Faq), Fraction of SO2 in gas
! where:
! Faq = Kheff * R * T * LWC,
! KHeff = Effective Henry's constant
!
! (c ) Then Calculate Aquous phase, S[IV] concentrations
! S[IV] = Kheff * P(so2(g) in atm) [M]
!
! (d ) The exact same procedure is applied to calculate H2O2(aq)
!
! NOTES:
! (1 ) Updated by Rokjin Park (rjp, bmy, 12/12/02)
!******************************************************************************
!
! Arguments
REAL*8, INTENT(IN) :: LWC, T, P, SO2, H2O2, O3, HPLUS
REAL*8, INTENT(OUT) :: KaqH2O2, KaqO3
! Local variables
REAL*8, PARAMETER :: R = 0.08205d0
REAL*8 :: KH2O2, RA, KS1, KS2, HCSO2
REAL*8 :: FHCSO2, XSO2G, SIV, HSO3, XSO2AQ
REAL*8 :: XHSO3, XSO3, KH1, HCH2O2, FHCH2O2
REAL*8 :: XH2O2G, H2O2aq, KO0, KO1, KO2
REAL*8 :: HCO3, XO3g, O3aq
!=================================================================
! AQCHEM_SO2 begins here!
!
! Aqueous reaction rate
! HSO3- + H2O2 + H+ => SO4-- + 2H+ + H2O [Jacob, 1986]
!=================================================================
! [Jacob, 1986]
KH2O2 = 6.31d14 * EXP( -4.76d3 / T )
! ! [Jacobson, 1999]
! KH2O2 = 7.45d07 * EXP( -15.96d0 * ( (298.15/T) - 1.) ) /
! & ( 1.d0 + 13.d0 * Hplus)
!=================================================================
! Equilibrium reaction of SO2-H2O
! SO2 + H2O = SO2(aq) (s0)
! SO2(ag) = HSO3- + H+ (s1)
! HSO3- = SO3-- + H+ (s2)
!
! Reaction constant for Aqueous chemistry -- No big difference
! between Jacob and Jacobson, choose one of them.
!
! Reaction rate dependent on Temperature is given
! H = A exp ( B (T./T - 1) )
!
! For equilibrium reactions of SO2:
! As1 Bs1 As2 Bs2
! Seinfeld 1.30d-2 7.02 6.60d-8 3.76 [1998]
! Jacob 1.30d-2 6.75 6.31d-8 5.05 [1986]
! Jacobson 1.71d-2 7.04 5.99d-8 3.74 [1996]
!=================================================================
Ks1 = 1.30d-2 * EXP( 6.75d0 * ( 298.15d0 / T - 1.d0 ) )
Ks2 = 6.31d-8 * EXP( 5.05d0 * ( 298.15d0 / T - 1.d0 ) )
! SIV Fraction
XSO2aq = 1.d0/(1.d0 + Ks1/Hplus + Ks1*Ks2/(Hplus*Hplus))
XHSO3 = 1.d0/(1.d0 + Hplus/Ks1 + Ks2/Hplus)
XSO3 = 1.d0/(1.d0 + Hplus/Ks2 + Hplus*Hplus/(Ks1*Ks2))
! Henry's constant [mol/l-atm] and Effective Henry's constant for SO2
HCSO2 = 1.22d0 * EXP( 10.55d0 * ( 298.15d0 / T - 1.d0) )
FHCSO2 = HCSO2 * (1.d0 + (Ks1/Hplus) + (Ks1*Ks2 / (Hplus*Hplus)))
XSO2g = 1.d0 / ( 1.d0 + ( FHCSO2 * R * T * LWC ) )
SIV = FHCSO2 * XSO2g * SO2 * P
! HSO3 = Ks1 * HCSO2 * XSO2g * SO2 * P
!=================================================================
! H2O2 equilibrium reaction
! H2O2 + H2O = H2O2.H2O
! H2O2.H2O = HO2- + H+ 1)
!
! Reaction rate dependent on Temperature is given
! H = A exp ( B (T./T - 1) )
!
! For equilibrium reactions of SO2
! Ah1 Bh1
! Jacob 1.58E-12 -12.49 [1986]
! Jacobson 2.20E-12 -12.52 [1996]
!=================================================================
Kh1 = 2.20d-12 * EXP( -12.52d0 * ( 298.15d0 / T - 1.d0 ) )
! Henry's constant [mol/l-atm] and Effective Henry's constant for H2O2
! [Seinfeld and Pandis, 1998]
! HCH2O2 = 7.45D4 * EXP( 24.48d0 * ( 298.15d0 / T - 1.d0) )
! [Jacobson,1999]
HCH2O2 = 7.45D4 * EXP( 22.21d0 * (298.15d0 / T - 1.d0) )
FHCH2O2 = HCH2O2 * (1.d0 + (Kh1 / Hplus))
XH2O2g = 1.d0 / ( 1.d0 + ( FHCH2O2 * R * T * LWC ) )
! H2O2aq = FHCH2O2 * XH2O2g * H2O2 * P
! Conversion rate from SO2 to SO4 via reaction with H2O2
KaqH2O2 = kh2o2 * Ks1 * FHCH2O2 * HCSO2 * XH2O2g * XSO2g
& * P * LWC * R * T ! [v/v/s]
!=================================================================
! Aqueous reactions of SO2 with O3
! SO2(aq) + O3 => (0)
! HSO3- + O3 => SO4-- + H+ + O2 (1)
! SO3-- + O3 => SO4-- + O2 (2)
!
! NOTE
! [Jacob, 1986]
! KO1 = 3.49E12 * EXP( -4.83E3 / T )
! KO2 = 7.32E14 * EXP( -4.03E3 / T )
!
! [Jacobson, 1999]
! KO0 = 2.40E+4
! KO1 = 3.70E+5 * EXP( -18.56 * ((298.15/T) - 1.))
! KO2 = 1.50E+9 * EXP( -17.72 * ((298.15/T) - 1.))
!
! Rate constants from Jacobson is larger than those of Jacob
! and results in faster conversion from S(IV) to S(VI)
! We choose Jacob 1) 2) and Jacobson 0) here
!=================================================================
KO0 = 2.40d+4
KO1 = 3.49d12 * EXP( -4.83d3 / T )
KO2 = 7.32d14 * EXP( -4.03d3 / T )
!=================================================================
! H2O2 equilibrium reaction
! O3 + H2O = O3.H2O
! HCO3 = 1.13E-2 * EXP( 8.51 * (298.15/T -1.) ), S & P
! HCO3 = 1.13E-2 * EXP( 7.72 * (298.15/T -1.) ), Jacobson
!=================================================================
! Calculate Henry's Law constant for atmospheric temperature
HCO3 = 1.13d-2 * EXP( 8.51d0 * ( 298.15d0 / T - 1.d0 ) )
XO3g = 1.d0 / ( 1.d0 + ( HCO3 * R * T * LWC ) )
! O3aq = HCO3 * XO3g * O3 * P
! Conversion rate from SO2 to SO4 via reaction with O3
KaqO3 = (KO0*XSO2AQ + KO1*XHSO3 + KO2*XSO3) * FHCSO2 * XSO2g
& * P * HCO3 * XO3g * LWC * R * T ! [v/v/s]
! Return to calling program
END SUBROUTINE AQCHEM_SO2
!------------------------------------------------------------------------------
SUBROUTINE CHEM_SO4 1,16
!
!******************************************************************************
! Subroutine CHEM_SO4 is the SO4 chemistry subroutine from Mian Chin's GOCART
! model, modified for the GEOS-CHEM model. Now also modified to account
! for production of crystalline & aqueous sulfur tracers.
! (rjp, bdf, cas, bmy, 5/31/00, 5/23/06)
!
! Module Variables Used:
! ============================================================================
! (1 ) PSO4_SO2 (REAL*8 ) : Array for P(SO4) from SO2 [v/v/timestep]
! (2 ) PSO4_ss (REAL*8 ) : Array for P(SO4) from SO2
! (coarse sea-salt aerosols) [v/v/timestep]
!
! Reaction List (by Mian Chin, chin@rondo.gsfc.nasa.gov)
! ============================================================================
! The Only production is from SO2 oxidation (save in CHEM_SO2), and the only
! loss is dry depsition here. Wet deposition will be treated in "wetdep.f".
!
! SO4 = SO4_0 * exp(-kt) + PSO4_SO2/kt * (1.-exp(-kt))
! where k = dry deposition.
!
! NOTES:
! (1 ) Now reference AD from "dao_mod.f". Added parallel DO-loops.
! Updated comments, cosmetic changes. (rjp, bdf, bmy, 9/16/02)
! (2 ) Now replace DXYP(JREF)*1d4 with routine GET_AREA_CM2 of "grid_mod.f"
! Now use function GET_TS_CHEM from "time_mod.f" (bmy, 3/27/03)
! (3 ) Now reference PBLFRAC from "drydep_mod.f". Now apply dry deposition
! to the entire PBL. (rjp, bmy, 8/1/03)
! (4 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (5 ) Now use parallel DO-loop to zero ND44_TMP (bmy, 4/14/04)
! (6 ) Now reference STT & TCVV from "tracer_mod.f" (bmy, 7/20/04)
! (7 ) Now references LCRYST from "logical_mod.f". Modified for crystalline
! and aqueous sulfate2 tracers: AS, AHS, LET, SO4aq. Also changed name
! of ND44_TMP to T44 to save space. (cas, bmy, 12/21/04)
! (8 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f" (bmy, 2/22/05)
! (9 ) Now remove reference to CMN, it's obsolete. Now reference
! ITS_IN_THE_STRAT from "tropopause_mod.f" (bmy, 8/22/05)
! (10) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
! (11) Rearrange error check to avoid SEG FAULTS (bmy, 5/23/06)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE LOGICAL_MOD, ONLY : LCRYST, LSSALT
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTSO4, IDTSO4s, IDTAS, IDTAHS
USE TRACERID_MOD, ONLY : IDTLET, IDTSO4aq, IDTNH4aq
USE TROPOPAUSE_MOD, ONLY : ITS_IN_THE_STRAT
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND44
! Local variables
INTEGER :: I, J, L, N, N_ND44
REAL*8 :: AS, AS0, AHS, AHS0, LET
REAL*8 :: LET0, SO4, SO40, SO4aq, SO4aq0
REAL*8 :: SO4s, SO40s, RKT, RKTs, E_RKT
REAL*8 :: E_RKTs, DTCHEM, AREA_CM2, FLUX
REAL*8 :: T44(IIPAR,JJPAR,LLTROP,6)
!=================================================================
! CHEM_SO4 begins here!
!=================================================================
! Return if tracers are not defined
IF ( IDTSO4 == 0 .or. IDTSO4s == 0 ) RETURN
IF ( DRYSO4 == 0 .or. DRYSO4s == 0 ) RETURN
! DTCHEM is the chemistry timestep in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Number of drydep tracers to save
N_ND44 = 2
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 4/13/05)
!%%% ! Number of drydep tracers to save
!%%% IF ( LCRYST ) THEN
!%%% N_ND44 = 6
!%%% ELSE
!%%% N_ND44 = 2
!%%% ENDIF
!------------------------------------------------------------------------------
! Zero T44 array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, N )
DO N = 1, N_ND44
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
T44(I,J,L,N) = 0d0
ENDDO
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Loop over tropospheric grid boxes
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% $OMP PARALLEL DO
!%%% $OMP+DEFAULT( SHARED )
!%%% $OMP+PRIVATE( I, J, L, AREA_CM2, RKT, RKTs )
!%%% $OMP+PRIVATE( E_RKT, E_RKTs, FLUX, SO4, SO4s, SO4aq )
!%%% $OMP+PRIVATE( AS, AHS, LET, SO40, SO40s, SO4aq0 )
!%%% $OMP+PRIVATE( AS0, AHS0, LET0 )
!%%% $OMP+SCHEDULE( DYNAMIC )
!------------------------------------------------------------------------------
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, AREA_CM2, RKT, RKTs, E_RKT )
!$OMP+PRIVATE( E_RKTs, FLUX, SO4, SO4s, SO40, SO40s )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
! Initialize for safety's sake
AREA_CM2 = 0d0
RKT = 0d0
RKTs = 0d0
E_RKT = 0d0
E_RKTs = 0d0
FLUX = 0d0
SO4 = 0d0
SO4s = 0d0
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% SO4aq = 0d0
!%%% AS = 0d0
!%%% AHS = 0d0
!%%% LET = 0d0
!------------------------------------------------------------------------------
! Skip stratospheric boxes
IF ( ITS_IN_THE_STRAT( I, J, L ) ) CYCLE
!==============================================================
! Initial concentrations before chemistry
!==============================================================
! SO4 [v/v]
SO40 = STT(I,J,L,IDTSO4)
! SO4 within coarse seasalt aerosol [v/v]
SO40s = STT(I,J,L,IDTSO4s)
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% ! SO4aq, AS, LET, AHS (if necessary) [v/v]
!%%% IF ( LCRYST ) THEN
!%%% SO4aq0 = STT(I,J,L,IDTSO4aq)
!%%% AS0 = STT(I,J,L,IDTAS)
!%%% AHS0 = STT(I,J,L,IDTAHS)
!%%% LET0 = STT(I,J,L,IDTLET)
!%%% ENDIF
!------------------------------------------------------------------------------
!==============================================================
! SO4 chemistry:
!
! (CASE 1) SO4 production from SO2 and loss by drydep
! --> see equation in header notes above
!
! (CASE 2) SO4 production from SO2 with no SO4 loss by drydep
!==============================================================
! SO4 drydep frequency [1/s]. Also accounts for the fraction
! of each vertical level that is located below the PBL top
RKT = DEPSAV(I,J,DRYSO4) * GET_FRAC_UNDER_PBLTOP( I, J, L )
! RKT > 0 denotes that we have SO4 drydep occurring
IF ( RKT > 0d0 ) THEN
!-----------------------------------------------------------
! CASE 1: SO4 production from SO2 and SO4 loss by drydep
!-----------------------------------------------------------
! Fraction of SO4 lost to drydep [unitless]
RKT = RKT * DTCHEM
! Pre-compute exponential term for use below
E_RKT = EXP( -RKT )
! Updated SO4 (gas phase) [v/v]
SO4 = ( SO40 * E_RKT ) +
& ( PSO4_SO2(I,J,L)/RKT * ( 1.d0 - E_RKT ) )
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% IF ( LCRYST ) THEN
!%%%
!%%% ! Updated SO4 (aqueous phase) [v/v]
!%%% SO4aq = ( SO4aq0 * E_RKT ) +
!%%% & ( PSO4_SO2(I,J,L)/RKT * ( 1.d0 - E_RKT ) )
!%%%
!%%% ! Updated AS, AHS, LET [v/v]
!%%% AS = AS0 * E_RKT
!%%% AHS = AHS0 * E_RKT
!%%% LET = LET0 * E_RKT
!%%%
!%%% ENDIF
!------------------------------------------------------------------------------
ELSE
!-----------------------------------------------------------
! CASE 2: Production of SO4 from SO2; no SO4 drydep loss
!-----------------------------------------------------------
!-----------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% IF ( LCRYST ) THEN
!%%%
!%%% ! SO4 production from SO2 (both gas-phase & aqueous)
!%%% SO4 = SO40 + ( 2d0 * PSO4_SO2(I,J,L) )
!%%% SO4aq = SO4aq0 + ( 2d0 * PSO4_SO2(I,J,L) )
!%%%
!%%% ! No production from AS, AHS, LET
!%%% AS = AS0
!%%% AHS = AHS0
!%%% LET = LET0
!%%%
!%%% ELSE
!%%%
!%%% ! SO4 production from SO2 [v/v/timestep]
!%%% SO4 = SO40 + PSO4_SO2(I,J,L)
!%%%
!%%% ENDIF
!-----------------------------------------------------------------------------
! SO4 production from SO2 [v/v/timestep]
SO4 = SO40 + PSO4_SO2(I,J,L)
ENDIF
!==============================================================
! SO4s (SO4 w/in seasalt aerosol) chemistry:
!
! (CASE 3) SO4s production from seasalt and loss by drydep
! --> see equation in header notes above
!
! (CASE 4) SO4s prod from seasalt w/ no SO4s loss by drydep
!==============================================================
! SO4s drydep frequency [1/s]. Also accounts for the fraction
! of each vertical level that is located below the PBL top
RKTs = DEPSAV(I,J,DRYSO4s) * GET_FRAC_UNDER_PBLTOP( I, J, L )
! RKTs > 0 indicates that SO4s drydep occurs
IF ( RKTs > 0d0 ) THEN
!-----------------------------------------------------------
! CASE 3: SO4s prod from seasalt SO4s loss by drydep
!-----------------------------------------------------------
! Fraction of SO4s lost to drydep [unitless]
RKTs = RKTs * DTCHEM
! Pre-compute exponential term for use below
E_RKTs = EXP( -RKTs )
! Updated SO4 (gas phase) [v/v]
SO4s = ( SO40s * E_RKTs ) +
& ( PSO4_ss(I,J,L)/RKTs * ( 1.d0 - E_RKTs ) )
ELSE
!--------------------------------------------------------
! CASE 4: Prod of SO4s from seasalt; no SO4s drydep loss
!--------------------------------------------------------
! SO4 production from SO2 [v/v/timestep]
SO4s = SO40s + PSO4_ss(I,J,L)
ENDIF
!==============================================================
! Final concentrations after chemistry
!==============================================================
! Error check
IF ( SO4 < SMALLNUM ) SO4 = 0d0
IF ( SO4s < SMALLNUM ) SO4s = 0d0
! Final concentrations [v/v]
STT(I,J,L,IDTSO4) = SO4
STT(I,J,L,IDTSO4s) = SO4s
!-----------------------------------------------------------------------------
!%%% Currently under development (bmy, 3/15/05)
!%%% ! SO4aq, AS, AHS, LET (if necessary)
!%%% IF ( LCRYST ) THEN
!%%%
!%%% ! Error check
!%%% IF ( SO4aq < SMALLNUM ) SO4aq = 0d0
!%%% IF ( AS < SMALLNUM ) AS = 0d0
!%%% IF ( AHS < SMALLNUM ) AHS = 0d0
!%%% IF ( LET < SMALLNUM ) LET = 0d0
!%%%
!%%% ! Final SO4aq, AS, AHS, LET [v/v]
!%%% STT(I,J,L,IDTSO4aq) = SO4aq
!%%% STT(I,J,L,IDTAS) = AS
!%%% STT(I,J,L,IDTAHS) = AHS
!%%% STT(I,J,L,IDTLET) = LET
!%%%
!%%% ENDIF
!-----------------------------------------------------------------------------
!==============================================================
! ND44 Diagnostic: Drydep flux of SO4 and the crystalline &
! aqueous tracers (AS, AHS, LET, SO4aq) in [molec/cm2/s]
!==============================================================
IF ( ND44 > 0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! SO4 drydep flux [molec/cm2/s]
FLUX = SO40 - SO4 + PSO4_SO2(I,J,L)
FLUX = FLUX * AD(I,J,L) / TCVV(IDTSO4)
FLUX = FLUX * XNUMOL(IDTSO4) / AREA_CM2 / DTCHEM
T44(I,J,L,1) = T44(I,J,L,1) + FLUX
! SO4s drydep flux [molec/cm2/s]
FLUX = SO40s - SO4s + PSO4_ss(I,J,L)
FLUX = FLUX * AD(I,J,L) / TCVV(IDTSO4s)
FLUX = FLUX * XNUMOL(IDTSO4s) / AREA_CM2 / DTCHEM
T44(I,J,L,2) = T44(I,J,L,2) + FLUX
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% ! SO4aq, AS, AHS, LET drydep fluxes (if necessary)
!%%% IF ( LCRYST ) THEN
!%%%
!%%% ! SO4aq drydep flux [molec/cm2/s]
!%%% FLUX = SO4aq0 - SO4aq + PSO4_SO2(I,J,L)
!%%% FLUX = FLUX * AD(I,J,L) / TCVV(IDTSO4aq)
!%%% FLUX = FLUX * XNUMOL(IDTSO4aq) / AREA_CM2 / DTCHEM
!%%% T44(I,J,L,3) = T44(I,J,L,3) + FLUX
!%%%
!%%% ! AS drydep flux [molec/cm2/s]
!%%% FLUX = AS0 - AS
!%%% FLUX = FLUX * AD(I,J,L) / TCVV(IDTAS)
!%%% FLUX = FLUX * XNUMOL(IDTAS) / AREA_CM2 / DTCHEM
!%%% T44(I,J,L,4) = T44(I,J,L,4) + FLUX
!%%%
!%%% ! AHS drydep flux [molec/cm2/s]
!%%% FLUX = AHS0 - AHS
!%%% FLUX = FLUX * AD(I,J,L) / TCVV(IDTAHS)
!%%% FLUX = FLUX * XNUMOL(IDTAHS) / AREA_CM2 / DTCHEM
!%%% T44(I,J,L,5) = T44(I,J,L,5) + FLUX
!%%%
!%%% ! LET drydep flux [molec/cm2/s]
!%%% FLUX = LET0 - LET
!%%% FLUX = FLUX * AD(I,J,L) / TCVV(IDTLET)
!%%% FLUX = FLUX * XNUMOL(IDTLET) / AREA_CM2 / DTCHEM
!%%% T44(I,J,L,6) = T44(I,J,L,6) + FLUX
!%%%
!%%% ENDIF
!------------------------------------------------------------------------------
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, LLTROP
! Sum SO4, SO4s drydep fluxes in the vertical [molec/cm2/s]
AD44(I,J,DRYSO4, 1) = AD44(I,J,DRYSO4, 1) + T44(I,J,L,1)
AD44(I,J,DRYSO4s,1) = AD44(I,J,DRYSO4s,1) + T44(I,J,L,2)
!------------------------------------------------------------------------------
!%%% Currently under development (rjp, bmy, 3/15/05)
!%%% ! Sum SO4aq, AS, AHS, LET drydep fluxes (if necessary)
!%%% IF ( LCRYST ) THEN
!%%% AD44(I,J,DRYSO4aq,1) = AD44(I,J,DRYSO4aq,1)+T44(I,J,L,3)
!%%% AD44(I,J,DRYAS, 1) = AD44(I,J,DRYAS, 1)+T44(I,J,L,4)
!%%% AD44(I,J,DRYAHS, 1) = AD44(I,J,DRYAHS, 1)+T44(I,J,L,5)
!%%% AD44(I,J,DRYLET, 1) = AD44(I,J,DRYLET, 1)+T44(I,J,L,6)
!%%% ENDIF
!------------------------------------------------------------------------------
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_SO4
!------------------------------------------------------------------------------
!SUBROUTINE PHASE_SO4
!
! *** Currently under development ***
!
!END SUBROUTINE PHASE_SO4
!------------------------------------------------------------------------------
!SUBROUTINE PHASE_RADIATIVE
!
! *** Currently under development ***
!
!END SUBROUTINE PHASE_RADIATIVE
!------------------------------------------------------------------------------
SUBROUTINE CHEM_MSA 1,12
!
!******************************************************************************
! Subroutine CHEM_MSA is the SO4 chemistry subroutine from Mian Chin's GOCART
! model, modified for the GEOS-CHEM model. (rjp, bdf, bmy, 5/31/00, 10/25/05)
!
! Module Variables Used:
! ============================================================================
! (1 ) PMSA_DMS (REAL*8 ) : Array for P(MSA) from DMS [v/v/timestep]
!
! Reaction List (by Mian Chin, chin@rondo.gsfc.nasa.gov)
! ============================================================================
! The Only production is from DMS oxidation (saved in CHEM_DMS), and the only
! loss is dry depsition here. Wet deposition will be treated in "wetdep.f".
!
! MSA = MSA_0 * exp(-dt) + PMSA_DMS/kt * (1.-exp(-kt))
! where k = dry deposition.
!
! NOTES:
! (1 ) Now reference AD from "dao_mod.f". Added parallel DO-loops.
! Updated comments, cosmetic changes. (rjp, bmy, bdf, 9/16/02)
! (2 ) Now replace DXYP(JREF)*1d4 with routine GET_AREA_CM2 of "grid_mod.f"
! Now use function GET_TS_CHEM from "time_mod.f" (bmy, 3/27/03)
! (3 ) Now reference PBLFRAC from "drydep_mod.f". Now apply dry deposition
! to the entire PBL. (rjp, bmy, 8/1/03)
! (4 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (5 ) Now use parallel DO-loop to zero ND44_TMP (bmy, 4/14/04)
! (6 ) Now references STT & TCVV from "tracer_mod.f" (bmy, 7/20/04)
! (7 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f". Also reference GET_PBL_MAX_L from "pbl_mix_mod.f"
! Vertical DO-loops can run up to PBL_MAX and not LLTROP. Also
! remove reference to header file CMN. (bmy, 2/22/05)
! (8 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP, GET_PBL_MAX_L
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTMSA
# include "CMN_SIZE" ! Size parameters
# include "CMN_GCTM" ! AIRMW
# include "CMN_DIAG" ! ND44
! Local variables
INTEGER :: I, J, L, PBL_MAX
REAL*8 :: DTCHEM, MSA0, MSA, RK
REAL*8 :: RKT, FLUX, AREA_CM2, F_UNDER_TOP
REAL*8 :: ND44_TMP(IIPAR,JJPAR,LLTROP)
!=================================================================
! CHEM_MSA begins here!
!=================================================================
IF ( IDTMSA == 0 .or. DRYMSA == 0 ) RETURN
! DTCHEM is the chemistry interval in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Model level where maximum PBL height occurs
PBL_MAX = GET_PBL_MAX_L()
! Zero ND44_TMP array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
ND44_TMP(I,J,L) = 0d0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Loop over tropospheric grid boxes
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, F_UNDER_TOP, MSA0, RKT, MSA, AREA_CM2, FLUX )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, PBL_MAX
DO J = 1, JJPAR
DO I = 1, IIPAR
! Fraction of box (I,J,L) underneath the PBL top [unitless]
F_UNDER_TOP = GET_FRAC_UNDER_PBLTOP( I, J, L )
! Only apply drydep loss to boxes w/in the PBL
IF ( F_UNDER_TOP > 0 ) THEN
! Initial MSA [v/v]
MSA0 = STT(I,J,L,IDTMSA)
! MSA drydep frequency [1/s]. Also accounts for the fraction
! of each grid box (I,J,L) that is located beneath the PBL top
RKT = DEPSAV(I,J,DRYMSA) * F_UNDER_TOP
! RKT > 0 denotes that we have drydep occurring
IF ( RKT > 0.d0 ) THEN
! Fraction of MSA lost to drydep [unitless]
RKT = RKT * DTCHEM
! Modified MSA concentration
MSA = ( MSA0 * EXP( -RKT ) ) +
& ( PMSA_DMS(I,J,L)/RKT * ( 1d0 - EXP( -RKT ) ) )
ELSE
! MSA production from DMS [v/v/timestep]
MSA = MSA0 + PMSA_DMS(I,J,L)
ENDIF
! Final MSA [v/v]
IF ( MSA < SMALLNUM ) MSA = 0d0
STT(I,J,L,IDTMSA) = MSA
!===========================================================
! ND44 Diagnostic: Drydep flux of MSA [molec/cm2/s]
!===========================================================
IF ( ND44 > 0 .and. RKT > 0d0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Convert [v/v/timestep] to [molec/cm2/s]
FLUX = MSA0 - MSA + PMSA_DMS(I,J,L)
FLUX = FLUX * AD(I,J,L) / TCVV(IDTMSA)
FLUX = FLUX * XNUMOL(IDTMSA) / AREA_CM2 / DTCHEM
! Store dryd flux in ND44_TMP as a placeholder
ND44_TMP(I,J,L) = ND44_TMP(I,J,L) + FLUX
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, PBL_MAX
AD44(I,J,DRYMSA,1) = AD44(I,J,DRYMSA,1) + ND44_TMP(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_MSA
!------------------------------------------------------------------------------
SUBROUTINE CHEM_NH3 1,12
!
!******************************************************************************
! Subroutine CHEM_NH3 removes NH3 from the surface via dry deposition.
! (rjp, bdf, bmy, 1/2/02, 10/25/05)
!
! Reaction List:
! ============================================================================
! (1 ) NH3 = NH3_0 * EXP( -dt ) where d = dry deposition rate [s-1]
!
! NOTES:
! (1 ) Now reference AD from "dao_mod.f". Added parallel DO-loops.
! Updated comments, cosmetic changes. (rjp, bmy, bdf, 9/16/02)
! (2 ) Now replace DXYP(J+J0)*1d4 with routine GET_AREA_CM2 from "grid_mod.f"
! Now use function GET_TS_CHEM from "time_mod.f" (bmy, 3/27/03)
! (3 ) Now reference PBLFRAC from "drydep_mod.f". Now apply dry deposition
! to the entire PBL. Added L and FREQ variables. Recode to avoid
! underflow from the EXP() function. (rjp, bmy, 8/1/03)
! (4 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (5 ) Now use parallel DO-loop to zero ND44_TMP (bmy, 4/14/04)
! (6 ) Now references STT & TCVV from "tracer_mod.f" Also remove reference to
! CMN, it's not needed(bmy, 7/20/04)
! (7 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f". Also reference GET_PBL_MAX_L from "pbl_mix_mod.f"
! Vertical DO-loops can run up to PBL_MAX and not LLTROP. (bmy, 2/22/05)
! (8 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP, GET_PBL_MAX_L
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTNH3
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND44
! Local variables
INTEGER :: I, J, L, PBL_MAX
REAL*8 :: DTCHEM, NH30, NH3
REAL*8 :: FREQ, AREA_CM2, FLUX, F_UNDER_TOP
REAL*8 :: ND44_TMP(IIPAR,JJPAR,LLTROP)
!=================================================================
! CHEM_NH3 begins here!
!=================================================================
IF ( IDTNH3 == 0 .or. DRYNH3 == 0 ) RETURN
! DTCHEM is the chemistry interval in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Model level where maximum PBL height occurs
PBL_MAX = GET_PBL_MAX_L()
! Zero ND44_TMP array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
ND44_TMP(I,J,L) = 0d0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, F_UNDER_TOP, FREQ, NH30, NH3, AREA_CM2, FLUX )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, PBL_MAX
DO J = 1, JJPAR
DO I = 1, IIPAR
! Fraction of box (I,J,L) underneath the PBL top [unitless]
F_UNDER_TOP = GET_FRAC_UNDER_PBLTOP( I, J, L )
! Only apply drydep to boxes w/in the PBL
IF ( F_UNDER_TOP > 0d0 ) THEN
! NH3 drydep frequency [1/s]. Also accounts for the fraction
! of each grid box (I,J,L) that is located beneath the PBL top
FREQ = DEPSAV(I,J,DRYNH3) * F_UNDER_TOP
! Only compute drydep loss if FREQ is nonzero
IF ( FREQ > 0d0 ) THEN
! Initial NH3 [v/v]
NH30 = STT(I,J,L,IDTNH3)
! Amount of NH3 lost to drydep [v/v]
NH3 = NH30 * ( 1d0 - EXP( -FREQ * DTCHEM ) )
! Prevent underflow condition
IF ( NH3 < SMALLNUM ) NH3 = 0d0
! Subtract NH3 lost to drydep from initial NH3 [v/v]
STT(I,J,L,IDTNH3) = NH30 - NH3
!========================================================
! ND44 diagnostic: Drydep flux of NH3 [molec/cm2/s]
!========================================================
IF ( ND44 > 0 .and. NH3 > 0d0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Convert drydep loss from [v/v/timestep] to [molec/cm2/s]
FLUX = NH3 * AD(I,J,L) / TCVV(IDTNH3)
FLUX = FLUX * XNUMOL(IDTNH3) / AREA_CM2 / DTCHEM
! Store dryd flx in ND44_TMP as a placeholder
ND44_TMP(I,J,L) = ND44_TMP(I,J,L) + FLUX
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, PBL_MAX
AD44(I,J,DRYNH3,1) = AD44(I,J,DRYNH3,1) + ND44_TMP(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_NH3
!------------------------------------------------------------------------------
SUBROUTINE CHEM_NH4 1,12
!
!******************************************************************************
! Subroutine CHEM_NH4 removes NH4 from the surface via dry deposition.
! (rjp, bdf, bmy, 1/2/02, 10/25/05)
!
! Reaction List:
! ============================================================================
! (1 ) NH4 = NH4_0 * EXP( -dt ) where d = dry deposition rate [s-1]
!
! NOTES:
! (1 ) Now reference AD from "dao_mod.f". Added parallel DO-loops.
! Updated comments, cosmetic changes. (rjp, bmy, bdf, 9/16/02)
! (2 ) Now replace DXYP(JREF)*1d4 with routine GET_AREA_CM2 of "grid_mod.f".
! Now use function GET_TS_CHEM from "time_mod.f" (bmy, 3/27/03)
! (3 ) Now reference PBLFRAC from "drydep_mod.f". Now apply dry deposition
! to the entire PBL. Added L and FREQ variables. Recode to avoid
! underflow from EXP(). (rjp, bmy, 8/1/03)
! (4 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (5 ) Now use parallel DO-loop to zero ND44_TMP (bmy, 4/14/04)
! (6 ) Now reference STT & TCVV from "tracer_mod.f". Also remove reference
! to CMN, it's not needed (bmy, 7/20/04)
! (7 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f". Also reference GET_PBL_MAX_L from "pbl_mix_mod.f"
! Vertical DO-loops can run up to PBL_MAX and not LLTROP. (bmy, 2/22/05)
! (8 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP, GET_PBL_MAX_L
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTNH4
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND44
! Local variables
INTEGER :: I, J, L, PBL_MAX
REAL*8 :: DTCHEM, NH4, NH40
REAL*8 :: FREQ, FLUX, AREA_CM2, F_UNDER_TOP
REAL*8 :: ND44_TMP(IIPAR,JJPAR,LLTROP)
!=================================================================
! CHEM_NH4 begins here!
!=================================================================
IF ( IDTNH4 == 0 .or. DRYNH4 == 0 ) RETURN
! DTCHEM is the chemistry interval in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Model level where maximum PBL height occurs
PBL_MAX = GET_PBL_MAX_L()
! Zero ND44_TMP array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
ND44_TMP(I,J,L) = 0d0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, F_UNDER_TOP, FREQ, NH40, NH4, AREA_CM2, FLUX )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, PBL_MAX
DO J = 1, JJPAR
DO I = 1, IIPAR
! Fraction of box (I,J,L) underneath the PBL top [unitless]
F_UNDER_TOP = GET_FRAC_UNDER_PBLTOP( I, J, L )
! Only apply drydep to boxes w/in the PBL
IF ( F_UNDER_TOP > 0d0 ) THEN
! NH4 drydep frequency [1/s]. Also accounts for the fraction
! of each grid box (I,J,L) that is located beneath the PBL top
FREQ = DEPSAV(I,J,DRYNH4) * F_UNDER_TOP
! Only apply drydep loss if FREQ is nonzero
IF ( FREQ > 0d0 ) THEN
! Initial NH4 [v/v]
NH40 = STT(I,J,L,IDTNH4)
! Amount of NH4 lost to drydep [v/v]
NH4 = NH40 * ( 1d0 - EXP( -FREQ * DTCHEM ) )
! Prevent underflow condition
IF ( NH4 < SMALLNUM ) NH4 = 0d0
! Subtract NH4 lost to drydep from initial NH4 [v/v]
STT(I,J,L,IDTNH4) = NH40 - NH4
!========================================================
! ND44 diagnostic: Drydep flux of NH4 [molec/cm2/s]
!========================================================
IF ( ND44 > 0 .and. NH4 > 0d0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Convert drydep loss from [v/v/timestep] to [molec/cm2/s]
FLUX = NH4 * AD(I,J,L) / TCVV(IDTNH4)
FLUX = FLUX * XNUMOL(IDTNH4) / AREA_CM2 / DTCHEM
! Store dryd flx in ND44_TMP as a placeholder
ND44_TMP(I,J,L) = ND44_TMP(I,J,L) + FLUX
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, PBL_MAX
AD44(I,J,DRYNH4,1) = AD44(I,J,DRYNH4,1) + ND44_TMP(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_NH4
!------------------------------------------------------------------------------
SUBROUTINE CHEM_NH4aq,12
!
!******************************************************************************
! Subroutine CHEM_NH4aq removes NH4aq from the surface via dry deposition.
! (cas, bmy, 1/6/05, 10/25/05)
!
! Reaction List:
! ============================================================================
! (1 ) NH4aq = NH4_0aq * EXP( -dt ) where d = dry deposition rate [s-1]
!
! NOTES:
! (1 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f". Also reference GET_PBL_MAX_L from "pbl_mix_mod.f"
! Vertical DO-loops can run up to PBL_MAX and not LLTROP. (bmy, 2/22/05)
! (31) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP, GET_PBL_MAX_L
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTNH4aq
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND44
! Local variables
INTEGER :: I, J, L, PBL_MAX
REAL*8 :: DTCHEM, NH4aq, NH4aq0
REAL*8 :: FREQ, FLUX, AREA_CM2, F_UNDER_TOP
REAL*8 :: T44(IIPAR,JJPAR,LLTROP)
!=================================================================
! CHEM_NH4 begins here!
!=================================================================
IF ( IDTNH4aq == 0 .or. DRYNH4aq == 0 ) RETURN
! DTCHEM is the chemistry interval in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Model level where maximum PBL height occurs
PBL_MAX = GET_PBL_MAX_L()
! Zero T44 array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
T44(I,J,L) = 0d0
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, F_UNDER_TOP, FREQ )
!$OMP+PRIVATE( NH4aq0, NH4aq, AREA_CM2, FLUX )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, PBL_MAX
DO J = 1, JJPAR
DO I = 1, IIPAR
! Fraction of box (I,J,L) underneath the PBL top [unitless]
F_UNDER_TOP = GET_FRAC_UNDER_PBLTOP( I, J, L )
! Only apply drydep to boxes w/in the PBL
IF ( F_UNDER_TOP > 0d0 ) THEN
! NH4 drydep frequency [1/s] -- PBLFRAC accounts for the fraction
! of each grid box (I,J,L) that is located beneath the PBL top
FREQ = DEPSAV(I,J,DRYNH4aq) * F_UNDER_TOP
! Only apply drydep loss if FREQ is nonzero
IF ( FREQ > 0d0 ) THEN
! Initial NH4 [v/v]
NH4aq0 = STT(I,J,L,IDTNH4aq)
! Amount of NH4 lost to drydep [v/v]
NH4aq = NH4aq0 * ( 1d0 - EXP( -FREQ * DTCHEM ) )
! Prevent underflow condition
IF ( NH4aq < SMALLNUM ) NH4aq = 0d0
! Subtract NH4 lost to drydep from initial NH4 [v/v]
STT(I,J,L,IDTNH4aq) = NH4aq0 - NH4aq
!========================================================
! ND44 diagnostic: Drydep flux of NH4 [molec/cm2/s]
!========================================================
IF ( ND44 > 0 .and. NH4aq > 0d0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Convert drydep loss from [v/v/timestep] to [molec/cm2/s]
FLUX = NH4aq * AD(I,J,L) / TCVV(IDTNH4aq)
FLUX = FLUX * XNUMOL(IDTNH4aq) / AREA_CM2 / DTCHEM
! Store dryd flx in ND44_TMP as a placeholder
T44(I,J,L) = T44(I,J,L) + FLUX
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, PBL_MAX
AD44(I,J,DRYNH4aq,1) = AD44(I,J,DRYNH4aq,1) + T44(I,J,L)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_NH4aq
!------------------------------------------------------------------------------
SUBROUTINE CHEM_NIT 1,13
!
!******************************************************************************
! Subroutine CHEM_NIT removes SULFUR NITRATES (NIT) from the surface
! via dry deposition. (rjp, bdf, bmy, 1/2/02, 5/23/06)
!
! Reaction List:
! ============================================================================
! (1 ) NIT = NIT_0 * EXP( -dt ) where d = dry deposition rate [s-1]
!
! NOTES:
! (1 ) Now reference AD from "dao_mod.f". Added parallel DO-loops.
! Updated comments, cosmetic changes. (rjp, bmy, bdf, 9/20/02)
! (2 ) Now replace DXYP(J+J0)*1d4 with routine GET_AREA_CM2 from "grid_mod.f".
! Now use function GET_TS_CHEM from "time_mod.f" (bmy, 3/27/03)
! (3 ) Now reference PBLFRAC from "drydep_mod.f". Now apply dry deposition
! to the entire PBL. Added L and FREQ variables. Recode to avoid
! underflow from EXP(). (rjp, bmy, 8/1/03)
! (4 ) Now use ND44_TMP array to store vertical levels of drydep flux, then
! sum into AD44 array. This preents numerical differences when using
! multiple processors. (bmy, 3/24/04)
! (5 ) Now use parallel DO-loop to zero ND44_TMP (bmy, 4/14/04)
! (6 ) Now reference STT & TCVV from "tracer_mod.f". Also remove reference
! to CMN, it's not needed anymore. (bmy, 7/20/04)
! (7 ) Replace PBLFRAC from "drydep_mod.f" with GET_FRAC_UNDER_PBLTOP from
! "pbl_mix_mod.f". Also reference GET_PBL_MAX_L from "pbl_mix_mod.f"
! Vertical DO-loops can run up to PBL_MAX and not LLTROP. (bmy, 2/22/05)
! (8 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
! (9 ) Rearrange error check to avoid SEG FAULTS (bmy, 5/23/06)
!******************************************************************************
!
! References to F90 modules
USE DAO_MOD, ONLY : AD
USE DIAG_MOD, ONLY : AD44
USE DRYDEP_MOD, ONLY : DEPSAV
USE GRID_MOD, ONLY : GET_AREA_CM2
USE LOGICAL_MOD, ONLY : LSSALT
USE PBL_MIX_MOD, ONLY : GET_FRAC_UNDER_PBLTOP, GET_PBL_MAX_L
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : STT, TCVV, XNUMOL
USE TRACERID_MOD, ONLY : IDTNIT, IDTNITs
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND44
! Local variables
INTEGER :: I, J, L, N, N_ND44, PBL_MAX
REAL*8 :: DTCHEM, NIT, NITs, NIT0, NIT0s, E_RKT
REAL*8 :: E_RKTs, FLUX, FREQ, FREQs, RKT, RKTs
REAL*8 :: AREA_CM2, F_UNDER_TOP
REAL*8 :: T44(IIPAR,JJPAR,LLTROP,2)
!=================================================================
! CHEM_NIT begins here!
!=================================================================
! Return if tracers are not defined
IF ( IDTNIT == 0 .or. IDTNITs == 0 ) RETURN
IF ( DRYNIT == 0 .or. DRYNITs == 0 ) RETURN
! DTCHEM is the chemistry interval in seconds
DTCHEM = GET_TS_CHEM() * 60d0
! Model level where maximum PBL height occurs
PBL_MAX = GET_PBL_MAX_L()
! Number of tracers for ND44
N_ND44 = 2
! Zero ND44 array
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, N )
DO N = 1, N_ND44
DO L = 1, LLTROP
DO J = 1, JJPAR
DO I = 1, IIPAR
T44(I,J,L,N) = 0d0
ENDDO
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, NIT0, NIT0s, NIT )
!$OMP+PRIVATE( NITs, FREQ, FREQs, F_UNDER_TOP, RKT, E_RKT )
!$OMP+PRIVATE( RKTs, E_RKTs, AREA_CM2, FLUX )
!$OMP+SCHEDULE( DYNAMIC )
DO L = 1, PBL_MAX
DO J = 1, JJPAR
DO I = 1, IIPAR
! Initial NITRATE [v/v]
NIT0 = STT(I,J,L,IDTNIT)
! Initial NITRATE w/in seasalt [v/v]
NIT0s = STT(I,J,L,IDTNITs)
! Initialize variables
NIT = 0d0
NITs = 0d0
FREQ = 0d0
FREQs = 0d0
! Fraction of box (I,J,L) underneath the PBL top [unitless]
F_UNDER_TOP = GET_FRAC_UNDER_PBLTOP( I, J, L )
! Only apply drydep to boxes w/in the PBL
IF ( F_UNDER_TOP > 0d0 ) THEN
!===========================================================
! NIT chemistry
!===========================================================
! NIT drydep frequency [1/s]. Also accounts for the fraction
! of each vertical level that is located below the PBL top
FREQ = DEPSAV(I,J,DRYNIT) * F_UNDER_TOP
! If there is drydep ...
IF ( FREQ > 0d0 ) THEN
! Fraction of NIT lost to drydep [unitless] (bec, 12/15/04)
RKT = FREQ * DTCHEM
! Pre-compute the exponential term
E_RKT = EXP( -RKT )
! Amount of NITRATE lost to drydep [v/v]
NIT = NIT0 * ( 1d0 - E_RKT )
! Prevent underflow condition
IF ( NIT < SMALLNUM ) NIT = 0d0
! Subtract NITRATE lost to drydep from initial NITRATE [v/v]
STT(I,J,L,IDTNIT) = NIT0 - NIT
ELSE
! No deposition occurs
NIT = 0d0
ENDIF
!===========================================================
! NITs chemistry
!===========================================================
! NITs drydep frequency [1/s]. Also accounts for the fraction
! of each vertical level that is located below the PBL top
FREQs = DEPSAV(I,J,DRYNITs) * F_UNDER_TOP
! If there is drydep ...
IF ( FREQs > 0d0 ) THEN
! Fraction of NIT lost to drydep [unitless] (bec, 12/15/04)
RKTs = FREQs * DTCHEM
! Pre-compute the exponential term
E_RKTs = EXP( -RKTs )
! Compute new NIT concentration [v/v],
! updated for seasalt chemistry
NITs = ( NIT0s * E_RKTs ) +
& ( PNITs(I,J,L)/RKTs * ( 1.d0 - E_RKTs ) )
ELSE
! NIT prod from HNO3 uptake on fine sea-salt [v/v/timestep]
NITs = NIT0s + PNITs(I,J,L)
ENDIF
! Store final concentration in STT [v/v]
STT(I,J,L,IDTNITs) = NITs
!========================================================
! ND44 Diagnostic: Drydep flux of NIT [molec/cm2/s]
!========================================================
IF ( ND44 > 0 ) THEN
! Surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
!-------------
! NIT drydep
!-------------
! If NIT drydep occurs ...
IF ( FREQ > 0d0 ) THEN
! Convert from [v/v/timestep] to [molec/cm2/s]
FLUX = NIT * AD(I,J,L) / TCVV(IDTNIT)
FLUX = FLUX * XNUMOL(IDTNIT) / AREA_CM2 / DTCHEM
! Store dryd flx in ND44_TMP as a placeholder
T44(I,J,L,1) = T44(I,J,L,1) + FLUX
ENDIF
!-------------
! NITs drydep
!-------------
! NOTE: if drydep doesn't occur then we still have
! production from seasalt (bec, bmy, 4/13/05)
! Convert from [v/v/timestep] to [molec/cm2/s]
FLUX = NIT0s - NITs + PNITs(I,J,L)
FLUX = FLUX * AD(I,J,L) / TCVV(IDTNITs)
FLUX = FLUX * XNUMOL(IDTNITs) / AREA_CM2 / DTCHEM
! Store dryd flx in ND44_TMP as a placeholder
T44(I,J,L,2) = T44(I,J,L,2) + FLUX
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!===============================================================
! ND44: Sum drydep fluxes by level into the AD44 array in
! order to ensure that we get the same results w/ sp or mp
!===============================================================
IF ( ND44 > 0 ) THEN
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L )
DO J = 1, JJPAR
DO I = 1, IIPAR
DO L = 1, PBL_MAX
AD44(I,J,DRYNIT, 1) = AD44(I,J,DRYNIT, 1) + T44(I,J,L,1)
AD44(I,J,DRYNITs,1) = AD44(I,J,DRYNITs,1) + T44(I,J,L,2)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE CHEM_NIT
!-----------------------------------------------------------------------------
SUBROUTINE EMISSSULFATE 2,29
!
!******************************************************************************
! Subroutine EMISSSULFATE is the interface between the GEOS-CHEM model and
! the sulfate emissions routines in "sulfate_mod.f" (bmy, 6/7/00, 10/3/05)
!
! NOTES:
! (1 ) BXHEIGHT is now dimensioned IIPAR,JJPAR,LLPAR (bmy, 9/26/01)
! (2 ) Removed obsolete commented out code from 9/01 (bmy, 10/24/01)
! (3 ) Now reference all arguments except FIRSTEMISS, LENV, LEEV from
! header files or F90 modules. Removed NSRCE, MONTH, JDAY,
! LWI, BXHEIGHT, DXYP, AD, PTOP, SIGE, PS, PBL, XTRA2, STT, DATA_DIR,
! JYEAR from the arg list. Now reference GET_PEDGE from F90 module
! "pressure_mod.f" to compute grid box edge pressures. Now uses
! GET_SEASON from "time_mod.f" to get the season. Now references
! IDTDMS, IDTSO2, etc from "tracerid_mod.f". Now make FIRSTEMISS
! a local SAVEd variable. Now call READ_BIOMASS_NH3 to read NH3
! biomass and biofuel emissions. (bmy, 12/13/02)
! (4 ) Now call READ_NATURAL_NH3 to read the NH3 source from natural
! emissions. (rjp, bmy, 3/23/03)
! (5 ) Now use functions GET_SEASON and GET_MONTH from the new "time_mod.f"
! (bmy, 3/27/03)
! (6 ) Added first-time printout message (bmy, 4/6/04)
! (7 ) Now references CMN_SETUP. Now read ship SO2 if LSHIPSO2=T. Also
! references ITS_A_NEW_MONTH from "time_mod.f". (bec, bmy, 5/20/04)
! (8 ) Now references STT and ITS_AN_AEROSOL_SIM from "tracer_mod.f".
! Now references LSHIPSO2 from "logical_mod.f" (bmy, 7/20/04)
! (9 ) Now references GET_YEAR from "time_mod.f". (bmy, 8/1/05)
! (10) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
!******************************************************************************
!
! References to F90 modules
USE LOGICAL_MOD, ONLY : LSHIPSO2
USE TIME_MOD, ONLY : GET_SEASON, GET_MONTH
USE TIME_MOD, ONLY : GET_YEAR, ITS_A_NEW_MONTH
USE TRACER_MOD, ONLY : STT, ITS_AN_AEROSOL_SIM
USE TRACERID_MOD, ONLY : IDTNITs, IDTSO4s
USE TRACERID_MOD, ONLY : IDTDMS, IDTSO2
USE TRACERID_MOD, ONLY : IDTSO4, IDTNH3
# include "CMN_SIZE" ! Size parameters
! Local variables
LOGICAL, SAVE :: FIRSTEMISS = .TRUE.
INTEGER :: NSEASON, MONTH, YEAR
!=================================================================
! EMISSSULFATE begins here!
!=================================================================
! Do only on the first timestep
IF ( FIRSTEMISS ) THEN
! Echo info
WRITE( 6, '(a)' ) REPEAT( '=', 79 )
WRITE( 6, 100 )
WRITE( 6, 110 )
WRITE( 6, 120 )
WRITE( 6, 130 )
WRITE( 6, '(a)' ) REPEAT( '=', 79 )
! FORMAT strings
100 FORMAT( 'S U L F A T E A E R O S O L E M I S S I O N S' )
110 FORMAT( 'Routines originally by Mian Chin''s GOCART model' )
120 FORMAT( 'Modified for GEOS-CHEM by R. Park and R. Yantosca' )
130 FORMAT( 'Last Modification Date: 4/6/04' )
! Initialize arrays
CALL INIT_SULFATE
! Read emissions from volcanoes
IF ( LENV ) CALL READ_NONERUP_VOLC
IF ( LEEV ) CALL READ_ERUP_VOLC
! We have now gone thru the first timestep
FIRSTEMISS = .FALSE.
ENDIF
! Get the season and month
NSEASON = GET_SEASON()
MONTH = GET_MONTH()
YEAR = GET_YEAR()
!=================================================================
! If this is a new month, read in the monthly mean quantities
!=================================================================
IF ( ITS_A_NEW_MONTH() ) THEN
! Read monthly mean data
CALL READ_SST( MONTH, YEAR )
CALL READ_OCEAN_DMS( MONTH )
CALL READ_AIRCRAFT_SO2( MONTH )
CALL READ_BIOMASS_SO2( MONTH )
CALL READ_ANTHRO_SOx( MONTH, NSEASON )
CALL READ_ANTHRO_NH3( MONTH )
CALL READ_BIOMASS_NH3( MONTH )
CALL READ_NATURAL_NH3( MONTH )
! Also read ship exhaust SO2 if necessary
IF ( LSHIPSO2 ) CALL READ_SHIP_SO2( MONTH )
! Read oxidants for the offline simulation only
IF ( ITS_AN_AEROSOL_SIM() ) CALL READ_OXIDANT( MONTH )
ENDIF
!=================================================================
! Add emissions into the STT tracer array
!=================================================================
IF ( IDTDMS /= 0 ) CALL SRCDMS( STT(:,:,:,IDTDMS) )
IF ( IDTSO2 /= 0 ) CALL SRCSO2( STT(:,:,:,IDTSO2), NSEASON )
IF ( IDTSO4 /= 0 ) CALL SRCSO4( STT(:,:,:,IDTSO4) )
IF ( IDTNH3 /= 0 ) CALL SRCNH3( STT(:,:,:,IDTNH3) )
! Return to calling program
END SUBROUTINE EMISSSULFATE
!-----------------------------------------------------------------------------
SUBROUTINE SRCDMS( TC ) 1,11
!
!*****************************************************************************
! Subroutine SRCDMS, from Mian Chin's GOCART model, add DMS emissions
! to the tracer array. Modified for use with the GEOS-CHEM model.
! (bmy, 6/2/00, 8/16/05)
!
! Arguments as Input/Output:
! ===========================================================================
! (1 ) TC (REAL*8 ) : Initial tracer mass [kg], plus DMS emissions
!
! NOTES:
! (1 ) Now reference NSRCE, LWI, DXYP, XTRA2 from either header files or
! F90 modules. Now use routines from "pressure_mod.f" to compute
! grid box surface pressures. (bmy, 9/18/02)
! (2 ) Now replace DXYP(J) with routine GET_AREA_M2 of "grid_mod.f"
! Now use routine GET_TS_EMIS from the new "time_mod.f". (bmy, 3/27/03)
! (3 ) For GEOS-4, convert PBL from [m] to [hPa] w/ the hydrostatic law.
! Now references SCALE_HEIGHT from "CMN_GCTM". Added BLTHIK variable
! for PBL thickness in [hPa]. (bmy, 1/15/04)
! (4 ) Remove reference to "pressure_mod.f". Now reference GET_FRAC_OF_PBL
! and GET_PBL_TOP_L from "pbl_mix_mod.f". (bmy, 2/22/05)
! (5 ) Switch from Liss & Merlivat to Nightingale formulation for DMS
! emissions. (swu, bmy, 8/16/05)
!******************************************************************************
!
! Reference to diagnostic arrays
USE DIAG_MOD, ONLY : AD13_DMS
USE DAO_MOD, ONLY : IS_WATER, LWI, PBL
USE GRID_MOD, ONLY : GET_AREA_M2
USE PBL_MIX_MOD, ONLY : GET_FRAC_OF_PBL, GET_PBL_TOP_L
USE TIME_MOD, ONLY : GET_TS_EMIS
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND13 (for now)
# include "CMN_GCTM" ! SCALE_HEIGHT
! Arguments
REAL*8, INTENT(INOUT) :: TC(IIPAR,JJPAR,LLPAR)
! Local variables
INTEGER :: I, J, L, NTOP
REAL*8 :: DTSRCE, SST, Sc, CONC, W10
REAL*8 :: ScCO2, AKw, ERATE, DMSSRC, FEMIS
! Molecular weight of DMS, kg/mole
REAL*8, PARAMETER :: DMS_MW = 62d0
! Ratio of molecular weights: S/DMS
REAL*8, PARAMETER :: S_DMS = 32d0 / 62d0
! External functions
REAL*8, EXTERNAL :: SFCWINDSQR
!=================================================================
! SRCDMS begins here!
!=================================================================
! Chemistry timestep in seconds
DTSRCE = GET_TS_EMIS() * 60d0
!=================================================================
! Compute DMS emissions = seawater DMS * transfer velocity
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, SST, Sc, CONC, W10, ScCO2 )
!$OMP+PRIVATE( AKw, ERATE, DMSSRC, NTOP, L, FEMIS )
!$OMP+SCHEDULE( DYNAMIC )
DO J = 1, JJPAR
DO I = 1, IIPAR
! Sea surface temperature in Celsius
SST = SSTEMP(I,J) - 273.15d0
! Only do the following for water boxes
IF ( IS_WATER(I,J) ) THEN
! Schmidt number for DMS (Saltzman et al., 1993)
Sc = 2674.0d0 - 147.12d0*SST +
& 3.726d0*(SST**2) - 0.038d0*(SST**3)
!===========================================================
! Calculate transfer velocity in cm/hr (AKw)
!
! Tans et al. transfer velocity (1990) for CO2 at
! 25oC (Erickson, 1993)
!
! Tans et al. assumed AKW=0 when W10<=3. I modified it
! to let DMS emit at low windseeds too. Chose 3.6m/s as
! the threshold.
!
! Schmidt number for CO2: Sc = 600 (20oC, fresh water)
! Sc = 660 (20oC, seawater)
! Sc = 428 (25oC, Erickson 93)
!===========================================================
CONC = DMSo(I,J)
W10 = SQRT( SFCWINDSQR(I,J) )
!-----------------------------------------------------------
! Tans et al. (1990)
!ScCO2 = 428.d0
!IF (W10 .LE. 3.6) THEN
! AKw = 1.0667d0 * W10
!ELSE
! AKw = 6.4d0 * (W10 - 3.d0)
!ENDIF
!-----------------------------------------------------------
! Wanninkhof (1992)
!ScCO2 = 660.d0
!AKw = 0.31d0 * W10**2
!-----------------------------------------------------------
!! Liss and Merlivat (1986)
!ScCO2 = 600.d0
!IF ( W10 <= 3.6d0 ) then
! AKw = 0.17d0 * W10
!
!ELSE IF ( W10 <= 13.d0 ) THEN
! AKw = 2.85d0 * W10 - 9.65d0
!
!ELSE
! AKw = 5.90d0 * W10 - 49.3d0
!
!ENDIF
!-----------------------------------------------------------
! NOTE: Also need to uncomment this section if using
! Tans, Wanninkhof, or Liss & Merlivat
!IF ( W10 <= 3.6d0 ) THEN
! AKw = AKw * ( (ScCO2/Sc)**0.667 )
!ELSE
! AKw = AKw * SQRT(ScCO2/Sc)
!ENDIF
!-----------------------------------------------------------
! Nightingale [2000] (swu, bmy, 8/16/05)
!
! Note that from Nightingale et al [2000a],
! the best fit formulation should be:
!
! AKw = ( 0.222*W10*W10 + 0.333*W10 ) * sqrt( ScCO2/Sc )
!
! But from Nightingale et al [2000b], which reported that
! more measurements were incorported, they claimed that
! the following is the best fit:
!
ScCO2 = 600.d0
AKw = ( 0.24d0*W10*W10 + 0.061d0*W10 ) * SQRT( ScCO2/Sc )
!-----------------------------------------------------------
!===========================================================
! Calculate emission flux in kg/box/timestep
!
! AKw is in cm/hr : AKw/100/3600 is m/sec.
! CONC is in nM/L (nM/dm3) : CONC*1E-12*1000 is kmole/m3.
! DMS_MW is in g DMS/mol = kg/kmole
! ERATE is in kg DMS/m2/timestep
! DMSSRC is in kg DMS/box/timestep
!===========================================================
ERATE = ( AKw / 100.d0 / 3600.d0 ) *
& ( CONC * 1.d-12 * 1000.d0 ) * DMS_MW * DTSRCE
DMSSRC = ERATE * GET_AREA_M2( J )
!===========================================================
! Add DMS emissions [kg DMS/box] into the tracer array
!===========================================================
! Top layer of the PBL
NTOP = CEILING( GET_PBL_TOP_L( I, J ) )
! Loop thru the boundary layer
DO L = 1, NTOP
! Fraction of PBL spanned by grid box (I,J,L) [unitless]
FEMIS = GET_FRAC_OF_PBL( I, J, L )
! DMS in box (I,J,L) plus emissions [kg]
TC(I,J,L) = TC(I,J,L) + ( FEMIS * DMSSRC )
ENDDO
ELSE
! If we are not over water, then there is no DMS source
DMSSRC = 0.d0
ENDIF
!==============================================================
! ND13 diagnostic: DMS emissions [kg S/box/timestep]
!==============================================================
IF ( ND13 > 0 ) THEN
AD13_DMS(I,J) = AD13_DMS(I,J) + ( DMSSRC * S_DMS ) ! / DTSRCE
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE SRCDMS
!------------------------------------------------------------------------------
SUBROUTINE SRCSO2( TC, NSEASON ) 1,28
!
!******************************************************************************
! Subroutine SRCSO2 (originally from Mian Chin) computes SO2 emissons from
! aircraft, biomass, and anthro sources. (rjp, bdf, bmy, 6/2/00, 6/26/06)
!
! Arguments as Input/Output:
! ===========================================================================
! (1 ) NSEASON (INTEGER) : Season number: 1=DJF; 2=MAM; 3=JJA; 4=SON
! (2 ) TC (REAL*8 ) : SO2 tracer mass [kg]
!
! NOTES:
! (1 ) Now reference NSRCE, JDAY, PBL, XTRA2, BXHEIGHT from either header
! files or F90 modules. Also use routines from "pressure_mod.f" to
! compute grid box pressures. (bmy, 9/18/02)
! (2 ) Now use routines GET_TS_EMIS and GET_DAY_OF_YEAR from the new
! "time_mod.f" (bmy, 3/27/03)
! (3 ) For GEOS-4, convert PBL from [m] to [hPa] w/ the hydrostatic law.
! Now references SCALE_HEIGHT from "CMN_GCTM". Added BLTHIK variable
! to hold PBL thickness in [hPa]. (bmy, 1/15/04)
! (4 ) Now references AD13_SO2_sh array from "diag_mod.f". Also references
! LSHIPSO2 from "CMN_SETUP" (bec, bmy, 5/20/04)
! (5 ) Now references LSHIPSO2 from "logical_mod.f" (bmy, 7/20/04)
! (6 ) Now references routines GET_EPA_ANTHRO and GET_USA_MASK from
! "epa_nei_mod.f". Now references GET_AREA_CM2 from "grid_mod.f".
! Now references GET_DAY_OF_WEEK from "time_mod.f" Now references
! LNEI99 from "logical_mod.f". Now can overwrite the anthro SOx
! emissions over the continental US if LNEI99=T. Now references IDTSO2
! from "tracerid_mod.f. (rch, rjp, bmy, 11/16/04)
! (7 ) Remove reference to "pressure_mod.f". Now reference GET_FRAC_OF_PBL
! and GET_PBL_TOP_L from "pbl_mix_mod.f". Removed reference to header
! file CMN. (bmy, 2/22/05)
! (8 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
! (9 ) Now references GET_BRAVO_ANTHRO and GET_BRAVO_MASK from "bravo_mod.f"
! for BRAVO Mexican emissions. (rjp, kfb, bmy, 6/26/06)
!******************************************************************************
!
! Reference to diagnostic arrays
USE BRAVO_MOD, ONLY : GET_BRAVO_ANTHRO, GET_BRAVO_MASK
USE DIAG_MOD, ONLY : AD13_SO2_an, AD13_SO2_ac
USE DIAG_MOD, ONLY : AD13_SO2_bb, AD13_SO2_nv
USE DIAG_MOD, ONLY : AD13_SO2_ev, AD13_SO2_bf
USE DIAG_MOD, ONLY : AD13_SO2_sh
USE DAO_MOD, ONLY : BXHEIGHT, PBL
USE EPA_NEI_MOD, ONLY : GET_EPA_ANTHRO, GET_EPA_BIOFUEL
USE EPA_NEI_MOD, ONLY : GET_USA_MASK
USE ERROR_MOD, ONLY : ERROR_STOP
USE GRID_MOD, ONLY : GET_AREA_CM2
USE LOGICAL_MOD, ONLY : LBRAVO, LNEI99, LSHIPSO2
USE PBL_MIX_MOD, ONLY : GET_FRAC_OF_PBL, GET_PBL_TOP_L
USE TIME_MOD, ONLY : GET_TS_EMIS, GET_DAY_OF_YEAR
USE TIME_MOD, ONLY : GET_DAY_OF_WEEK
USE TRACER_MOD, ONLY : XNUMOL
USE TRACERID_MOD, ONLY : IDTSO2
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND13, LD13 (for now)
# include "CMN_GCTM" ! SCALE_HEIGHT
! Arguments
INTEGER, INTENT(IN) :: NSEASON
REAL*8, INTENT(INOUT) :: TC(IIPAR,JJPAR,LLPAR)
! Local variables
LOGICAL :: WEEKDAY
INTEGER :: I, J, K, L, LV1, LV2, NTOP, JDAY
INTEGER :: DAY_NUM
REAL*8 :: ZH(0:LLPAR), DZ(LLPAR), SO2(LLPAR)
REAL*8 :: DTSRCE, HGHT, SO2SRC
REAL*8 :: SLAB, SLAB1
REAL*8 :: TSO2, FEMIS
REAL*8 :: AREA_CM2, AN, BF
REAL*8 :: SO2an(IIPAR,JJPAR,2)
REAL*8 :: SO2bf(IIPAR,JJPAR)
! Ratio of molecular weights: S/SO2
REAL*8, PARAMETER :: S_SO2 = 32d0 / 64d0
!=================================================================
! SRCSO2 begins here!
!================================================================
! DTSRCE is the emission timestep in seconds
DTSRCE = GET_TS_EMIS() * 60d0
! JDAY is the day of year (0-365 or 0-366)
JDAY = GET_DAY_OF_YEAR()
! Get current day of the week
DAY_NUM = GET_DAY_OF_WEEK()
! Is it a weekday?
WEEKDAY = ( DAY_NUM > 0 .and. DAY_NUM < 6 )
!=================================================================
! SO2 emissions from non-eruptive volcanoes [kg SO2/box/s].
! Assume that emission only occurs at the crater altitude.
!=================================================================
IF ( LENV ) THEN
! Initialize
ESO2_nv = 0.d0
! Loop thru each non-erupting volcano
DO K = 1, NNVOL
! Elevation of volcano crater
HGHT = DBLE( IELVn(k) )
! Altitude of crater from the ground
ZH(0) = 0.d0
! Loop over levels
DO L = 1, LLPAR
! Thickness of layer [m] w/ crater
DZ(L) = BXHEIGHT(INV(K),JNV(K),L)
! Increment altitude
ZH(L) = ZH(L-1) + DZ(L)
! If we are at the crater altitude, add emissions and exit
IF ( ZH(L-1) <= HGHT .AND. ZH(L) > HGHT ) THEN
ESO2_nv(INV(K),JNV(K),L) =
& ESO2_nv(INV(K),JNV(K),L) + Env(K)
EXIT
ENDIF
ENDDO
ENDDO
ENDIF
!=================================================================
! Calculate eruptive volcanic emission of SO2.
!=================================================================
IF ( LEEV ) THEN
! Initialize
ESO2_ev = 0.D0
! Loop thru each erupting volcano
DO K = 1, NEVOL
! Test to see if the volcano is erupting
IF ( JDAY < IDAYS(K) .OR. JDAY > IDAYe(K) ) GOTO 20
!===========================================================
! Define a slab at the top 1/3 of the volcano plume.
!===========================================================
HGHT = DBLE( IHGHT(K) )
! slab bottom height
SLAB1 = HGHT - ( HGHT - DBLE ( IELVe(K) ) ) / 3.d0
! Slab thickness
SLAB = HGHT - SLAB1
ZH(0) = 0.d0
! Loop over each level
DO L = 1, LLPAR
! DZ is the thickness of level L [m]
DZ(L) = BXHEIGHT(IEV(K),JEV(K),L)
! ZH is the height of the top edge of
! level L, measured from the ground up [m]
ZH(L) = ZH(L-1) + DZ(L)
! max model erup.height
IF ( L == LLPAR .AND. HGHT > ZH(L) ) THEN
LV2 = LLPAR
!HGHT = ZH(L)
!SLAB1 = SLAB1 - ( HGHT - ZH(L) )
ENDIF
!========================================================
! If Zh(l) <= bottom of the slab, go to next level.
!========================================================
IF ( ZH(L) <= SLAB1 ) GOTO 22
!========================================================
! If the slab is only in current level: CASE 1
! ---------------------------------- ZH(L)
! HGHT ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! SLAB1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! ---------------------------------- ZH(L-1)
!========================================================
IF ( ZH(L-1) <= SLAB1 .AND. ZH(L) >= HGHT ) THEN
LV1 = L
LV2 = L
DZ(L) = SLAB
!========================================================
! The slab extends more then one level. Find the
! lowest (lv1) and the highest (lv2) levels:
! --------------------------------- ZH(L)
! HGHT ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! --------------------------------- ZH(L-1)
!
! --------------------------------- ZH(L)
! SLAB1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! --------------------------------- ZH(L-1)
!========================================================
ELSE IF (ZH(L-1) <= SLAB1 .AND. ZH(L) > SLAB1) THEN
LV1 = L
DZ(L) = ZH(L) - SLAB1
ELSE IF (ZH(L-1) < HGHT .AND. ZH(L) > HGHT ) THEN
LV2 = L
DZ(L) = HGHT - ZH(L-1)
EXIT ! do 20
ENDIF
! Go to next level
22 CONTINUE
ENDDO
!===========================================================
! Calculate SO2 emission in the levels between LV1 and LV2.
! Convert Eev from [kg SO2/box/event] to [kg SO2/box/s].
! ESO2_ev is distributed evenly with altitude among the slab.
!===========================================================
DO L = LV1, LV2
ESO2_ev(IEV(K),JEV(K),L) = ESO2_ev(IEV(K),JEV(K),L) +
& EEV(K) / ( (IDAYe(K)-IDAYs(K)+1) * 24.d0 * 3600.d0 )
& * DZ(L) / SLAB
ENDDO
! Go to next volcano
20 CONTINUE
ENDDO
ENDIF
!=================================================================
! Overwrite USA w/ EPA/NEI99 (anthro+biofuel) SO2 emissions
! Overwrite MEXICO w/ BRAVO (anthro only ) SO2 emissions
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, AREA_CM2, AN, BF )
!$OMP+SCHEDULE( DYNAMIC )
DO J = 1, JJPAR
! Initialize
AN = 0d0
BF = 0d0
! Grid box surface area [cm2]
AREA_CM2 = GET_AREA_CM2(J)
DO I = 1, IIPAR
!-----------------------------------------------------------
! If we are using EPA/NEI99 (anthro + biofuel) ...
!-----------------------------------------------------------
IF ( LNEI99 ) THEN
! If we are over the USA ...
IF ( GET_USA_MASK( I, J ) > 0d0 ) THEN
! Read EPA/NEI99 SO2 emissions in [molec/cm2/s]
AN = GET_EPA_ANTHRO( I, J, IDTSO2, WEEKDAY )
BF = GET_EPA_BIOFUEL(I, J, IDTSO2, WEEKDAY )
! Convert anthro SO2 from [molec/cm2/s] to [kg/box/s]
! Place all anthro SO2 into surface layer
SO2an(I,J,1) = AN * AREA_CM2 / XNUMOL(IDTSO2)
SO2an(I,J,2) = 0d0
! Convert anthro SO2 from [molec/cm2/s] to [kg/box/s]
SO2bf(I,J) = BF * AREA_CM2 / XNUMOL(IDTSO2)
ELSE
! If we are not over the USA, then just use the regular
! emissions from ESO2_an and ESO2_bf (bmy, 11/16/04)
SO2an(I,J,1) = ESO2_an(I,J,1)
SO2an(I,J,2) = ESO2_an(I,J,2)
SO2bf(I,J) = ESO2_bf(I,J)
ENDIF
ELSE
! If we are not using EPA/NEI99 emissions, then just copy
! ESO2_an and ESO2_bf into local arrays (bmy, 11/16/04)
SO2an(I,J,1) = ESO2_an(I,J,1)
SO2an(I,J,2) = ESO2_an(I,J,2)
SO2bf(I,J) = ESO2_bf(I,J)
ENDIF
!-----------------------------------------------------------
! If we are using BRAVO emissions over Mexico ...
!-----------------------------------------------------------
IF ( LBRAVO ) THEN
! If we are over Mexico ...
IF ( GET_BRAVO_MASK( I, J ) > 0d0 ) THEN
! Read BRAVO SO2 emissions in [molec/cm2/s]
AN = GET_BRAVO_ANTHRO( I, J, IDTSO2 )
! Convert anthro SO2 from [molec/cm2/s] to [kg/box/s]
! Place all anthro SO2 into surface layer
SO2an(I,J,1) = AN * AREA_CM2 / XNUMOL(IDTSO2)
SO2an(I,J,2) = 0d0
ELSE
! If we are not over MEXICO, then just use
! the regular emissions from ESO2_an
SO2an(I,J,1) = ESO2_an(I,J,1)
SO2an(I,J,2) = ESO2_an(I,J,2)
ENDIF
ELSE
! If we are not using BRAVO emissions, then just copy
! ESO2_an and ESO2_bf into local arrays
SO2an(I,J,1) = ESO2_an(I,J,1)
SO2an(I,J,2) = ESO2_an(I,J,2)
SO2bf(I,J) = ESO2_bf(I,J)
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Add SO2 emissions into model levels
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, NTOP, L, SO2, TSO2, FEMIS, SO2SRC )
!$OMP+SCHEDULE( DYNAMIC )
DO J = 1, JJPAR
DO I = 1, IIPAR
! Top of the boundary layer
NTOP = CEILING( GET_PBL_TOP_L( I, J ) )
! Zero SO2 array
DO L = 1, LLPAR
SO2(L) = 0d0
ENDDO
! Sum of anthro (surface + 100m), biomass, biofuel SO2 at (I,J)
TSO2 = SUM( SO2an(I,J,:) ) + ESO2_bb(I,J) + SO2bf(I,J)
! Also add SO2 from ship exhaust if necessary (bec, bmy, 5/20/04)
IF ( LSHIPSO2 ) TSO2 = TSO2 + ESO2_sh(I,J)
! Zero SO2SRC
SO2SRC = 0d0
!===============================================================
! Partition the total anthro and biomass SO2 emissions thru
! the entire boundary layer (if PBL top is higher than level 2)
!===============================================================
IF ( NTOP > 2 ) THEN
! Loop thru levels in the PBL
DO L = 1, NTOP
! Fraction of PBL spanned by grid box (I,J,L) [unitless]
FEMIS = GET_FRAC_OF_PBL( I, J, L )
! Partition total SO2 into level K
SO2(L) = FEMIS * TSO2
ENDDO
!===============================================================
! If PBL top occurs lower than or close to the top of level 2,
! then then surface SO2 goes into level 1 and the smokestack
! stack SO2 goes into level 2.
!===============================================================
ELSE
SO2(1) = SO2an(I,J,1) + ESO2_bb(I,J) + SO2bf(I,J)
SO2(2) = SO2an(I,J,2)
! Also add ship exhaust SO2 into surface if necessary
! (bec, bmy, 5/20/04)
IF ( LSHIPSO2 ) SO2(1) = SO2(1) + ESO2_sh(I,J)
ENDIF
! Error check
IF ( ABS( SUM( SO2 ) - TSO2 ) > 1.D-5 ) THEN
!$OMP CRITICAL
PRINT*, '### ERROR in SRCSO2!'
PRINT*, '### I, J, L, : ', I, J, L
PRINT*, '### SUM(SO2) : ', SUM( SO2 )
PRINT*, '### TSO2 : ', TSO2
!$OMP END CRITICAL
CALL ERROR_STOP( 'Check SO2 redistribution!',
& 'SRCSO2 (sulfate_mod.f)' )
ENDIF
!==============================================================
! Add anthro SO2, aircraft SO2, volcano SO2, and biomass SO2
! Convert from [kg SO2/box/s] -> [kg SO2/box/timestep]
!==============================================================
DO L = 1, LLPAR
! SO2 emissions [kg/box/s]
SO2SRC = SO2(L) + ESO2_ac(I,J,L) +
& ESO2_nv(I,J,L) + ESO2_ev(I,J,L)
! Add SO2 to TC array [kg/box/timestep]
TC(I,J,L) = TC(I,J,L) + ( SO2SRC * DTSRCE )
ENDDO
!==============================================================
! ND13 Diagnostic: SO2 emissions in [kg S/box/timestep]
!==============================================================
IF ( ND13 > 0 ) THEN
! Anthropogenic SO2 -- Levels 1-2
DO L = 1, 2
AD13_SO2_an(I,J,L) = AD13_SO2_an(I,J,L) +
& ( SO2an(I,J,L) * S_SO2 * DTSRCE )
ENDDO
! SO2 from biomass burning
AD13_SO2_bb(I,J) = AD13_SO2_bb(I,J) +
& ( ESO2_bb(I,J) * S_SO2 * DTSRCE )
! SO2 from biofuel burning
AD13_SO2_bf(I,J) = AD13_SO2_bf(I,J) +
& ( SO2bf(I,J) * S_SO2 * DTSRCE )
! SO2 from ship emissions (bec, bmy, 5/20/04)
IF ( LSHIPSO2 ) THEN
AD13_SO2_sh(I,J) = AD13_SO2_sh(I,J) +
& ( ESO2_sh(I,J) * S_SO2 * DTSRCE )
ENDIF
! Loop thru LD13 levels
DO L = 1, LD13
! SO2 from aircraft emissions
AD13_SO2_ac(I,J,L) = AD13_SO2_ac(I,J,L) +
& ( ESO2_ac(I,J,L) * S_SO2 * DTSRCE )
! SO2 from non-eruptive volcanoes
AD13_SO2_nv(I,J,L) = AD13_SO2_nv(I,J,L) +
& ( ESO2_nv(I,J,L) * S_SO2 * DTSRCE )
! SO2 from eruptive volcanoes
AD13_SO2_ev(I,J,L) = AD13_SO2_ev(I,J,L) +
& ( ESO2_ev(I,J,L) * S_SO2 * DTSRCE )
ENDDO
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE SRCSO2
!-----------------------------------------------------------------------------
SUBROUTINE SRCSO4( TC ) 1,20
!
!******************************************************************************
! Subroutine SRCSO4 (originally from Mian Chin) computes SO4 emissions from
! anthropogenic sources (rjp, bdf, bmy, 6/2/00, 10/25/05)
!
! Arguments as Input/Output:
! ===========================================================================
! (2) TC (REAL*8 ) : Array for SO4 mass [kg]
!
! NOTES:
! (1 ) Emission of SO4 is read in SULFATE_READYR, in [kg/box/s].
! It is converted to [kg/box/timestep] here.
! (2 ) Now use routine GET_TS_EMIS from the new "time_mod.f" (bmy, 3/27/03)
! (3 ) For GEOS-4, convert PBL from [m] to [hPa] w/ the barometric law.
! Now references SCALE_HEIGHT from "CMN_GCTM". Added BLTHIK variable
! to hold PBL thickness in [hPa]. (bmy, 1/15/04)
! (4 ) Now references GET_EPA_ANTHRO, GET_EPA_BIOFUEL, and GET_USA_MASK from
! "epa_nei_mod.f". Now references AD13_SO4_bf from "diag_mod.f". Now
! references GET_AREA_CM2 from "grid_mod.f". Now references
! GET_DAY_OF_WEEK from "time_mod.f". Now references LNEI99 from
! "logical_mod.f". Now can overwrite the anthro SOx emissions over
! the continental US if LNEI99=T. Now references IDTSO4 from
! "tracerid_mod.f". (rch, rjp, bmy, 11/16/04)
! (5 ) Remove reference to "pressure_mod.f". Now reference GET_FRAC_OF_PBL
! and GET_PBL_TOP_L from "pbl_mix_mod.f". Removed reference to header
! file CMN. (bmy, 2/22/05)
! (6 ) Now references XNUMOL & XNUMOLAIR from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! Reference to diagnostic arrays
USE DAO_MOD, ONLY : PBL
USE DIAG_MOD, ONLY : AD13_SO4_an, AD13_SO4_bf
USE EPA_NEI_MOD, ONLY : GET_EPA_ANTHRO, GET_EPA_BIOFUEL
USE EPA_NEI_MOD, ONLY : GET_USA_MASK
USE ERROR_MOD, ONLY : ERROR_STOP
USE GRID_MOD, ONLY : GET_AREA_CM2
USE LOGICAL_MOD, ONLY : LNEI99
USE PBL_MIX_MOD, ONLY : GET_FRAC_OF_PBL, GET_PBL_TOP_L
USE TIME_MOD, ONLY : GET_DAY_OF_WEEK, GET_TS_EMIS
USE TRACER_MOD, ONLY : XNUMOL
USE TRACERID_MOD, ONLY : IDTSO4
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND13 (for now)
# include "CMN_GCTM" ! SCALE_HEIGHT
! Arguments
REAL*8, INTENT(INOUT) :: TC(IIPAR,JJPAR,LLPAR)
! Local variables
LOGICAL :: WEEKDAY
INTEGER :: I, J, K, L, DAY_NUM, NTOP
REAL*8 :: SO4(LLPAR), DTSRCE
REAL*8 :: TSO4, FEMIS
REAL*8 :: AREA_CM2, EPA_AN, EPA_BF
REAL*8 :: SO4an(IIPAR,JJPAR,2)
REAL*8 :: SO4bf(IIPAR,JJPAR)
! Ratio of molecular weights: S/SO4
REAL*8, PARAMETER :: S_SO4 = 32d0 / 96d0
!=================================================================
! SRCSO4 begins here!
!=================================================================
! DTSRCE is the emission timestep in seconds
DTSRCE = GET_TS_EMIS() * 60d0
! Get current day of the week
DAY_NUM = GET_DAY_OF_WEEK()
! Is it a weekday?
WEEKDAY = ( DAY_NUM > 0 .and. DAY_NUM < 6 )
!=================================================================
! Overwrite USA w/ EPA/NEI99 SO4 emissions (if necessary)
! Store emissions into local arrays SO4an, SO4bf
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, AREA_CM2, EPA_AN, EPA_BF )
DO J = 1, JJPAR
! Grid box surface area [cm2]
AREA_CM2 = GET_AREA_CM2(J)
DO I = 1, IIPAR
! If we are using EPA/NEI99 emissions
IF ( LNEI99 ) THEN
! If we are over the USA ...
IF ( GET_USA_MASK( I, J ) > 0d0 ) THEN
! Read SO4 emissions in [molec/cm2/s]
EPA_AN = GET_EPA_ANTHRO( I, J, IDTSO4, WEEKDAY )
EPA_BF = GET_EPA_BIOFUEL(I, J, IDTSO4, WEEKDAY )
! Convert anthro SO4 from [molec/cm2/s] to [kg/box/s]
! Place all EPA/NEI99 anthro SO4 into surface layer
SO4an(I,J,1) = EPA_AN * AREA_CM2 / XNUMOL(IDTSO4)
SO4an(I,J,2) = 0d0
! Convert biofuel SO4 from [molec/cm2/s] to [kg/box/s]
SO4bf(I,J) = EPA_BF * AREA_CM2 / XNUMOL(IDTSO4)
ELSE
! If we are not over the USA, then just use the regular
! emissions from ESO4_an. Also set biofuel SO4 to zero
! since we currently don't read this in. (bmy, 11/16/04)
SO4an(I,J,1) = ESO4_an(I,J,1)
SO4an(I,J,2) = ESO4_an(I,J,2)
SO4bf(I,J) = 0d0
ENDIF
ELSE
! If we are not using EPA/NEI99 emissions, then just copy
! ESO4_an into ESO4 array. Also set biofuel SO4 to zero
! since we currently don't read this in. (bmy, 11/16/04)
SO4an(I,J,1) = ESO4_an(I,J,1)
SO4an(I,J,2) = ESO4_an(I,J,2)
SO4bf(I,J) = 0d0
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Compute SO4 emissions
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, NTOP, SO4, TSO4, L, FEMIS )
!$OMP+SCHEDULE( DYNAMIC )
DO J = 1, JJPAR
DO I = 1, IIPAR
! Top level of boundary layer at (I,J)
NTOP = CEILING( GET_PBL_TOP_L( I, J ) )
! Zero SO4 array at all levels
DO L = 1, LLPAR
SO4(L) = 0.0
ENDDO
! Compute total anthro SO4 (surface + 100m) plus biofuel SO4
TSO4 = SUM( SO4an(I,J,:) ) + SO4bf(I,J)
!==============================================================
! Partition the total anthro SO4 emissions thru the entire
! boundary layer (if PBL top is higher than level 2)
!==============================================================
IF ( NTOP > 2 ) THEN
! Loop thru boundary layer
DO L = 1, NTOP
! Fraction of PBL spanned by grid box (I,J,L) [unitless]
FEMIS = GET_FRAC_OF_PBL( I, J, L )
! Fraction of total SO4 in layer L
SO4(L) = FEMIS * TSO4
ENDDO
!==============================================================
! If PBL height is low and lower or similar to the second
! model layer then surface emission is emitted to the first
! model layer and the stack emission goes to the second model
! layer. Also add biofuel SO4 into the surface layer.
!==============================================================
ELSE
SO4(1) = SO4an(I,J,1) + SO4bf(I,J)
SO4(2) = SO4an(I,J,2)
ENDIF
IF ( ABS( SUM( SO4 ) - TSO4 ) > 1.D-5 ) THEN
!$OMP CRITICAL
PRINT*, '### ERROR in SRCSO4!'
PRINT*, '### I, J, L, : ', I, J, L
PRINT*, '### SUM(SO4) : ', SUM( SO4 )
PRINT*, '### TSO4 : ', TSO4
!$OMP END CRITICAL
CALL ERROR_STOP( 'Check SO4 redistribution',
& 'SRCSO4 (sulfate_mod.f)' )
ENDIF
!=============================================================
! Add SO4 emissions to tracer array
! Convert from [kg SO4/box/s] -> [kg SO4/box/timestep]
!=============================================================
DO L = 1, LLPAR
TC(I,J,L) = TC(I,J,L) + ( SO4(L) * DTSRCE )
ENDDO
!==============================================================
! ND13 Diagnostic: SO4 emission in [kg S/box/timestep]
!==============================================================
IF ( ND13 > 0 ) THEN
! Anthro SO4
DO L = 1, 2
AD13_SO4_an(I,J,L) = AD13_SO4_an(I,J,L) +
& ( SO4an(I,J,L) * S_SO4 * DTSRCE )
ENDDO
! Biofuel SO4
AD13_SO4_bf(I,J) = AD13_SO4_bf(I,J) +
& ( SO4bf(I,J) * S_SO4 * DTSRCE )
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE SRCSO4
!------------------------------------------------------------------------------
SUBROUTINE SRCNH3( TC ) 1,16
!
!******************************************************************************
! Subroutine SRCNH3 handles NH3 emissions into the GEOS-CHEM tracer array.
! (rjp, bmy, 12/17/01, 2/22/05)
!
! Arguments as Input/Output
! ============================================================================
! (1 ) TC (REAL*8 ) : Array for NH3 tracer mass in kg
!
! NOTES:
! (1 ) Now save NH3 emissions to ND13 diagnostic (bmy, 12/13/02)
! (2 ) Now reference AD13_NH3_na from "diag_mod.f", and archive natural
! source NH3 diagnostics for ND13. Also consider natural source NH3
! when partitioning by level into the STT array. (rjp, bmy, 3/23/03)
! (3 ) Now use routine GET_TS_EMIS from the new "time_mod.f" (bmy, 3/27/03)
! (4 ) For GEOS-4, convert PBL from [m] to [hPa] w/ the barometric law.
! Now references SCALE_HEIGHT from "CMN_GCTM". Added BLTHIK variable
! to hold PBL thickness in [hPa]. (bmy, 1/15/04)
! (5 ) Now references GET_EPA_ANTHRO, GET_EPA_BIOFUEL, and GET_USA_MASK from
! "epa_nei_mod.f". Now references GET_DAY_OF_WEEK from "time_mod.f".
! Now references LNEI99 from "logical_mod.f". Now references
! GET_AREA_CM2 from "grid_mod.f". Now references IDTNH3 from
! "tracerid_mod.f". Now references XNUMOL from CMN_O3. Now can
! overwrite the anthro & biofuel NH3 emissions over the continental US
! if LNEI99=T. Now references IDTNH3 from "tracerid_mod.f".
! (rjp, rch, bmy, 11/16/04)
! (6 ) Remove reference to "pressure_mod.f". Now reference GET_FRAC_OF_PBL
! and GET_PBL_TOP_L from "pbl_mix_mod.f". Removed reference to header
! file CMN. (bmy, 2/22/05)
! (7 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE DIAG_MOD, ONLY : AD13_NH3_an, AD13_NH3_bb
USE DIAG_MOD, ONLY : AD13_NH3_bf, AD13_NH3_na
USE DAO_MOD, ONLY : PBL
USE GRID_MOD, ONLY : GET_AREA_CM2
USE EPA_NEI_MOD, ONLY : GET_EPA_ANTHRO, GET_EPA_BIOFUEL
USE EPA_NEI_MOD, ONLY : GET_USA_MASK
USE ERROR_MOD, ONLY : ERROR_STOP
USE LOGICAL_MOD, ONLY : LNEI99
USE PBL_MIX_MOD, ONLY : GET_FRAC_OF_PBL, GET_PBL_TOP_L
USE TIME_MOD, ONLY : GET_DAY_OF_WEEK, GET_TS_EMIS
USE TRACER_MOD, ONLY : XNUMOL
USE TRACERID_MOD, ONLY : IDTNH3
# include "CMN_SIZE" ! Size parameters
# include "CMN_DIAG" ! ND13
# include "CMN_GCTM" ! SCALE_HEIGHT
! Argumetns
REAL*8, INTENT(INOUT) :: TC(IIPAR,JJPAR,LLPAR)
! Local variables
LOGICAL :: WEEKDAY
INTEGER :: I, J, L, K, NTOP, DAY_NUM
REAL*8 :: FEMIS, DTSRCE, TNH3
REAL*8 :: AREA_CM2, EPA_AN, EPA_BF
REAL*8 :: NH3an(IIPAR,JJPAR)
REAL*8 :: NH3bf(IIPAR,JJPAR)
!=================================================================
! SRCNH3 begins here!
!=================================================================
! Emission timestep [s]
DTSRCE = GET_TS_EMIS() * 60d0
! Get current day of the week
DAY_NUM = GET_DAY_OF_WEEK()
! Is it a weekday?
WEEKDAY = ( DAY_NUM > 0 .and. DAY_NUM < 6 )
!=================================================================
! Overwrite USA with EPA/NEI NH3 emissions (if necessary)
! Store emissions into local arrays NH3an, NH3bf
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!-----------------------------------------------------------------------------
! NOTE: There seems to be some problems with the EPA/NEI NH3 emissions.
! Therefore we will use the existing emissions for NH3 until further notice.
! Comment out the lines below until further notice. (bmy, 11/17/04)
!!$OMP+PRIVATE( I, J, AREA_CM2, EPA_AN, EPA_BF )
!-----------------------------------------------------------------------------
!$OMP+PRIVATE( I, J )
DO J = 1, JJPAR
!-------------------------------------------------------------------
! NOTE: There seems to be some problems with the EPA/NEI NH3
! emissions. Therefore we will use the existing emissions for NH3
! until further notice. Comment out the lines below until
! further notice. (bmy, 11/17/04)
!! Grid box surface area [cm2]
!AREA_CM2 = GET_AREA_CM2( J )
!-------------------------------------------------------------------
DO I = 1, IIPAR
!-----------------------------------------------------------------
! NOTE: There seems to be some problems with the EPA/NEI NH3
! emissions. Therefore we will use the existing emissions for
! NH3 until further notice. Comment out the lines below until
! further notice. (bmy, 11/17/04)
!
!! If we are using EPA/NEI99 emissions ...
!IF ( LNEI99 ) THEN
!
! ! If we are over the USA ...
! IF ( GET_USA_MASK( I, J ) > 0d0 ) THEN
!
! ! Read NH3 anthro emissions in [molec NH3/cm2/s]
! EPA_AN = GET_EPA_ANTHRO( I, J, IDTNH3, WEEKDAY )
! EPA_BF = GET_EPA_BIOFUEL( I, J, IDTNH3, WEEKDAY )
!
! ! Convert from [molec NH3/cm2/s] to [kg NH3/box/sec]
! NH3an(I,J) = EPA_AN * AREA_CM2 / XNUMOL(IDTNH3)
! NH3bf(I,J) = EPA_BF * AREA_CM2 / XNUMOL(IDTNH3)
!
! ELSE
!
! ! If we are not over the USA, just use the regular
! ! emissions in NH3_an and NH3bf (bmy, 11/16/04)
! NH3an(I,J) = ENH3_an(I,J)
! NH3bf(I,J) = ENH3_bf(I,J)
!
! ENDIF
!
!ELSE
!-----------------------------------------------------------------
! If we are not using the EPA/NEI emissions, just copy the
! regular ENH3_an and ENH3_bf to local arrays. (bmy, 11/16/04)
NH3an(I,J) = ENH3_an(I,J)
NH3bf(I,J) = ENH3_bf(I,J)
!-----------------------------------------------------------------
! NOTE: There seems to be some problems with the EPA/NEI NH3
! emissions. Therefore we will use the existing emissions for
! NH3 until further notice. Comment out the lines below until
! further notice. (bmy, 11/17/04)
!ENDIF
!-----------------------------------------------------------------
ENDDO
ENDDO
!$OMP END PARALLEL DO
!=================================================================
! Partition NH3 emissions into the STT tracer array
!=================================================================
! Loop over surface grid boxes
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, NTOP, TNH3, L, FEMIS )
!$OMP+SCHEDULE( DYNAMIC )
DO J = 1, JJPAR
DO I = 1, IIPAR
! Layer where the PBL top happens
NTOP = CEILING( GET_PBL_TOP_L( I, J ) )
! Sum all types of NH3 emission [kg/box/s]
TNH3 = NH3an(I,J) + ENH3_bb(I,J) +
& NH3bf(I,J) + ENH3_na(I,J)
!==============================================================
! Add NH3 emissions [kg NH3/box] into the tracer array
! Partition total NH3 throughout the entire boundary layer
!==============================================================
! Loop over all levels in the boundary layer
DO L = 1, NTOP
! Fraction of PBL spanned by grid box (I,J,L) [unitless]
FEMIS = GET_FRAC_OF_PBL( I, J, L )
! Add NH3 emissions into tracer array [kg NH3/timestep]
TC(I,J,L) = TC(I,J,L) + ( TNH3 * FEMIS * DTSRCE )
ENDDO
!============================================================
! ND13 diagnostics: NH3 emissions [kg NH3/box/timestep]
!============================================================
IF ( ND13 > 0 ) THEN
! Anthro NH3
AD13_NH3_an(I,J) = AD13_NH3_an(I,J) +
& ( NH3an(I,J) * DTSRCE )
! Biomass NH3
AD13_NH3_bb(I,J) = AD13_NH3_bb(I,J) +
& ( ENH3_bb(I,J) * DTSRCE )
! Biofuel NH3
AD13_NH3_bf(I,J) = AD13_NH3_bf(I,J) +
& ( NH3bf(I,J) * DTSRCE )
! Natural source NH3
AD13_NH3_na(I,J) = AD13_NH3_na(I,J) +
& ( ENH3_na(I,J) * DTSRCE )
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE SRCNH3
!------------------------------------------------------------------------------
FUNCTION GET_OH( I, J, L ) RESULT( OH_MOLEC_CM3 ) 5,26
!
!******************************************************************************
! Function GET_OH returns OH from SMVGEAR's CSPEC array (for coupled runs)
! or monthly mean OH (for offline runs). Imposes a diurnal variation on
! OH for offline simulations. (bmy, 12/16/02, 7/20/04)
!
! Arguments as Input:
! ============================================================================
! (1-3) I, J, L (INTEGER) : Grid box indices for lon, lat, vertical level
!
! NOTES:
! (1 ) We assume SETTRACE has been called to define IDOH (bmy, 11/1/02)
! (2 ) Now use function GET_TS_CHEM from "time_mod.f" (bmy, 3/27/03)
! (3 ) Now reference ITS_A_FULLCHEM_SIM, ITS_AN_AEROSOL_SIM from
! "tracer_mod.f". Also replace IJSURF w/ an analytic function.
! (bmy, 7/20/04)
!******************************************************************************
!
! References to F90 modules
USE COMODE_MOD, ONLY : CSPEC, JLOP
USE DAO_MOD, ONLY : SUNCOS
USE ERROR_MOD, ONLY : ERROR_STOP
USE GLOBAL_OH_MOD, ONLY : OH
USE TIME_MOD, ONLY : GET_TS_CHEM
USE TRACER_MOD, ONLY : ITS_A_FULLCHEM_SIM, ITS_AN_AEROSOL_SIM
USE TRACERID_MOD, ONLY : IDOH
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: I, J, L
! Local variables
INTEGER :: JLOOP
REAL*8 :: OH_MOLEC_CM3
!=================================================================
! GET_OH begins here!
!=================================================================
IF ( ITS_A_FULLCHEM_SIM() ) THEN
!---------------------
! Coupled simulation
!---------------------
! JLOOP = SMVGEAR 1-D grid box index
JLOOP = JLOP(I,J,L)
! Take OH from the SMVGEAR array CSPEC
! OH is defined only in the troposphere
IF ( JLOOP > 0 ) THEN
OH_MOLEC_CM3 = CSPEC(JLOOP,IDOH)
ELSE
OH_MOLEC_CM3 = 0d0
ENDIF
ELSE IF ( ITS_AN_AEROSOL_SIM() ) THEN
!---------------------
! Offline simulation
!---------------------
! 1-D grid box index for SUNCOS
JLOOP = ( (J-1) * IIPAR ) + I
! Test for sunlight...
IF ( SUNCOS(JLOOP) > 0d0 .and. TCOSZ(I,J) > 0d0 ) THEN
! Impose a diurnal variation on OH during the day
OH_MOLEC_CM3 = OH(I,J,L) *
& ( SUNCOS(JLOOP) / TCOSZ(I,J) ) *
& ( 1440d0 / GET_TS_CHEM() )
! Make sure OH is not negative
OH_MOLEC_CM3 = MAX( OH_MOLEC_CM3, 0d0 )
ELSE
! At night, OH goes to zero
OH_MOLEC_CM3 = 0d0
ENDIF
ELSE
!---------------------
! Invalid simulation
!---------------------
CALL ERROR_STOP( 'Invalid NSRCX!', 'GET_OH (sulfate_mod.f)')
ENDIF
! Return to calling program
END FUNCTION GET_OH
!------------------------------------------------------------------------------
SUBROUTINE SET_OH( I, J, L, OH ) 1,2
!
!******************************************************************************
! Function SET_OH saves the modified OH value back to SMVGEAR's CSPEC array
! for coupled sulfate/aerosol simulations. (bmy, 12/16/02)
!
! Arguments as Input:
! ============================================================================
! (1-3) I, J, L (INTEGER) : Grid box indices for lon, lat, vertical level
! (4 ) OH (REAL*8 ) : OH at grid box (I,J,L) to be saved into CSPEC
!
! NOTES:
! (1 ) We assume SETTRACE has been called to define IDOH (bmy, 12/16/02)
!******************************************************************************
!
! References to F90 modules
USE COMODE_MOD, ONLY : CSPEC, JLOP
USE TRACERID_MOD, ONLY : IDOH
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: I, J, L
REAL*8, INTENT(IN) :: OH
! Local variables
INTEGER :: JLOOP
!=================================================================
! SET_OH begins here!
!=================================================================
! JLOOP = SMVGEAR 1-D grid box index
JLOOP = JLOP(I,J,L)
! Replace OH into CSPEC(troposphere only)
IF ( JLOOP > 0 ) THEN
CSPEC(JLOOP,IDOH) = OH
ENDIF
! Return to calling program
END SUBROUTINE SET_OH
!------------------------------------------------------------------------------
FUNCTION GET_NO3( I, J, L ) RESULT( NO3_MOLEC_CM3 ) 2,20
!
!******************************************************************************
! Function GET_NO3 returns NO3 from SMVGEAR's CSPEC array (for coupled runs)
! or monthly mean OH (for offline runs). For offline runs, the concentration
! of NO3 is set to zero during the day. (rjp, bmy, 12/16/02)
!
! Arguments as Input:
! ============================================================================
! (1-3) I, J, L (INTEGER) : Grid box indices for lon, lat, vertical level
!
! NOTES:
! (1 ) Now references ERROR_STOP from "error_mod.f". We also assume that
! SETTRACE has been called to define IDNO3. Now also set NO3 to
! zero during the day. (rjp, bmy, 12/16/02)
! (2 ) Now reference ITS_A_FULLCHEM_SIM and ITS_AN_AEROSOL_SIM from
! "tracer_mod.f". Also remove reference to CMN. Also replace
! IJSURF with an analytic function. (bmy, 7/20/04)
!******************************************************************************
!
! References to F90 modules
USE COMODE_MOD, ONLY : CSPEC, JLOP
USE DAO_MOD, ONLY : AD, SUNCOS
USE ERROR_MOD, ONLY : ERROR_STOP
USE GLOBAL_NO3_MOD, ONLY : NO3
USE TRACER_MOD, ONLY : ITS_A_FULLCHEM_SIM, ITS_AN_AEROSOL_SIM
USE TRACERID_MOD, ONLY : IDNO3
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: I, J, L
! Local variables
INTEGER :: JLOOP
REAL*8 :: NO3_MOLEC_CM3
! External functions
REAL*8, EXTERNAL :: BOXVL
!=================================================================
! GET_NO3 begins here!
!=================================================================
IF ( ITS_A_FULLCHEM_SIM() ) THEN
!--------------------
! Coupled simulation
!--------------------
! 1-D SMVGEAR grid box index
JLOOP = JLOP(I,J,L)
! Take NO3 from the SMVGEAR array CSPEC
! NO3 is defined only in the troposphere
IF ( JLOOP > 0 ) THEN
NO3_MOLEC_CM3 = CSPEC(JLOOP,IDNO3)
ELSE
NO3_MOLEC_CM3 = 0d0
ENDIF
ELSE IF ( ITS_AN_AEROSOL_SIM() ) THEN
!==============================================================
! Offline simulation: Read monthly mean GEOS-CHEM NO3 fields
! in [v/v]. Convert these to [molec/cm3] as follows:
!
! vol NO3 moles NO3 kg air kg NO3/mole NO3
! ------- = --------- * -------- * ---------------- = kg NO3
! vol air moles air 1 kg air/mole air
!
! And then we convert [kg NO3] to [molec NO3/cm3] by:
!
! kg NO3 molec NO3 mole NO3 1 molec NO3
! ------ * --------- * -------- * ----- = ---------
! 1 mole NO3 kg NO3 cm3 cm3
! ^ ^
! |____________________|
! this is XNUMOL_NO3
!
! If at nighttime, use the monthly mean NO3 concentration from
! the NO3 array of "global_no3_mod.f". If during the daytime,
! set the NO3 concentration to zero. We don't have to relax to
! the monthly mean concentration every 3 hours (as for HNO3)
! since NO3 has a very short lifetime. (rjp, bmy, 12/16/02)
!==============================================================
! 1-D grid box index for SUNCOS
JLOOP = ( (J-1) * IIPAR ) + I
! Test if daylight
IF ( SUNCOS(JLOOP) > 0d0 ) THEN
! NO3 goes to zero during the day
NO3_MOLEC_CM3 = 0d0
ELSE
! At night: Get NO3 [v/v] and convert it to [kg]
NO3_MOLEC_CM3 = NO3(I,J,L) * AD(I,J,L) * ( 62d0/28.97d0 )
! Convert NO3 from [kg] to [molec/cm3]
NO3_MOLEC_CM3 = NO3_MOLEC_CM3 * XNUMOL_NO3 / BOXVL(I,J,L)
ENDIF
! Make sure NO3 is not negative
NO3_MOLEC_CM3 = MAX( NO3_MOLEC_CM3, 0d0 )
ELSE
!--------------------
! Invalid simulation
!--------------------
CALL ERROR_STOP( 'Invalid NSRCX!','GET_NO3 (sulfate_mod.f)')
ENDIF
! Return to calling program
END FUNCTION GET_NO3
!------------------------------------------------------------------------------
SUBROUTINE SET_NO3( I, J, L, NO3 ) 1,2
!
!******************************************************************************
! Function SET_NO3 saves the modified NO3 value back to SMVGEAR's CSPEC array
! for coupled lfate/aerosol simulations. (rjp, bmy, 12/16/02, 7/20/04)
!
! Arguments as Input:
! ============================================================================
! (1-3) I, J, L (INTEGER) : Grid box indices for lon, lat, vertical level
! (4 ) NO3 (REAL*8 ) : OH at grid box (I,J,L) to be saved into CSPEC
!
! NOTES:
! (1 ) We assume SETTRACE has been called to define IDNO3. (bmy, 12/16/02)
! (2 ) Remove references to "error_mod.f" and CMN (bmy, 7/20/04)
!******************************************************************************
!
! References to F90 modules
USE COMODE_MOD, ONLY : CSPEC, JLOP
USE TRACERID_MOD, ONLY : IDNO3
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: I, J, L
REAL*8, INTENT(IN) :: NO3
! Local variables
INTEGER :: JLOOP
!=================================================================
! SET_NO3 begins here!
!=================================================================
! 1-D grid box index for CSPEC
JLOOP = JLOP(I,J,L)
! Replace OH into CSPEC (troposphere only)
IF ( JLOOP > 0 ) THEN
CSPEC(JLOOP,IDNO3) = NO3
ENDIF
! Return to calling program
END SUBROUTINE SET_NO3
!------------------------------------------------------------------------------
FUNCTION GET_O3( I, J, L ) RESULT( O3_VV ) 3,25
!
!******************************************************************************
! Function GET_O3 returns monthly mean O3 for offline sulfate aerosol
! simulations. (bmy, 12/16/02, 10/25/05)
!
! Arguments as Input:
! ============================================================================
! (1-3) I, J, L (INTEGER) : Grid box indices for lon, lat, vertical level
!
! NOTES:
! (1 ) We assume SETTRACE has been called to define IDO3. (bmy, 12/16/02)
! (2 ) Now reference inquiry functions from "tracer_mod.f" (bmy, 7/20/04)
! (3 ) Now remove reference to CMN, it's obsolete. (bmy, 8/22/05)
! (4 ) Now references XNUMOLAIR from "tracer_mod.f" (bmy, 10/25/05)
!******************************************************************************
!
! References to F90 modules
USE COMODE_MOD, ONLY : CSPEC, JLOP, VOLUME
USE DAO_MOD, ONLY : AIRDEN
USE ERROR_MOD, ONLY : ERROR_STOP
USE TRACER_MOD, ONLY : ITS_A_FULLCHEM_SIM, ITS_AN_AEROSOL_SIM
USE TRACER_MOD, ONLY : XNUMOLAIR
USE TRACERID_MOD, ONLY : IDO3
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: I, J, L
! Local variables
INTEGER :: JLOOP
REAL*8 :: O3_VV
!=================================================================
! GET_O3 begins here!
!=================================================================
IF ( ITS_A_FULLCHEM_SIM() ) THEN
!--------------------
! Coupled simulation
!--------------------
! JLOOP = SMVGEAR 1-D grid box index
JLOOP = JLOP(I,J,L)
! Get O3 from CSPEC [molec/cm3] and convert it to [v/v]
! O3 data will only be defined below the tropopause
IF ( JLOOP > 0 ) THEN
O3_VV = ( CSPEC(JLOOP,IDO3) * 1d6 ) /
& ( AIRDEN(L,I,J) * XNUMOLAIR )
ELSE
O3_VV = 0d0
ENDIF
ELSE IF ( ITS_AN_AEROSOL_SIM() ) THEN
!--------------------
! Offline simulation
!--------------------
! Get O3 [v/v] for this gridbox & month
! O3 data will only be defined below the tropopause
IF ( L <= LLTROP ) THEN
O3_VV = O3m(I,J,L)
ELSE
O3_VV = 0d0
ENDIF
ELSE
!--------------------
! Invalid simulation
!--------------------
CALL ERROR_STOP( 'Invalid NSRCX!', 'GET_OH (sulfate_mod.f)')
ENDIF
! Return to calling program
END FUNCTION GET_O3
!------------------------------------------------------------------------------
SUBROUTINE READ_NONERUP_VOLC 1,7
!
!******************************************************************************
! Subroutine READ_NONERUP_VOLC reads SO2 emissions from non-eruptive
! volcanoes. (rjp, bdf, bmy, 9/19/02, 10/3/05)
!
! NOTES:
! (1 ) Split off from old module routine "sulfate_readyr" (bmy, 9/19/02)
! (2 ) Now references DATA_DIR from "directory_mod.f" (bmy, 7/20/04)
! (3 ) Now read files from "sulfate_sim_200508/" (bmy, 7/28/05)
! (4 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
!******************************************************************************
!
! References to F90 modules
USE BPCH2_MOD, ONLY : GET_RES_EXT, GET_TAU0, READ_BPCH2
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE FILE_MOD, ONLY : IU_FILE, IOERROR
# include "CMN_SIZE" ! Size parameters
! Local variables
INTEGER :: I, IOS, J, K, L
REAL*8 :: EE
CHARACTER(LEN=255) :: FILENAME
!=================================================================
! READ_NONERUP_VOLC begins here!
!=================================================================
! Initialize
K = 1
Env = 0.d0
FILENAME = TRIM( DATA_DIR ) //
& 'sulfate_sim_200508/volcano.con.' // GET_RES_EXT()
!=================================================================
! Read NON-eruptive volcanic SO2 emission (GEIA) into Env.
! Convert Env from [Mg SO2/box/day] to [kg SO2/box/s].
!=================================================================
! Fancy output
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - READ_NONERUP_VOLC: Reading ', a )
! Open file
OPEN( IU_FILE, FILE=TRIM( FILENAME ), STATUS='OLD', IOSTAT=IOS )
IF ( IOS > 0 ) CALL IOERROR( IOS, IU_FILE, 'read_nonerup_volc:1' )
! Read header lines
DO L = 1, 2
READ( IU_FILE, '(a)', IOSTAT=IOS )
IF ( IOS > 0 ) THEN
CALL IOERROR( IOS, IU_FILE, 'read_nonerup_volc:2' )
ENDIF
ENDDO
! Read data values
DO
READ( IU_FILE, '(49x,i4,e11.3,1x,2i4)', IOSTAT=IOS )
& IELVn(k), EE, INV(K), JNV(k)
! Check for EOF
IF ( IOS < 0 ) EXIT
! Trap I/O error
IF ( IOS > 0 ) THEN
CALL IOERROR( IOS, IU_FILE, 'read_nonerup_volc:3' )
ENDIF
! Unit conversion: [Mg SO2/box/day] -> [kg SO2/box/s]
Env(k) = EE * 1000.d0 / ( 24.d0 * 3600.d0 )
! Increment counter
K = K + 1
ENDDO
! Close file
CLOSE( IU_FILE )
! NNVOL = Number of non-eruptive volcanoes
NNVOL = K - 1
! Return to calling program
END SUBROUTINE READ_NONERUP_VOLC
!------------------------------------------------------------------------------
SUBROUTINE READ_ERUP_VOLC 1,7
!
!*****************************************************************************
! Subroutine READ_ERUP_VOLC reads SO2 emissions from eruptive
! volcanoes. (rjp, bdf, bmy, 9/19/02, 10/3/05)
!
! NOTES:
! (1 ) Split off from old module routine "sulfate_readyr" (bmy, 9/19/02)
! (2 ) Now references DATA_DIR from "directory_mod.f" (bmy, 7/20/04)
! (3 ) Now read files from "sulfate_sim_200508/" (bmy, 7/28/05)
! (4 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
!*****************************************************************************
!
! References to F90 modules
USE BPCH2_MOD, ONLY : GET_RES_EXT, GET_TAU0, READ_BPCH2
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE FILE_MOD, ONLY : IU_FILE, IOERROR
# include "CMN_SIZE" ! Size parameters
! Local variables
INTEGER :: I, IOS, IUNIT, J, K, L, M
REAL*8 :: A, B, Fe, X, EE
CHARACTER(LEN=255) :: FILENAME
!==================================================================
! READ_ERUP_VOLC begins here
!==================================================================
! Initialize
K = 1
Eev(:) = 0.d0
FILENAME = TRIM( DATA_DIR ) //
& 'sulfate_sim_200508/volcano.erup.1990.' //
& GET_RES_EXT()
!==================================================================
! Read eruptive volcanic SO2 emission (based on Smithsonian data
! base, SO2 emission and cloud height are a function of VEI.
! Data are over-written if TOMS observations are available.
! Also define a slab with a thickness of 1/3 of the cloud column,
! and SO2 are emitted uniformely within the slab.
!
! Convert Ee from [kton SO2] to [kg SO2/box] and store in Eev.
! ESO2_ev(i,j,l) in [kg SO2/box/s] will be calculated in SRCSO2.
!==================================================================
! Fancy output
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - READ_ERUP_VOLC: Reading ', a )
! Open file
OPEN( IU_FILE, FILE=TRIM( FILENAME ), STATUS='OLD', IOSTAT=IOS )
IF ( IOS > 0 ) CALL IOERROR( IOS, IU_FILE, 'read_erup_volc:1' )
! Read header lines
DO L = 1, 2
READ( IU_FILE, '(a)', IOSTAT=IOS )
IF ( IOS > 0 ) THEN
CALL IOERROR( IOS, IU_FILE, 'read_erup_volc:2' )
ENDIF
ENDDO
! Read data values
DO
READ( IU_FILE, '(47x,3i6,6x,i6,es11.3,1x,2i4)', IOSTAT=IOS )
& IELVe(K), IDAYs(K), IDAYe(K), IHGHT(K),
& Ee, IEV(K), JEV(K)
! Check for EOF
IF ( IOS < 0 ) EXIT
! Trap I/O error
IF ( IOS > 0 ) THEN
CALL IOERROR( IOS, IU_FILE, 'sulfate_readyr:6' )
ENDIF
! Unit conversion: [kton SO2/box/event] -> [kg SO2/box/event]
Eev(k) = Ee * 1.d6
! Increment count
K = K + 1
ENDDO
! Close file
CLOSE( IU_FILE )
! NEVOL = Number of eruptive volcanoes
NEVOL = K - 1
! Return to calling program
END SUBROUTINE READ_ERUP_VOLC
!------------------------------------------------------------------------------
SUBROUTINE READ_ANTHRO_SOx( THISMONTH, NSEASON ) 1,33
!
!******************************************************************************
! Suborutine READ_ANTHRO_SOx reads the anthropogenic SOx from disk,
! and partitions it into anthropogenic SO2 and SO4.
! (rjp, bdf, bmy, 9/20/02, 8/17/06)
!
! Arguments as Input:
! ============================================================================
! (1 ) THISMONTH (INTEGER) : Current month number (1-12)
!
! NOTES:
! (1 ) Now use functions GET_XMID and GET_YMID to compute lon and lat
! centers of grid box (I,J). Now replace DXYP(JREF)*1d4 with routine
! GET_AREA_CM2 of "grid_mod.f". Now use functions GET_MONTH and
! GET_YEAR of time_mod.f". Now call READ_BPCH2 with QUIET=.TRUE.
! (bmy, 3/27/03)
! (2 ) Now references DATA_DIR from "directory_mod.f". Also removed
! reference to CMN, it's not needed. (bmy, 7/20/04)
! (3 ) Now read files from "sulfate_sim_200508/". Now read data for both
! GCAP and GEOS grids (bmy, 8/16/05)
! (4 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
! (5 ) Now references XNUMOL from "tracer_mod.f" (bmy, 10/25/05)
! (6 ) Now computes future SOx emissions (swu, bmy, 5/30/06)
! (7 ) Now can read either EDGAR or GEIA emissions (avd, bmy, 7/14/06)
! (8 ) Now overwrite David Streets' SO2, if necessary (yxw, bmy, 8/14/06)
!******************************************************************************
!
! References to F90 modules
USE BPCH2_MOD, ONLY : GET_NAME_EXT_2D, GET_RES_EXT
USE BPCH2_MOD, ONLY : GET_TAU0, READ_BPCH2
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE EDGAR_MOD, ONLY : GET_EDGAR_ANTH_SO2
USE FUTURE_EMISSIONS_MOD, ONLY : GET_FUTURE_SCALE_SO2ff
USE GRID_MOD, ONLY : GET_XMID, GET_YMID
USE GRID_MOD, ONLY : GET_AREA_CM2
USE LOGICAL_MOD, ONLY : LFUTURE, LEDGARSOx
USE LOGICAL_MOD, ONLY : LSTREETS
USE STREETS_ANTHRO_MOD, ONLY : GET_SE_ASIA_MASK
USE STREETS_ANTHRO_MOD, ONLY : GET_STREETS_ANTHRO
USE TIME_MOD, ONLY : GET_YEAR
USE TRACER_MOD, ONLY : XNUMOL
USE TRACERID_MOD, ONLY : IDTSO2, IDTSO4
USE TRANSFER_MOD, ONLY : TRANSFER_2D
# include "CMN_SIZE" ! Size parameters
# include "CMN_O3" ! FSCALYR
! Arguments
INTEGER, INTENT(IN) :: THISMONTH, NSEASON
! Local variables
INTEGER :: I, J, L, IX, JX, IOS
INTEGER, SAVE :: LASTYEAR = -99
REAL*4 :: E_SOx(IGLOB,JGLOB,2)
REAL*4 :: ARRAY(IGLOB,JGLOB,1)
REAL*8 :: XTAU, Fe, X, Y, AREA_CM2
REAL*8 :: EDG_SO2
CHARACTER(LEN=4) :: SYEAR
CHARACTER(LEN=255) :: FILENAME
!=================================================================
! READ_ANTHRO_SOx begins here!
!=================================================================
IF ( LEDGARSOx ) THEN
!==============================================================
! Use EDGAR SOx emissions
!
! Partition SOx into SO2 and SO4, according to the following
! fractions (cf Chin et al, 2000):
!
! Europe [ 36N-78N, 12.5W-60.0E ]: 5.0% of SOx is SO4
! 95.0% of SOx is SO2
!
! N. America [ 26N-74N, 167.5W-52.5W ]: 1.4% of SOx is SO4
! 98.6% of SOx is SO2
!
! Everywhere else: 3.0% of SOx is SO4
! 97.0% of SOx is SO2
!==============================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, X, Y, EDG_SO2, Fe )
! Loop over latitudes
DO J = 1, JJPAR
! Latitude [degrees]
Y = GET_YMID( J )
! Loop over longitudes
DO I = 1, IIPAR
! Longitude [degrees]
X = GET_XMID( I )
! Get EDGAR SO2 emissions [kg/s]
! NOTE: Future emissions are already applied!
EDG_SO2 = GET_EDGAR_ANTH_SO2( I, J, KG_S=.TRUE. )
! If we are using David Streets' emissions ...
IF ( LSTREETS ) THEN
! If we are over the SE Asia region ...
IF ( GET_SE_ASIA_MASK( I, J ) > 0d0 ) THEN
! Overwrite EDGAR SO2 w/ David Streets' [kg SO2/s]
EDG_SO2 = GET_STREETS_ANTHRO( I, J,
& 26, KG_S=.TRUE. )
ENDIF
ENDIF
! Compute SO4/SOx fraction for EUROPE
IF ( ( X >= -12.5 .and. X <= 60.0 ) .and.
& ( Y >= 36.0 .and. Y <= 78.0 ) ) THEN
Fe = 0.05d0
! Compute SO4/SOx fraction for NORTH AMERICA
ELSE IF ( ( X >= -167.5 .and. X <= -52.5 ) .and.
& ( Y >= 26.0 .and. Y <= 74.0 ) ) THEN
Fe = 0.014d0
! Compute SO4/SOx fraction for EVERYWHERE ELSE
ELSE
Fe = 0.03d0
ENDIF
! Store SO2 emission [kg SO2/s]
ESO2_an(I,J,1) = EDG_SO2
ESO4_an(I,J,2) = 0d0
! Compute SO4 from SO2 [kg SO4/s]
ESO4_an(I,J,1) = EDG_SO2 * Fe / ( 1.d0 - Fe )
ESO4_an(I,J,2) = 0d0
ENDDO
ENDDO
!$OMP END PARALLEL DO
ELSE
!==============================================================
! Use GEIA SOx emissions
!==============================================================
! Define filename
FILENAME = TRIM( DATA_DIR ) //
& 'fossil_200104/merge_nobiofuels.' //
& GET_NAME_EXT_2D() // '.' // GET_RES_EXT()
! Echo output
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - READ_ANTHRO_SOx: Reading ', a )
! Pick the right TAU value for the given season
! Seasons: 1=DJF, 2=MAM, 3=JJA, 4=SON
SELECT CASE ( NSEASON )
CASE ( 1 )
XTAU = -744d0
CASE ( 2 )
XTAU = 1416d0
CASE ( 3 )
XTAU = 3624d0
CASE ( 4 )
XTAU = 5832d0
END SELECT
! Read anthropogenic SOx [molec/cm2/s]
CALL READ_BPCH2( FILENAME, 'ANTHSRCE', 27,
& XTAU, IGLOB, JGLOB,
& 2, E_SOx, QUIET=.TRUE. )
!=================================================================
! Read in yearly SO2 scale factors here
! (For now we only have 1998, deal w/ other years later)
!=================================================================
IF ( LASTYEAR < 0 ) THEN
! put in SOX scale year here (hardwired to 1998 for now)
SYEAR = '1998'
FILENAME = TRIM( DATA_DIR ) //
& 'sulfate_sim_200508/scalefoss.SOx.' //
& GET_RES_EXT() // '.' // SYEAR
! Echo output
WRITE( 6, 100 ) TRIM( FILENAME )
! Get TAU value (use Jan 1, 1998 for scale factors)
XTAU = GET_TAU0( 1, 1, 1998 )
! Read anthropogenic SOx [molec/cm2/s]
CALL READ_BPCH2( FILENAME, 'SCALFOSS', 3,
& XTAU, IGLOB, JGLOB,
& 1, ARRAY, QUIET=.TRUE. )
! Cast from REAL*4 to REAL*8
CALL TRANSFER_2D( ARRAY(:,:,1), SOx_SCALE )
! Reset LASTYEAR
LASTYEAR = GET_YEAR()
ENDIF
!==============================================================
! Partition SOx into SO2 and SO4, according to the following
! fractions (cf Chin et al, 2000):
!
! Europe [ 36N-78N, 12.5W-60.0E ]: 5.0% of SOx is SO4
! 95.0% of SOx is SO2
!
! N. America [ 26N-74N, 167.5W-52.5W ]: 1.4% of SOx is SO4
! 98.6% of SOx is SO2
!
! Everywhere else: 3.0% of SOx is SO4
! 97.0% of SOx is SO2
!==============================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, L, AREA_CM2, Y, X, Fe )
DO L = 1, 2
DO J = 1, JJPAR
! Grid box surface area [cm2]
AREA_CM2 = GET_AREA_CM2( J )
! Latitude [degrees]
Y = GET_YMID( J )
DO I = 1, IIPAR
! Longitude [degrees]
X = GET_XMID( I )
! First scale SOx to the given fossil fuel year
E_SOx(I,J,L) = E_SOx(I,J,L) * SOx_SCALE(I,J)
! Compute future SOx emissions (if necessary)
IF ( LFUTURE ) THEN
E_SOx(I,J,L) = E_SOx(I,J,L) *
& GET_FUTURE_SCALE_SO2ff( I, J )
ENDIF
! Compute SO4/SOx fraction for EUROPE
IF ( ( X >= -12.5 .and. X <= 60.0 ) .and.
& ( Y >= 36.0 .and. Y <= 78.0 ) ) THEN
Fe = 0.05d0
! Compute SO4/SOx fraction for NORTH AMERICA
ELSE IF ( ( X >= -167.5 .and. X <= -52.5 ) .and.
& ( Y >= 26.0 .and. Y <= 74.0 ) ) THEN
Fe = 0.014d0
! Compute SO4/SOx fraction for EVERYWHERE ELSE
ELSE
Fe = 0.03d0
ENDIF
! Compute SO2 (tracer #2) from SOx
! Convert from [molec SOx/cm2/s] to [kg SO2/box/s]
ESO2_an(I,J,L) = E_SOx(I,J,L) * ( 1.d0 - Fe ) *
& AREA_CM2 / XNUMOL(IDTSO2)
! If we are using David Streets' emissions
IF ( LSTREETS ) THEN
! If we are over the SE Asia region
IF ( GET_SE_ASIA_MASK( I, J ) > 0d0 ) THEN
! Overwrite GEIA SO2 w/ David Streets' SO2 [kg SO2/s]
ESO2_an(I,J,1) = GET_STREETS_ANTHRO( I, J, 26,
& KG_S=.TRUE. )
! Zero 2nd level of emissions
ESO2_an(I,J,2) = 0d0
ENDIF
ENDIF
! Compute SO4 (tracer #3) from SOx
! Convert from [molec SOx/cm2/s] to [kg SO4/box/s]
ESO4_an(I,J,L) = E_SOx(I,J,L) * Fe *
& AREA_CM2 / XNUMOL(IDTSO4)
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
ENDIF
! Return to calling program
END SUBROUTINE READ_ANTHRO_SOx
!------------------------------------------------------------------------------
SUBROUTINE READ_OCEAN_DMS( THISMONTH ) 1,8
!
!*****************************************************************************
! Subroutine READ_OCEAN_DMS reads seawater concentrations of DMS (nmol/L).
! (rjp, bdf, bmy, 9/20/02, 10/3/05)
!
! Arguments as input:
! ===========================================================================
! (1 ) THISMONTH (INTEGER) : Current month number (1-12)
!
! NOTES:
! (1 ) Extracted from old module routine SULFATE_READMON (bmy, 9/18/02)
! (2 ) Now call READ_BPCH2 with QUIET=.TRUE. (bmy, 3/27/03)
! (3 ) Now references DATA_DIR from "directory_mod.f" (bmy, 7/20/04)
! (4 ) Now read files from "sulfate_sim_200508/". Now read data for both
! GCAP and GEOS grids (bmy, 8/16/05)
! (5 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
!*****************************************************************************
!
! References to F90 modules
USE BPCH2_MOD, ONLY : GET_NAME_EXT_2D, GET_RES_EXT
USE BPCH2_MOD, ONLY : GET_TAU0, READ_BPCH2
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE TRANSFER_MOD, ONLY : TRANSFER_2D
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: THISMONTH
! Local variables
REAL*4 :: ARRAY(IGLOB,JGLOB,1)
REAL*8 :: XTAU
CHARACTER(LEN=255) :: FILENAME
!==================================================================
! READ_OCEAN_DMS begins here!
!==================================================================
! File name
FILENAME = TRIM( DATA_DIR ) //
& 'sulfate_sim_200508/DMS_seawater.' //
& GET_NAME_EXT_2D() // '.' // GET_RES_EXT()
! Echo output
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - READ_OCEAN_DMS: Reading ', a )
! TAU value (use generic year 1985)
XTAU = GET_TAU0( THISMONTH, 1, 1985 )
! Read ocean DMS [nmol/L]
CALL READ_BPCH2( FILENAME, 'BIOGSRCE', 25,
& XTAU, IGLOB, JGLOB,
& 1, ARRAY(:,:,1), QUIET=.TRUE. )
! Cast from REAL*4 to REAL*8
CALL TRANSFER_2D( ARRAY(:,:,1), DMSo )
! Return to calling program
END SUBROUTINE READ_OCEAN_DMS
!------------------------------------------------------------------------------
SUBROUTINE READ_SST( THISMONTH, THISYEAR ) 1,11
!
!*****************************************************************************
! Subroutine READ_SST reads monthly mean sea surface temperatures.
! (rjp, bdf, bmy, 9/18/02, 11/17/05)
!
! Arguments as input:
! ===========================================================================
! (1 ) THISMONTH (INTEGER) : Current month number (1-12)
! (2 ) THISYEAR (INTEGER) : Current 4-digit year
!
! NOTES:
! (1 ) Extracted from old module routine SULFATE_READMON (bmy, 9/18/02)
! (2 ) Now call READ_BPCH2 with QUIET=.TRUE. (bmy, 3/27/03)
! (3 ) Now references DATA_DIR from "directory_mod.f" (bmy, 7/20/04)
! (4 ) Now use interannual SST data from NOAA if present; otherwise use
! climatological SST data. Now read data for both GCAP and GEOS
! grids (bmy, 8/16/05)
! (5 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
! (6 ) Now read int'annual SST data on the GEOS 1x1 grid (bmy, 11/17/05)
!*****************************************************************************
!
! References to F90 modules
USE BPCH2_MOD, ONLY : GET_NAME_EXT_2D, GET_RES_EXT
USE BPCH2_MOD, ONLY : GET_TAU0, READ_BPCH2
USE DIRECTORY_MOD, ONLY : DATA_DIR, DATA_DIR_1x1
USE REGRID_1x1_MOD, ONLY : DO_REGRID_1x1
USE TRANSFER_MOD, ONLY : TRANSFER_2D
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: THISMONTH, THISYEAR
! Local variables
REAL*4 :: ARRAY(IGLOB,JGLOB,1)
REAL*4 :: ARRAY1(I1x1,J1x1,1)
REAL*8 :: XTAU
CHARACTER(LEN=4) :: SYEAR
CHARACTER(LEN=255) :: FILENAME
!==================================================================
! READ_SST begins here!
!==================================================================
IF ( THISYEAR >= 1985 .and. THISYEAR <= 2004 ) THEN
!------------------------------------
! Use interannual SST data from NOAA
! Data exists for 1985 - 2004,
! Add other years as necessary
!------------------------------------
! Make a string for THISYEAR
WRITE( SYEAR, '(i4)' ) THISYEAR
! File name
FILENAME = TRIM( DATA_DIR_1x1 ) //
& 'SST_200508/SST.geos.1x1.' // SYEAR
! Echo output
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - READ_SST: Reading ', a )
! TAU value (use this year)
XTAU = GET_TAU0( THISMONTH, 1, THISYEAR )
! Read sea surface temperature [K]
CALL READ_BPCH2( FILENAME, 'GMAO-2D', 69,
& XTAU, I1x1, J1x1,
& 1, ARRAY1(:,:,1), QUIET=.TRUE. )
! Regrid from 1x1 and cast to REAL*8
CALL DO_REGRID_1x1( 'K', ARRAY1, SSTEMP )
ELSE
!-------------------------------
! Use climatological SST data
!-------------------------------
! File name
FILENAME = TRIM( DATA_DIR ) //
& 'sulfate_sim_200508/SST.' // GET_NAME_EXT_2D() //
& '.' // GET_RES_EXT()
! Echo output
WRITE( 6, 100 ) TRIM( FILENAME )
! TAU value (use generic year 1985)
XTAU = GET_TAU0( THISMONTH, 1, 1985 )
! Read sea surface temperature [K]
CALL READ_BPCH2( FILENAME, 'DAO-FLDS', 5,
& XTAU, IGLOB, JGLOB,
& 1, ARRAY(:,:,1), QUIET=.TRUE. )
! Cast from REAL*4 to REAL*8
CALL TRANSFER_2D( ARRAY(:,:,1), SSTEMP )
ENDIF
! Return to calling program
END SUBROUTINE READ_SST
!------------------------------------------------------------------------------
SUBROUTINE READ_BIOMASS_SO2( THISMONTH ) 1,17
!
!******************************************************************************
! Subroutine READ_BIOMASS_SO2 reads monthly mean biomass burning
! emissions for SO2. SOx is read in, and converted to SO2.
! (rjp, bdf, bmy, 1/16/03, 10/3/06)
!
! Arguments as input:
! ===========================================================================
! (1 ) THISMONTH (INTEGER) : Current month number (1-12)
!
! NOTES:
! (1 ) Extracted from old module routine SULFATE_READMON (bmy, 9/18/02)
! (2 ) Modified molar ratio of biomass burning SO2 per CO. Added SO2
! emission from biofuel burning. (rjp, bmy, 1/16/03)
! (3 ) Now replace DXYP(J+J0)*1d4 with routine GET_AREA_CM2 of "grid_mod.f"
! Now replace MONTH with the argument THISMONTH. Now call READ_BPCH2
! with QUIET=.TRUE. (bmy, 3/27/03)
! (4 ) Now references DATA_DIR from "directory_mod.f". Also removed
! references to CMN and CMN_SETUP. (bmy, 7/20/04)
! (5 ) Now can read either seasonal or interannual biomass burning emissions.
! Now references routines from both "logical_mod.f" and "time_mod.f".
! Now reads SO2 biomass emissions directly rather than computing
! it by mole fraction from CO. (rjp, bmy, 1/11/05)
! (6 ) Now read data for both GCAP and GEOS grids (bmy, 8/16/05)
! (7 ) Now make sure all USE statements are USE, ONLY (bmy, 10/3/05)
! (8 ) Now computes future biomass emissions, if necessary (swu, bmy, 5/30/06)
! (9 ) Now only read the biofuel, we have moved the biomass-reading code
! to "gc_biomass_mod.f" for compatibility with GFED2 biomass emissions
! (bmy, 9/27/06)
! (10) Now prevent seg fault if BIOMASS emissions are turned off.
! (bmy, 10/3/06)
!******************************************************************************
!
! References to F90 modules
USE BIOMASS_MOD, ONLY : BIOMASS, IDBSO2
USE BPCH2_MOD, ONLY : GET_NAME_EXT_2D, GET_RES_EXT
USE BPCH2_MOD, ONLY : GET_TAU0, READ_BPCH2
USE DIRECTORY_MOD, ONLY : DATA_DIR
USE FUTURE_EMISSIONS_MOD, ONLY : GET_FUTURE_SCALE_SO2bf
USE GRID_MOD, ONLY : GET_AREA_CM2
USE LOGICAL_MOD, ONLY : LBIOMASS, LFUTURE
USE TIME_MOD, ONLY : ITS_A_LEAPYEAR
USE TRACER_MOD, ONLY : XNUMOL
USE TRACERID_MOD, ONLY : IDTSO2
USE TRANSFER_MOD, ONLY : TRANSFER_2D
# include "CMN_SIZE" ! Size parameters
! Arguments
INTEGER, INTENT(IN) :: THISMONTH
! Local variables
INTEGER :: I, J, THISYEAR
REAL*4 :: ARRAY(IGLOB,JGLOB,1)
REAL*8 :: BIOCO(IIPAR,JJPAR)
REAL*8 :: CONV, XTAU
CHARACTER(LEN=4 ) :: CYEAR
CHARACTER(LEN=255) :: FILENAME
! Days per month
REAL*8 :: NMDAY(12) = (/ 31d0, 28d0, 31d0, 30d0, 31d0, 30d0,
& 31d0, 31d0, 30d0, 31d0, 30d0, 31d0/)
!=================================================================
! READ_BIOMASS_SO2 begins here!
!=================================================================
!=================================================================
! Compute biofuel SO2 from biofuel CO. Use a molar
! ratio of 0.0015 moles SO2/mole CO from biofuel burning.
! (Table 2, [Andreae and Merlet, 2001])
!=================================================================
! File name for biofuel burning file
FILENAME = TRIM( DATA_DIR ) //
& 'biofuel_200202/biofuel.' // GET_NAME_EXT_2D() //
& '.' // GET_RES_EXT()
! Echo info
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( ' - READ_BIOMASS_SO2: Reading ', a )
! Get TAU0 value (use generic year 1985)
XTAU = GET_TAU0( 1, 1, 1985 )
! Read Biofuel burning of CO [kg/yr]
CALL READ_BPCH2( FILENAME, 'BIOFSRCE', 4,
& XTAU, IGLOB, JGLOB,
& 1, ARRAY(:,:,1), QUIET=.TRUE. )
! Cast from REAL*4 to REAL*8
CALL TRANSFER_2D( ARRAY(:,:,1), BIOCO )
!=================================================================
! Unit conversion to [kg SO2/s]
!=================================================================
!$OMP PARALLEL DO
!$OMP+DEFAULT( SHARED )
!$OMP+PRIVATE( I, J, CONV )
! Loop over longitudes
DO J = 1, JJPAR
! Conversion factor for [cm2 * kg/molec]
CONV = GET_AREA_CM2( J ) / XNUMOL(IDTSO2)
! Loop over latitudes
DO I = 1, IIPAR
! Convert biomass SO2 from [molec SO2/cm2/s] -> [kg SO2/s]
! NOTE: Future scale has already been applied by this point
IF ( LBIOMASS ) THEN
ESO2_bb(I,J) = BIOMASS(I,J,IDBSO2) * CONV
ELSE
ESO2_bb(I,J) = 0d0
ENDIF
! Convert biofuel SO2 from [kg CO/yr] to [kg SO2/s]
ESO2_bf(I,J) = ( BIOCO(I,J) * 64d-3 * 0.0015d0 /
& ( 28d-3 * 86400.d0 * 365.25d0 ) )
! Apply future emissions to biofuel SO2, if necessary
IF ( LFUTURE ) THEN
ESO2_bf(I,J) = ESO2_bf(I,J) *
& GET_FUTURE_SCALE_SO2bf( I, J )
ENDIF
ENDDO
ENDDO
!$OMP END PARALLEL DO
! Return to calling program
END SUBROUTINE READ_BIOMASS_SO2
!------------------------------------------------------------------------------
SUBROUTINE READ_AIRCRAFT_SO2( THISMONTH ) 1,8
!
!******************************************************************************
! Subroutine READ_AIRCRAFT_SO2 reads monthly mean aircraft fuel emissions
! and converts them to SO2 emissions.