C*********************************************************************** C Portions of Models-3/CMAQ software were developed or based on * C information from various groups: Federal Government employees, * C contractors working on a United States Government contract, and * C non-Federal sources (including research institutions). These * C research institutions have given the Government permission to * C use, prepare derivative works, and distribute copies of their * C work in Models-3/CMAQ to the public and to permit others to do * C so. EPA therefore grants similar permissions for use of the * C Models-3/CMAQ software, but users are requested to provide copies * C of derivative works to the Government without restrictions as to * C use by others. Users are responsible for acquiring their own * C copies of commercial software associated with Models-3/CMAQ and * C for complying with vendor requirements. Software copyrights by * C the MCNC Environmental Modeling Center are used with their * C permissions subject to the above restrictions. * C*********************************************************************** C RCS file, release, date & time of last delta, author, state, [locker] C $Header: /project/work/rep/CCTM/src/aero/aero4/aero_driver.F,v 1.6 2006/05/30 14:49:31 sjr Exp $ C what(1) key, module and SID; SCCS file; date and time of last delta: C %W% %P% %G% %U% C >>> 08/04/2000 Changes necessary to be able to read and process C two different types of emissions files. C the first type is the existing opperational PM2.5 & PM10 unspeciated C file. The new file format has speciated emissions. C >>> This version uses the FORTRAN 90 feature for runtime memory C allocation. C 1/12/99 David Wong at LM: C -- introduce new variable MY_NUMBLKS (eliminate NUMBLKS) C -- re-calculate NOXYZ accordingly C FSB Updated for inclusion of surface area / second moment C 25 Sep 00 (yoj) various bug fixes, cleanup to coding standards C Jeff - Dec 00 - move CGRID_MAP into f90 module C FSB/Jeff - May 01 - optional emissions processing C Jerry Gipson - Jun 01 - added SOA linkages for saprc99 C Bill Hutzell - Jun 01 - simplified CBLK mapping C Jerry Gipson - Jun 03 - modified for new soa treatment C Jerry Gipson - Aug 03 - removed SOA prod form alkenes & added C emission adjustment factors for ALK & TOL ( RADM2 & SAPRC99 only) C Shawn Roselle - Jan 04 C - removed SOA from transported aerosol surface area C - fixed bug in calculation of wet parameters. Previously, DRY aerosol C parameters were being written to the AERDIAG files and mislabeled C as WET. C Prakash Bhave - May 04 C - changed AERODIAG species (added RH; removed M0 & M2dry) C J.Young 31 Jan 05: dyn alloc - establish both horizontal & vertical C domain specifications in one module c Uma Shankar and Prakash Bhave - Jun 05 c - added code to handle the following species: ANAI, ANAJ, ANAK, ACLI, c ACLJ, ACLK, ASO4K, AH2OK, ANO3K, and HCL; removed code for ASEAS c - removed obsolete MW variables C Prakash Bhave - Jul 05 - added PM25 mass-fraction calculations C J.Young 16 Feb 06: trap fractional humidity above 0.005 C Prakash Bhave - Apr 06 - added GAMMA_N2O5 to the AEROPROC call vector C and the aerosol diagnostic file C Prakash Bhave - May 06 - changed units of DG variables from m to um in C the aerosol diagnostic file as suggested by Dr. Bill Hutzell C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE AERO ( JDATE, JTIME, TSTEP ) USE CGRID_DEFN ! inherits GRID_CONF and CGRID_SPCS USE VIS_DEFN USE AERO_INFO_AE4 ! replaces aero include files IMPLICIT NONE C *** includes: ! INCLUDE SUBST_HGRD_ID ! horizontal dimensioning parameters ! INCLUDE SUBST_VGRD_ID ! vertical dimensioning parameters INCLUDE SUBST_GC_SPC ! gas chemistry species table INCLUDE SUBST_GC_EMIS ! gas chem emis surrogate names and map ! table INCLUDE SUBST_AE_SPC ! aerosol species table ! INCLUDE SUBST_AE_EMIS ! aerosol emis surrogate names and map ! table INCLUDE SUBST_NR_SPC ! non-reactive species table INCLUDE SUBST_GC_G2AE ! gas chem aerosol species and map table INCLUDE SUBST_NR_N2AE ! non-react aerosol species and map table INCLUDE SUBST_RXCMMN ! to get mech name ! INCLUDE SUBST_CONST ! constants INCLUDE SUBST_IOPARMS ! I/O parameters definitions INCLUDE SUBST_IOFDESC ! file header data structure #include SUBST_IODECL # I/O definitions and declarations INCLUDE SUBST_FILES_ID ! file name parameters ! INCLUDE SUBST_BLKPRM ! sets BLKSIZE AND MXBLKS ! INCLUDE SUBST_COORD_ID ! coord. and domain definitions ! (req IOPARMS) ! INCLUDE SUBST_PACTL_ID ! process analysis C *** arguments: INTEGER JDATE ! Current model date , coded YYYYDDD INTEGER JTIME ! Current model time , coded HHMMSS INTEGER TSTEP( 2 ) ! time step vector (HHMMSS) ! TSTEP(1) = local output step ! TSTEP(2) = sciproc sync. step (chem) REAL, PARAMETER :: CONMIN = 1.0E-30 ! concentration lower limit INTEGER, SAVE :: LOGDEV ! unit number for the log file C *** local variables: CHARACTER( 16 ), SAVE :: PNAME = 'AERO_DRIVER' CHARACTER( 16 ) :: VNAME ! variable name CHARACTER( 96 ) :: XMSG = ' ' INTEGER MDATE, MTIME, MSTEP ! julian date, time and ! timestep in sec INTEGER C, R, L, V, N ! loop counters INTEGER SPC ! species loop counter INTEGER STRT, FINI ! loop induction variables INTEGER ALLOCSTAT ! memory allocation status INTEGER LAYER ! model layer LOGICAL LERROR ! Error flag C *** External Functions not previously declared in IODECL3.EXT: INTEGER, EXTERNAL :: SECSDIFF, SEC2TIME, TIME2SEC, INDEX1, SETUP_LOGDEV LOGICAL, EXTERNAL :: ENVYN ! get environment variable as boolean C *** Grid description REAL DX1 ! Cell x-dimension REAL DX2 ! Cell y-dimension C *** Variable to set time step for writing visibility file INTEGER, SAVE :: WSTEP = 0 ! local write counter LOGICAL, SAVE :: WRITETIME = .FALSE. ! local write flag C *** meteorological variables C *** Reciprocal (air density X Jacobian, where Jacobian = sq. root of the C determinant of the metric tensor) at midlayer -inverted after read REAL PRES ( NCOLS,NROWS,NLAYS ) ! Atmospheric pressure [ Pa ] REAL TA ( NCOLS,NROWS,NLAYS ) ! Air temperature [ K ] REAL DENS ( NCOLS,NROWS,NLAYS ) ! Air density [ kg/m**-3 ] REAL QV ( NCOLS,NROWS,NLAYS ) ! Water vapor mixing ratio [ kg/kg ] C *** variables computed and output but not carried in CGRID C *** visibility variables ! INTEGER, PARAMETER :: N_AE_VIS_SPC = 4 C visual range in deciview (Mie) ! INTEGER, PARAMETER :: IDCVW1 = 1 C extinction [ 1/km ] (Mie) ! INTEGER, PARAMETER :: IBEXT1 = 2 INTEGER, PARAMETER :: IBEXT1 = 1 C visual range in deciview (Reconstructed) ! INTEGER, PARAMETER :: IDCVW2 = 3 C extinction [ 1/km ] (Reconstructed) ! INTEGER, PARAMETER :: IBEXT2 = 4 INTEGER, PARAMETER :: IBEXT2 = 2 ! REAL VIS_SPC( NCOLS,NROWS,N_AE_VIS_SPC ) ! Visual range information C *** aerosol distribution variables ! REAL DIAM_SPC( NCOLS,NROWS,NLAYS,19 ) ! aerosol size distribution variables C *** grid variables C *** information about blocks C *** pointers to gas (vapor) phase species and production rates in CGRID INTEGER, SAVE :: LSULF, LSULFP, LHNO3, LNH3, LHCL, LN2O5 INTEGER, SAVE :: LTOLAER, LXYLAER, LCSLAER INTEGER, SAVE :: LALKAER INTEGER, SAVE :: LOLIAER, LTERPAER, LTERP C *** meteorological information: REAL BLKPRS ! Air pressure in [ Pa ] REAL BLKTA ! Air temperature [ K ] REAL BLKDENS ! Air density [ kg m^-3 ] REAL BLKDENS1 ! Reciprocal of air density REAL BLKQA ! Water vapor mixing ratio REAL BLKESAT ! Saturation water vapor pressure [Pa] REAL BLKEVAP ! Ambient water vapor pressure [Pa] REAL BLKRH ! Fractional relative humidity C *** chemical production rates: [ ug / m**3 s ] REAL SO4RATE ! sulfate gas-phase production rate C *** new information f or secondary organic aerosols INTEGER VV ! loop index for vapors INTEGER VVLOC ! generic index INTEGER, PARAMETER :: NPSPCS = 6 INTEGER, PARAMETER :: NCVAP = 10 CHARACTER( 16 ), SAVE :: VAPNAMES ( NCVAP ) REAL VAPORS( NCVAP ) ! condensible secondary organic vapors INTEGER, SAVE :: LOCVAP( NCVAP ) REAL ORGPROD( NPSPCS ) c..The following paramters derived form annual emission inv estimates c..They represent the fraction of lumped compunds that are SOA precursors REAL, PARAMETER :: SPALK = 0.57 ! Frac of SAPRC99 ALK5 producing SOA REAL, PARAMETER :: SPTOL = 0.93 ! Frac of SAPRC99 ARO1 producing SOA REAL, PARAMETER :: RDALK = 0.56 ! Frac of RADM2 HC8 producing SOA REAL, PARAMETER :: RDTOL = 0.94 ! Frac of RADM2 TOL producing SOA C *** Current implementation ! ORGPROD(1) -> "long" alkanes ( alk5 in saprc99) ! ORGPROD(2) -> internal alkenes ( cyclohexene. ole2 in saprc99 ) ! ORGPROD(3) -> aromatics like xylene (aro2 in saprc99) ! ORGPROD(4) -> aromatics like cresol (cres in saprc99) ! ORGPROD(5) -> aromatics like toluene (aro1 in saprc99) ! ORGPROD(6) -> monoterpenes (trp1 in saprc99) C *** atmospheric properties REAL XLM ! atmospheric mean free path [ m ] REAL AMU ! atmospheric dynamic viscosity [ kg/m s ] C *** modal diameters [ m ] REAL DGATK ! Aitken mode geometric mean diameter [ m ] REAL DGACC ! accumulation geometric mean diameter [ m ] REAL DGCOR ! coarse mode geometric mean diameter [ m ] C *** log of modal geometric standard deviation REAL XXLSGAT ! Aitken mode REAL XXLSGAC ! accumulation mode C *** aerosol properties: C *** modal mass concentrations [ ug m**3 ] REAL PMASSAT ! mass concentration in Aitken mode REAL PMASSAC ! mass concentration in accumulation mode REAL PMASSCO ! mass concentration in coarse mode C *** average modal particle densities [ kg/m**3 ] REAL PDENSAT ! average particle density in Aitken mode REAL PDENSAC ! average particle density in ! accumulation mode REAL PDENSCO ! average particle density in coarse mode C *** N2O5 heterogeneous reaction probability [ ] REAL GAMMA_N2O5 ! computed in SUBROUTINE EQL3 C *** mass fraction of each mode less than 2.5um aerodynamic diameter REAL PM25AT ! fine fraction of Aitken mode REAL PM25AC ! fine fraction of accumulation mode REAL PM25CO ! fine fraction of coarse mode C *** visual range information REAL BLKDCV1 ! block deciview (Mie) REAL BLKEXT1 ! block extinction [ km**-1 ] (Mie) REAL BLKDCV2 ! block deciview (Reconstructed) REAL BLKEXT2 ! block extinction [ km**-1 ] (Reconstructed) C *** molecular weights of gas-phase species REAL, SAVE :: MWH2SO4, MWHNO3, MWNH3, MWHCL, MWN2O5 C *** conversion factors for unit conversions between ppm and ugm**-3 REAL, SAVE :: H2SO4CONV, HNO3CONV, NH3CONV, HCLCONV, N2O5CONV REAL, SAVE :: H2SO4CONV1, HNO3CONV1, NH3CONV1, HCLCONV1, N2O5CONV1 C FSB PM emissions now done in vertical diffusion INTEGER GXOFF, GYOFF ! global origin offset from file C for INTERPX INTEGER, SAVE :: STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3 C *** other internal aerosol variables INTEGER IND ! index to be used with INDEX1 C *** synchronization time step [ s ] REAL DT C *** variables to set up for "dry transport " REAL M3_WET, M3_DRY ! third moment with and without water REAL M2_WET, M2_DRY ! second moment with and without water C flag to include water in the 3rd moment calculation LOGICAL, PARAMETER :: M3_WET_FLAG = .FALSE. C FSB the following variables are defined and set in AERO_INFO_AE4 C 11/09/2001 ! REAL, PARAMETER :: TWO3 = 2.0 / 3.0 ! 2/3 ! REAL, PARAMETER :: F6DPI = 6.0 / PI ! 6/PI ! REAL, PARAMETER :: F6DPIM9 = 1.0E-9 * F6DPI ! 1.0e-9 * 6/PI ! REAL, PARAMETER :: RHOH2O = 1.0E+3 ! bulk density of aerosol water ! REAL, PARAMETER :: H2OFAC = F6DPIM9 / RHOH2O C *** variables aerosol diagnostic file flag INTEGER STATUS ! ENV... status CHARACTER( 80 ) :: VARDESC ! environment variable description C *** environment variable for AERDIAG file ! CHARACTER( 16 ), SAVE :: CTM_AERDIAG = 'CTM_AERDIAG' C *** flag for AERDIAG file [F], default ! LOGICAL, SAVE :: AERDIAG C *** number of species (gas + aerosol) includes 2nd,3rd Moments & 4 gases INTEGER, PARAMETER :: NSPCSDA = N_AE_SPC + 14 ! N2O5 added - FSB ! HCl added - US C *** main array of variables by species REAL CBLK( NSPCSDA ) C *** map of aerosol species INTEGER, SAVE :: N_AE_MAP INTEGER, SAVE :: CBLK_MAP( N_AE_SPCD ) = 0 LOGICAL, SAVE :: FIRSTIME = .TRUE. C *** ratio of molecular weights of water vapor to dry air = 0.622015 REAL, PARAMETER :: EPSWATER = MWWAT / MWAIR C *** mechanism name CHARACTER( 16 ), SAVE :: MECH C *** Statement Function ************** REAL ESATL ! arithmetic statement function for vapor pressure [Pa] REAL TT C *** Coefficients for the equation, ESATL defining saturation vapor pressure REAL, PARAMETER :: AL = 610.94 REAL, PARAMETER :: BL = 17.625 REAL, PARAMETER :: CL = 243.04 C *** values of AL, BL, and CL are from: C Alduchov and Eskridge, "Improved Magnus Form Approximations of C Saturation Vapor Pressure," C Jour. of Applied Meteorology, vol. 35, C pp 601-609, April, 1996. ESATL( TT ) = AL * EXP( BL * ( TT - 273.15 ) / ( TT - 273.15 + CL ) ) C *** End Statement Function ******** C ------------------ begin body of AERO_DRIVER ------------------------- IF ( FIRSTIME ) THEN FIRSTIME = .FALSE. ! LOGDEV = INIT3() LOGDEV = SETUP_LOGDEV() C*** Make sure an ae4 version of the mechanism is being used IF ( INDEX ( MECHNAME, 'AE4' ) .LE. 0 ) THEN XMSG = 'AERO4 requires an AE4 version of chemical mechanism' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF C *** Get aerosol diagnostic file flag. ! AERDIAG = .FALSE. ! default ! VARDESC = 'Flag for writing the aerosol diagnostic file' ! AERDIAG = ENVYN( CTM_AERDIAG, VARDESC, AERDIAG, STATUS ) ! IF ( STATUS .NE. 0 ) WRITE( LOGDEV, '(5X, A)' ) VARDESC ! IF ( STATUS .EQ. 1 ) THEN ! XMSG = 'Environment variable improperly formatted' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 ) ! ELSE IF ( STATUS .EQ. -1 ) THEN ! XMSG = ! & 'Environment variable set, but empty ... Using default:' ! WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME ! ELSE IF ( STATUS .EQ. -2 ) THEN ! XMSG = 'Environment variable not set ... Using default:' ! WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME ! END IF C *** Set up names and indices. C *** Initialize species indices for CBLK. C *** Get CGRID offsets. CALL CGRID_MAP ( NSPCSD, GC_STRT, AE_STRT, NR_STRT, TR_STRT ) C *** Determine CGRID species map from AE_SPC.EXT. V = 0 VNAME = 'ASO4J' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VSO4AJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ASO4I' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VSO4AI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANH4J' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNH4AJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANH4I' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNH4AI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANO3J' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNO3AJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANO3I' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNO3AI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AORGAJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VORGAJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AORGAI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VORGAI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AORGPAJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VORGPAJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AORGPAI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VORGPAI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AORGBJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VORGBAJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AORGBI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VORGBAI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AECJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VECJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AECI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VECI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'A25J' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VP25AJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'A25I' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VP25AI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ACORS' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VANTHA = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ASOIL' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VSOILA = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'NUMATKN' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VAT0 = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'NUMACC' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VAC0 = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'NUMCOR' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VCOR0 = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'SRFATKN' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VSURFAT = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'SRFACC' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VSURFAC = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AH2OJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VH2OAJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AH2OI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VH2OAI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANAJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNAJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANAI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNAI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ACLJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VCLJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ACLI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VCLI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANAK' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNAK = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ACLK' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VCLK = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ASO4K' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VSO4K = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'ANO3K' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VNO3K = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'AH2OK' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VH2OK = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF C define additional AE species for Concentration of HPLUS in Aitken and C Accumulation modes VNAME = 'HPLUSJ' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VHPLUSJ = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) CALL M3MESG ( XMSG ) VHPLUSJ = 0 END IF VNAME = 'HPLUSI' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 VHPLUSI = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) CALL M3MESG ( XMSG ) VHPLUSI = 0 END IF N_AE_MAP = V C additional AE species for mean geometric diameters in Aitken and C accumulation modes VNAME = 'DGATKN_DRY' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 DGATKN_DRY = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) CALL M3MESG ( XMSG ) DGATKN_DRY = 0 END IF VNAME = 'DGACC_DRY' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 DGACC_DRY = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) CALL M3MESG ( XMSG ) DGACC_DRY = 0 END IF VNAME = 'DGATKN_WET' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 DGATKN_WET = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) CALL M3MESG ( XMSG ) DGATKN_WET = 0 END IF VNAME = 'DGACC_WET' N = INDEX1( VNAME, N_AE_SPC, AE_SPC ) IF ( N .NE. 0 ) THEN V = V + 1 DGACC_WET = V CBLK_MAP( V ) = AE_STRT - 1 + N ELSE XMSG = 'Could not find ' // VNAME // 'in aerosol table' ! CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) CALL M3MESG ( XMSG ) DGACC_DRY = 0 END IF C *** Set additional species contained only in CBLK not needed in CBLK_MAP. V = N_AE_SPC + 1 VSGAT = V V = V + 1 VSGAC = V V = V + 1 VDGAT = V V = V + 1 VDGAC = V V = V + 1 VAT2 = V V = V + 1 VAC2 = V V = V + 1 VAT3 = V V = V + 1 VAC3 = V V = V + 1 VCOR3 = V V = V + 1 VSULF = V V = V + 1 VHNO3 = V V = V + 1 VNH3 = V V = V + 1 VN2O5 = V V = V + 1 VHCL = V C *** Get generic mechanism name IF ( INDEX ( MECHNAME, 'CB4' ) .GT. 0 ) THEN MECH = 'CB4' ELSE IF ( INDEX ( MECHNAME, 'CB05' ) .GT. 0 ) THEN MECH = 'CB05' ELSE IF ( INDEX ( MECHNAME, 'RADM2' ) .GT. 0 ) THEN MECH = 'RADM2' ELSE IF ( INDEX ( MECHNAME, 'SAPRC99' ) .GT. 0 ) THEN MECH = 'SAPRC99' ELSE XMSG = 'Base chemical mechanism must CB4, CB05, RADM2, or SAPRC99' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF C *** Set pointers for gas (vapor) phase species and production rates in CGRID. VNAME = 'SULF' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LSULF = GC_STRT - 1 + GC_G2AE_MAP( N ) MWH2SO4 = GC_MOLWT( GC_G2AE_MAP( N ) ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'HNO3' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LHNO3 = GC_STRT - 1 + GC_G2AE_MAP( N ) MWHNO3 = GC_MOLWT( GC_G2AE_MAP( N ) ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF VNAME = 'N2O5' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LN2O5 = GC_STRT - 1 + GC_G2AE_MAP( N ) MWN2O5 = GC_MOLWT( GC_G2AE_MAP( N ) ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF IF ( MECH .EQ. 'RADM2' .OR. MECH .EQ. 'CB4' & .OR. MECH .EQ. 'CB05') THEN VNAME = 'TERPSP' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LTERP = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE LTERP = 0 XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG ( 'Terpene species not modeled' ) END IF END IF VNAME = 'ALKRXN' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LALKAER = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG & ( 'No production of organic aerosols from alkanes' ) LALKAER = 0 END IF VNAME = 'OLIRXN' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LOLIAER = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas Chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG & ( 'No production of organic aerosols from olefins' ) LOLIAER = 0 END IF VNAME = 'TOLRXN' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LTOLAER = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG & ( 'No production of organic aerosols from toluene' ) LTOLAER = 0 END IF VNAME = 'XYLRXN' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LXYLAER = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG & ( 'No production of organic aerosols from xylene' ) LXYLAER = 0 END IF VNAME = 'CSLRXN' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LCSLAER = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG & ( 'No production of organic aerosols from cresol' ) LCSLAER = 0 END IF VNAME = 'TERPRXN' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LTERPAER = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3MESG ( XMSG ) CALL M3MESG & ( 'No production of organic aerosols from terpenes' ) LTERPAER = 0 END IF VNAME = 'SULPRD' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LSULFP = GC_STRT - 1 + GC_G2AE_MAP( N ) ELSE XMSG = & 'Could not find ' // VNAME // 'in gas chem aerosol table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF C *** For NH3 and HCl, search gas-phase species list first. Then try non- C reactive list. If not found in either list, quit with error message. VNAME = 'NH3' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LNH3 = GC_STRT - 1 + GC_G2AE_MAP( N ) MWNH3 = GC_MOLWT( GC_G2AE_MAP( N ) ) ELSE N = INDEX1( VNAME, N_NR_N2AE, NR_N2AE ) IF ( N .NE. 0 ) THEN LNH3 = NR_STRT - 1 + NR_N2AE_MAP( N ) MWNH3 = NR_MOLWT( NR_N2AE_MAP( N ) ) ELSE XMSG = 'Could not find ' // VNAME // & 'in gas chem aerosol or non-reactives table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF END IF VNAME = 'HCL' N = INDEX1( VNAME, N_GC_G2AE, GC_G2AE ) IF ( N .NE. 0 ) THEN LHCL = GC_STRT - 1 + GC_G2AE_MAP( N ) MWHCL = GC_MOLWT( GC_G2AE_MAP( N ) ) ELSE N = INDEX1( VNAME, N_NR_N2AE, NR_N2AE ) IF ( N .NE. 0 ) THEN LHCL = NR_STRT - 1 + NR_N2AE_MAP( N ) MWHCL = NR_MOLWT( NR_N2AE_MAP( N ) ) ELSE XMSG = 'Could not find ' // VNAME // & 'in gas chem aerosol or non-reactives table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF END IF C *** Set VAPNAMES values VAPNAMES( 1 ) = 'SGTOT_ALK' VAPNAMES( 2 ) = 'SGTOT_OLI_1' VAPNAMES( 3 ) = 'SGTOT_OLI_2' VAPNAMES( 4 ) = 'SGTOT_XYL_1' VAPNAMES( 5 ) = 'SGTOT_XYL_2' VAPNAMES( 6 ) = 'SGTOT_CSL' VAPNAMES( 7 ) = 'SGTOT_TOL_1' VAPNAMES( 8 ) = 'SGTOT_TOL_2' VAPNAMES( 9 ) = 'SGTOT_TRP_1' VAPNAMES( 10 ) = 'SGTOT_TRP_2' C *** Fetch the indices for organic vapors from the non-reactive table DO VV = 1, NCVAP VNAME = VAPNAMES( VV ) N = INDEX1( VNAME, N_NR_N2AE, NR_N2AE ) IF ( N .NE. 0 ) THEN LOCVAP( VV ) = NR_STRT - 1 + NR_N2AE_MAP( N ) ELSE XMSG = 'Could not find ' // VNAME // & 'in non-reactives table' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 ) END IF END DO C *** In the following conversion factors, the 1.0e3 factor C is to convert density from kg m**-3 to g m**-3. C *** factors for converting from ppm to ug m-**3 H2SO4CONV = 1.0E3 * MWH2SO4 / MWAIR HNO3CONV = 1.0E3 * MWHNO3 / MWAIR NH3CONV = 1.0E3 * MWNH3 / MWAIR HCLCONV = 1.0E3 * MWHCL / MWAIR N2O5CONV = 1.0E3 * MWN2O5 / MWAIR C *** reciprocals for converting from ug m**-3 to ppm H2SO4CONV1 = 1.0 / H2SO4CONV HNO3CONV1 = 1.0 / HNO3CONV NH3CONV1 = 1.0 / NH3CONV HCLCONV1 = 1.0 / HCLCONV N2O5CONV1 = 1.0 / N2O5CONV C *** Open the met files. IF ( .NOT. OPEN3( MET_CRO_3D, FSREAD3, PNAME ) ) THEN XMSG = 'Could not open MET_CRO_3D file ' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) END IF IF ( .NOT. OPEN3( MET_CRO_2D, FSREAD3, PNAME ) ) THEN XMSG = 'Could not open MET_CRO_2D file ' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) END IF C *** Set up file structure for visibility file. It has two variables, C visual range in deciview units (dimensionless) and extinction in C units of (1/km) and is for layer 1 only. ! IF ( MYPE .EQ. 0 ) CALL OPVIS ( JDATE, JTIME, TSTEP( 1 ) ) C *** Open the aerosol parameters file (diameters and standard deviations). ! IF ( AERDIAG .AND. ! & MYPE .EQ. 0 ) CALL OPDIAM ( JDATE, JTIME, TSTEP( 1 ) ) C Get domain decomp info from the MET_CRO_3D file CALL SUBHFILE ( MET_CRO_3D, GXOFF, GYOFF, & STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3 ) END IF ! FIRSTIME C ------------ Begin Interpolation of Meteorological Variables --------- MDATE = JDATE MTIME = JTIME MSTEP = TIME2SEC( TSTEP( 2 ) ) CALL NEXTIME ( MDATE, MTIME, SEC2TIME( MSTEP / 2 ) ) WSTEP = WSTEP + TIME2SEC( TSTEP( 2 ) ) IF ( WSTEP .GE. TIME2SEC( TSTEP( 1 ) ) ) WRITETIME = .TRUE. C *** Set floating point synchronization time step: DT = FLOAT( MSTEP ) ! set time step in seconds C *** layered variables PRES, TA, QV, DENS: C *** air density X square root of the determinant of the metric tensor C at midlayer C *** pressure (Pa) VNAME = 'PRES' IF ( .NOT. INTERPX( MET_CRO_3D, VNAME, PNAME, & STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3, 1,NLAYS, & MDATE, MTIME, PRES ) ) THEN XMSG = 'Could not interpolate PRES from ' // MET_CRO_3D CALL M3EXIT ( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) END IF C *** temperature (K) VNAME = 'TA' IF ( .NOT. INTERPX( MET_CRO_3D, VNAME, PNAME, & STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3, 1,NLAYS, & MDATE, MTIME, TA ) ) THEN XMSG = 'Could not interpolate '// VNAME // ' from MET_CRO_3D ' CALL M3EXIT ( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) END IF C *** specific humidity (g H2O / g air) VNAME = 'QV' IF ( .NOT. INTERPX( MET_CRO_3D, VNAME, PNAME, & STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3, 1,NLAYS, & MDATE, MTIME, QV ) ) THEN XMSG = 'Could not interpolate specific humidity from MET_CRO_3D' CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) END IF C *** air density (kg/m3) VNAME = 'DENS' IF ( .NOT. INTERPX( MET_CRO_3D, VNAME, PNAME, & STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3, 1,NLAYS, & MDATE, MTIME, DENS ) ) THEN XMSG = 'Could not interpolate '// VNAME // ' from MET_CRO_3D ' CALL M3EXIT ( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) END IF C *** end interpolation of meteorological variables C --------------------- Begin loops over grid cells -------------------------- C *** initialize CBLK CBLK = 0.0 DO L = 1, NLAYS DO R = 1, MY_NROWS DO C = 1, MY_NCOLS C *** Fetch the grid cell meteorological data. LAYER = L BLKTA = TA ( C,R,L ) BLKPRS = PRES ( C,R,L ) ! Note pascals BLKQA = QV ( C,R,L ) BLKDENS = DENS ( C,R,L ) BLKDENS1 = 1.0 / BLKDENS BLKESAT = ESATL( BLKTA ) BLKEVAP = BLKPRS * BLKQA / ( EPSWATER + BLKQA ) ! BLKRH = MIN( 0.99, BLKEVAP / BLKESAT ) BLKRH = MAX( 0.005, MIN( 0.99, BLKEVAP / BLKESAT ) ) C *** Transfer CGRID values to CBLK (limit to mapped species) DO SPC = 1, N_AE_MAP V = CBLK_MAP( SPC ) CBLK( SPC ) = MAX ( CGRID( C,R,L,V ), CONMIN ) END DO C *** Add gas and vapor phase species to CBLK CBLK( VSULF ) = MAX( CONMIN, & H2SO4CONV * BLKDENS * CGRID( C,R,L,LSULF ) ) CBLK( VHNO3 ) = MAX( CONMIN, & HNO3CONV * BLKDENS * CGRID( C,R,L,LHNO3 ) ) CBLK( VNH3 ) = MAX( CONMIN, & NH3CONV * BLKDENS * CGRID( C,R,L,LNH3 ) ) CBLK( VHCL ) = MAX( CONMIN, & HCLCONV * BLKDENS * CGRID( C,R,L,LHCL ) ) CBLK( VN2O5 ) = MAX( CONMIN, & N2O5CONV * BLKDENS * CGRID( C,R,L,LN2O5 ) ) C *** Fetch gas-phase production rates. C *** sulfate SO4RATE = H2SO4CONV * BLKDENS & * CGRID( C,R,L,LSULFP ) / DT C *** secondary organics DO SPC = 1, NPSPCS ORGPROD( SPC ) = 0.0 END DO IF ( LALKAER .GT. 0 ) & ORGPROD( 1 ) = CGRID( C,R,L,LALKAER ) c..SOA production forom alkenes eliminated c IF ( MECH .EQ. 'RADM2' ) THEN c IF ( LOLIAER .GT. 0 .AND. LTERPAER .GT. 0 ) c & ORGPROD( 2 ) = MAX( 0.0, ( CGRID( C,R,L,LOLIAER ) c & - CGRID( C,R,L,LTERPAER ) ) ) c ELSE IF ( MECH .EQ. 'SAPRC99' ) THEN c IF ( LOLIAER .GT. 0 ) c & ORGPROD( 2 ) = CGRID( C,R,L,LOLIAER ) c END IF IF ( LXYLAER .GT. 0 ) & ORGPROD( 3 ) = CGRID( C,R,L,LXYLAER ) IF ( LCSLAER .GT. 0 ) & ORGPROD( 4 ) = CGRID( C,R,L,LCSLAER ) IF ( LTOLAER .GT. 0 ) & ORGPROD( 5 ) = CGRID( C,R,L,LTOLAER ) IF ( LTERPAER .GT. 0 ) & ORGPROD(6) = CGRID( C,R,L,LTERPAER ) c..Adjust TOL & ALK orgprods for fraction of precursor that produces SOA IF( MECH .EQ. 'SAPRC99' ) THEN ORGPROD( 1 ) = SPALK * ORGPROD( 1 ) ORGPROD( 5 ) = SPTOL * ORGPROD( 5 ) ELSEIF( MECH .EQ. 'RADM2' ) THEN ORGPROD( 1 ) = RDALK * ORGPROD( 1 ) ORGPROD( 5 ) = RDTOL * ORGPROD( 5 ) ENDIF c *** Transfer SVOC concentrations from CGRID to VAPORS DO VV = 1, NCVAP VAPORS( VV ) = 0.0 VVLOC = LOCVAP( VV ) IF ( VVLOC .GT. 0 ) VAPORS( VV ) = CGRID( C,R,L,VVLOC ) END DO ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c call aerosol process routines CALL AEROPROC( NSPCSDA, & CBLK, DT, LAYER, & BLKTA, BLKPRS, BLKDENS, BLKRH, & SO4RATE, & ORGPROD, NPSPCS, VAPORS, NCVAP, & XLM, AMU, & DGATK, DGACC, DGCOR, & XXLSGAT, XXLSGAC, & PMASSAT, PMASSAC, PMASSCO, & PDENSAT, PDENSAC, PDENSCO, & GAMMA_N2O5, & LOGDEV ) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C *** Transfer new aerosol information from CBLK to CGRID c (limit to mapped species) DO SPC = 1, N_AE_MAP V = CBLK_MAP( SPC ) CGRID( C,R,L,V ) = MAX( CONMIN, CBLK( SPC ) ) END DO C *** Transfer new gas/vapor concentrations from CBLK to CGRID CGRID( C,R,L,LSULF ) = MAX( CONMIN, & H2SO4CONV1 * BLKDENS1 * CBLK( VSULF ) ) CGRID( C,R,L,LHNO3 ) = MAX( CONMIN, & HNO3CONV1 * BLKDENS1 * CBLK( VHNO3 ) ) CGRID( C,R,L,LNH3 ) = MAX( CONMIN, & NH3CONV1 * BLKDENS1 * CBLK( VNH3 ) ) CGRID( C,R,L,LHCL ) = MAX( CONMIN, & HCLCONV1 * BLKDENS1 * CBLK( VHCL ) ) CGRID( C,R,L,LN2O5 ) = MAX( CONMIN, & N2O5CONV1 * BLKDENS1 * CBLK( VN2O5 ) ) C *** Zero out species representing the contributions to production C of aerosols. CGRID( C,R,L,LSULFP ) = 0.0 IF ( LALKAER .GT. 0 ) CGRID( C,R,L,LALKAER ) = 0.0 IF ( LOLIAER .GT. 0 ) CGRID( C,R,L,LOLIAER ) = 0.0 IF ( LXYLAER .GT. 0 ) CGRID( C,R,L,LXYLAER ) = 0.0 IF ( LCSLAER .GT. 0 ) CGRID( C,R,L,LCSLAER ) = 0.0 IF ( LTOLAER .GT. 0 ) CGRID( C,R,L,LTOLAER ) = 0.0 IF ( LTERPAER .GT. 0 ) CGRID( C,R,L,LTERPAER ) = 0.0 C *** Transfer new SVOC concentrations from VAPORS to CGRID DO VV = 1, NCVAP VVLOC = LOCVAP( VV ) IF ( VVLOC .GT. 0 ) & CGRID( C,R,L,VVLOC ) = MAX ( CONMIN, VAPORS( VV ) ) END DO C *** Calculate volume fraction of each mode < 2.5um aerodynamic diameter XXLSGCO = LOG( SGINICO ) CALL INLET25 ( DGATK, XXLSGAT, PDENSAT, PM25AT ) CALL INLET25 ( DGACC, XXLSGAC, PDENSAT, PM25AC ) CALL INLET25 ( DGCOR, XXLSGCO, PDENSCO, PM25CO ) C *** Write aerosol extinction coefficients and deciviews to visibility C diagnostic array (lowest vertical layer only) ! IF ( WRITETIME .AND. L .EQ. 1 ) THEN IF ( WRITETIME ) THEN CALL GETVISBY ( NSPCSDA, & CBLK, BLKRH, & BLKDCV1, BLKEXT1, BLKDCV2, BLKEXT2, & DGATK, DGACC, DGCOR, & XXLSGAT, XXLSGAC, & PMASSAT, PMASSAC, PMASSCO ) ! VIS_SPC( C,R,L,IDCVW1 ) = BLKDCV1 ! visual range [ deciview ] ! (Mie) VIS_SPC( C,R,L,IBEXT1 ) = BLKEXT1 ! aerosol extinction [ 1/km ] ! (Mie) ! VIS_SPC( C,R,L,IDCVW2 ) = BLKDCV2 ! visual range [ deciview ] ! (Reconstructed) VIS_SPC( C,R,L,IBEXT2 ) = BLKEXT2 ! aerosol extinction [ 1/km ] ! (Reconstructed) END IF C *** Write wet diameters, 2nd, and 3rd moments to aerosol diagnostic array C This assumes that GETPAR was last called with M3_WET_FLAG = .TRUE. IF ( WRITETIME .AND. AER_DIAG ) THEN DIAM_SPC( C,R,L, 5 ) = DGATK*1.0e6 ! wet i-mode diameter DIAM_SPC( C,R,L, 6 ) = DGACC*1.0e6 ! wet j-mode diameter DIAM_SPC( C,R,L, 8 ) = CBLK( VAT2 ) ! wet i-mode 2nd moment DIAM_SPC( C,R,L, 9 ) = CBLK( VAC2 ) ! wet j-mode 2nd moment DIAM_SPC( C,R,L,12 ) = CBLK( VAT3 ) ! wet i-mode 3rd moment DIAM_SPC( C,R,L,13 ) = CBLK( VAC3 ) ! wet j-mode 3rd moment DIAM_SPC( C,R,L,16 ) = PM25AT ! i-mode fine fraction DIAM_SPC( C,R,L,17 ) = PM25AC ! j-mode fine fraction DIAM_SPC( C,R,L,18 ) = PM25CO ! coarse-mode fine fraction DIAM_SPC( C,R,L,19 ) = GAMMA_N2O5 ! N2O5 heterorxn probability END IF ! WRITETIME .AND. AER_DIAG C *** Add wet mean diameters to CGRID CGRID( C,R,L,CBLK_MAP( DGATKN_WET ) ) = DGATK CGRID( C,R,L,CBLK_MAP( DGACC_WET ) ) = DGACC C *** Calculate 2nd and 3rd moments of the "dry" aerosol distribution C NOTE! "dry" aerosol excludes both H2O and SOA (January 2004 --SJR) C Aitken mode. M3_WET = CBLK( VAT3 ) M3_DRY = M3_WET & - H2OFAC * CBLK( VH2OAI ) & - ORGFAC * CBLK( VORGAI ) & - ORGFAC * CBLK( VORGBAI ) M2_WET = CBLK( VAT2 ) M2_DRY = M2_WET * ( M3_DRY / M3_WET ) ** TWO3 CBLK( VAT3 ) = M3_DRY CBLK( VAT2 ) = M2_DRY C accumulation mode. M3_WET = CBLK( VAC3 ) M3_DRY = M3_WET & - H2OFAC * CBLK( VH2OAJ ) & - ORGFAC * CBLK( VORGAJ ) & - ORGFAC * CBLK( VORGBAJ ) M2_WET = CBLK( VAC2 ) M2_DRY = M2_WET * ( M3_DRY / M3_WET ) ** TWO3 CBLK( VAC3 ) = M3_DRY CBLK( VAC2 ) = M2_DRY C *** Calculate geometric mean diameters and standard deviations of the C "dry" size distribution CALL GETPAR( NSPCSDA, & CBLK, & PMASSAT, PMASSAC, PMASSCO, & PDENSAT, PDENSAC, PDENSCO, & DGATK, DGACC, DGCOR, & XXLSGAT, XXLSGAC, & M3_WET_FLAG ) C *** Write dry aerosol distribution parameters to aerosol diagnostic array IF ( WRITETIME .AND. AER_DIAG ) THEN DIAM_SPC( C,R,L, 1 ) = EXP( XXLSGAT ) DIAM_SPC( C,R,L, 2 ) = EXP( XXLSGAC ) DIAM_SPC( C,R,L, 3 ) = DGATK*1.0e6 ! dry i-mode diameter DIAM_SPC( C,R,L, 4 ) = DGACC*1.0e6 ! dry j-mode diameter DIAM_SPC( C,R,L, 7 ) = DGCOR*1.0e6 ! dry coarse-mode diam. DIAM_SPC( C,R,L,10 ) = CBLK( VAT3 ) ! dry i-mode 3rd moment DIAM_SPC( C,R,L,11 ) = CBLK( VAC3 ) ! dry j-mode 3rd moment DIAM_SPC( C,R,L,14 ) = CBLK( VCOR3 ) ! dry coarse-mode 3rd ! moment DIAM_SPC( C,R,L,15 ) = BLKRH ! relative humidity END IF ! WRITETIME .AND. AER_DIAG C *** Calculate aerosol surface area from the dry 2nd moment. Dry value is C used in transport routines. CGRID( C,R,L,CBLK_MAP( VSURFAT ) ) = PI * CBLK( VAT2 ) CGRID( C,R,L,CBLK_MAP( VSURFAC ) ) = PI * CBLK( VAC2 ) C *** Add dry mean diameters to CGRID CGRID( C,R,L,CBLK_MAP( DGATKN_DRY ) ) = DGATK CGRID( C,R,L,CBLK_MAP( DGACC_DRY ) ) = DGACC END DO ! loop on MY_COLS END DO ! loop on MY_ROWS END DO ! loop on NLAYS C *** end of loops over grid cells C ------------ Write diagnostic information to output files -------------- C *** If last call this hour, write visibility information. IF ( WRITETIME ) THEN MDATE = JDATE MTIME = JTIME CALL NEXTIME ( MDATE, MTIME, TSTEP( 2 ) ) WSTEP = 0 WRITETIME = .FALSE. ! IF ( .NOT. WRITE3( CTM_VIS_1, ALLVAR3, ! & MDATE, MTIME, VIS_SPC ) ) THEN ! XMSG = 'Could not write ' // CTM_VIS_1 // ' file' ! CALL M3EXIT ( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) ! END IF ! WRITE( LOGDEV, '( /5X, 3( A, :, 1X ), I8, ":", I6.6 )' ) ! & 'Timestep written to', CTM_VIS_1, ! & 'for date and time', MDATE, MTIME C *** Write data to the aerosol parameters file. ! IF ( AERDIAG ) THEN ! IF ( .NOT. WRITE3( CTM_DIAM_1, ALLVAR3, ! & MDATE, MTIME, DIAM_SPC ) ) THEN ! XMSG = 'Could not write ' // CTM_DIAM_1 // ' file' ! CALL M3EXIT ( PNAME, MDATE, MTIME, XMSG, XSTAT1 ) ! END IF ! WRITE( LOGDEV, '( /5X, 3( A, :, 1X ), I8, ":", I6.6 )' ) ! & 'Timestep written to', CTM_DIAM_1, ! & 'for date and time', MDATE, MTIME ! END IF END IF RETURN END