subroutine chemdriv(m1,m2,m3,i0,j0,ia,iz,ja,jz,igrd, & ncol,nrow,nlay,dt,itzon,idfin,dx,dy, & cldtrns,water,tempk,tsurf,press,height,uwind, & vwind,cwc,pwr,cph,aph,ajno2, & cigwtr,conc,cellat,cellon, & snow,snowrat,snowfrc,snowalb,topo,fsurf,albedo, & lrdlai,lai,ldark,solmas,vegmas,reemis,l3dflag, & iptr2d,iptrsa,ipa_cel,iproc_id) use filunit use chmstry use o3colmap use bndary use camxcom use procan use tracer use rtracchm use rtcmcchm c c----CAMx v7Beta6 190902 c c CHEMDRIV performs chemistry on the current grid for one time step. c It calls one of two ODE solvers: EBISOLV or LSCALL. c EBISOLV is the driver for the EBI solver with Hertel's solutions c LSCALL is the driver for the LSODE solver c c Mechanism specific subroutines are passed to the driver routines. c c Local array element con(nspec+1) is used to store the concentrations c of non-used species, i.e., species that are in the chem solver but c not on the species list for this run c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 1/9/02 A minor bug fix related to fast aero routine c (units conversion) c 1/15/02 Added code to handle RTRAC c 01/30/02 Added code for RTRAC probing tool c 10/18/02 Updated AEROCHEM routine; added CWC for aqueous c PM chemistry c 03/14/03 Implemented Hg chemistry (Prakash Karamchandani, AER) c 03/21/03 Removed the OSAT technology type OPPAT c 04/2/03 Removed option for UAM-V type cloud adjustment c 07/21/03 Added optional snow effect on UV albedo c 04/21/04 Added CMU sectional aerosol chemistry c 10/14/04 Added Mechanism 10 c 10/05/05 Added PA calls for Mech 4 c 10/06/05 Removed Mech 2 c 10/07/05 Cloud water content below 273 K is partitioned to c liquid fraction like in wet dep c 10/27/05 Added cloud water pH as a 3-D argument c 01/12/06 Aqueous PM now called each timestep c 11/28/06 Fixed bug in assigning ldchm flag c 07/04/07 Added heterogeneous hydrolysis of N2O5 c 07/11/07 Added code for RTCMC probing tool c 12/26/06 -bkoo- Removed AEROCHEM_AQ (merged into AEROCHEM_CF/CMU) c 12/29/06 -bkoo- Added SOACHEM c 01/08/07 -bkoo- Added Mechanism 6 (CB05) c Fixed titration rxn # for SAPRC99 (PA) c 07/04/07 -bkoo- Added heterogeneous hydrolysis of N2O5 c 07/16/07 -bkoo- Revised for HDDM c 04/24/08 -gyarwood- Added EBI solver c 06/01/08 -bkoo- Modified call to EBI for DDM c 06/11/08 -bkoo- Revised for rate constant sensitivity c 07/16/08 -bkoo- Added DDM turn-off flag c 01/30/09 -bkoo- Removed the CMC fast solver c 12/08/09 -jjung- Added new RTRAC surface model and c partitioning algorithm c 05/27/09 -jjung- Added Mechanism 7 (CB6) c 7/14/10 -cemery- Added code for in-line TUV cloud adjustment c 06/02/11 -pkaramchandani- Updated Hg chemistry code c 04/02/12 -cemery- Removed RADM cloud adjustment option, cloud/aerosol c adjustments now always done with in-line TUV; AHO c file is now just ozone column; replaced haze c dimension with terrain height in photo file c 04/04/12 -cemery- Introduced 4-dim interpolation of J values over c zenith, height, albedo, terrain c 4/12/12 -cemery- Added T and P adjustments to photolysis rates c 07/23/12 -unopmongcol- Added Mechanism 1 (CB6-I) c 10/08/12 -jjung- Added Mechanism 2 (CB6r1) c 04/29/13 -cemery- Added surface model c 05/24/13 -cemery- RTCMC photolysis can point to core reactions c 06/28/13 -bkoo- Changed the last argument to LSCALL c 08/23/13 -bkoo- Revised for DDM-PM c 09/21/13 -bkoo- Updated DDM-PM for KHETERO c 10/14/13 -bkoo- Added pressure to KHETERO call c 01/03/14 -bkoo- Increased SDDM array size for organic nitrate c hydrolysis c 03/18/14 -bkoo- Revised for benzene SOA c 06/27/14 -bkoo- Added Mechanism 3 (CB6r2h) c 08/11/14 -cemery- Albedo now includes effects of snow c 08/22/14 -cemery- Extended Zhang LU to surface model, added snow c 10/09/14 -cemery- Added subgrid convective model c 11/11/14 -bkoo- Added Mechanism 4 (CB6r3) c 12/07/14 -bkoo- Revised for VBS c 09/11/15 -bkoo- Revised for SA v3 c 01/08/16 -bkoo- Updated for DDMRATE c Revised VBS c Removal of IEH solver c 02/29/16 -bkoo- Revised SDDM array for hetero N2O5 + HCL rxn c 07/05/16 -gyarwood- Added HR_XO c 07/20/16 -bkoo- Revised SDDM array for hetero INTR hydrolysis c 08/25/16 -bkoo- Updated for new SOAP c 09/02/16 -bkoo- Set SDDM array size with # of het rxns c 09/16/16 -pk- Added LDARK to AEROCHEM_CF & AEROCHEM_CMU calls c 11/28/16 -bkoo- Revised VBS precursors c 10/13/17 -bkoo- Added APH & AERPH to track aerosol pH c 10/16/17 -bkoo- Added AJNO2 to track time-weighted NO2 photolysis rate c Expanded DELCON to store SOA-specific changes c 01/12/18 -bkoo- Removed BNZA/TOLA/XYLA/ISP/TRP c 08/09/18 -bkoo- Revised for rate term sensitivity c c Input arguments: c igrd grid index c ncol number of columns c nrow number of rows c nlay number of layers c dt timestep (s) c itzon time zone c idfin map of nested grids in this grid c dx cell size in x-direction (m) c dy cell size in y-direction (m) c cldtrns energy transmission coefficient (fraction) c water water vapor (ppm) c tempk temperature (K) c tsurf surface temperature (K) c press pressure (mb) c height layer interface height (m) c uwind u-component wind (m/s) c vwind v-component wind (m/s) c cwc cloud water content (g/m3) c pwr precipitation water content (g/m3) c cph cloud water pH c aph aerosol pH c ajno2 time-weighted accumulated NO2 photolysis rate ([1/hr]*[hr]) c cigwtr subgrid cloud water (g/m3) c conc species concentration (umol/m3,ug/m3) c cellat cell centroid latitude (deg) c cellon cell centroid longitude (deg) c snow snow cover water equivalent (m) c snowrat snow water accumulation rate (m/hr) c snowfrc snow cover fraction c snowalb snow albedo c topo cell terrain altitude (m MSL) c fsurf landuse fraction c albedo surface albedo c lrdlai LAI was read c lai gridded LAI c ldark darkness flag (T=dark) c solmas surface model soil mass (mol/ha) c vegmas surface model vegetation mass (mol/ha) c l3dflag save 3-D average concentrations c iptr2d pointers into vectors for 2-D fields c iptrsa pointers into vectors for tracer conc c ipa_cel gridded array to identify if cell is c in a IPRM sub-domain c iproc_id process ID for this slice c c Output arguments: c conc species concentration (umol/m3,ug/m3) c solmas surface model soil mass (mol/ha) c vegmas surface model vegetation mass (mol/ha) c reemiss surface model reemission rate (mol/s) c c Routines called: c CHEM10 c KTHERM c KHETERO c GETZNTH c KPHOTO c RECALIB c LOADDM c HDDMADJ c EBISOLV c LSCALL c HDDMCHEM c VBSCHEM c SOACHEM c LOADDMPM c AEROCHEM_CF c AEROCHEM_CMU c CALDATE c SRFMOD c HGGASCHEM c HGAQSCHEM c HGADSORB c O3PRDSA c OSATSA c PSATSA c REPARTSA c CHEMRT c PARTITIONRT c SRFMODRT c KRTC c KSCI c DLSDRV c SLSDRV c RBKDRV c CPAMECH c CPARAD c UPTIME c CAMXERR c c Called by: c CHEMRXN c implicit none include "camx.prm" include 'camx_aero.inc' include 'soap.inc' include 'vbs.inc' include 'deposit.inc' include 'flags.inc' include 'chmdbg.inc' include 'rtracsrf.inc' c c----Argument declarations c integer :: m1,m2,m3,i0,j0,ia,iz,ja,jz integer igrd integer ncol integer nrow integer nlay real dt integer itzon integer idfin(m1,m2) integer iproc_id c real, dimension(m1,m2,m3) :: cldtrns,cph,aph,ajno2, & water,tempk,press,cwc,pwr,uwind,vwind,cigwtr real, dimension(m1,m2,m3) :: height c real conc(m1,m2,m3,nspec) c real cellat(m1,m2) real cellon(m1,m2) real snow(m1,m2) real snowrat(m1,m2) real snowfrc(m1,m2) real snowalb(m1,m2) real topo(m1,m2) real tsurf(m1,m2) real fsurf(m1,m2,NLU) real albedo(m1,m2) real lai(m1,m2) real dx(nrow) real dy c real solmas(m1,m2,nsmspc) real vegmas(m1,m2,nsmspc) real reemis(m1,m2,nsmspc) c logical ldark(m1,m2) logical lrdlai logical l3dflag integer iptr2d integer*8 iptrsa integer ipa_cel(ncol,nrow,nlay) c c======================== Process Analysis Begin ======================= c logical ldoirr, ldoipr, ltrapirr real titrt real chg0, chg2, thg20 c real rrxn_irr(MXREACT) real patmp(MXTRSP) c c======================== Process Analysis End ========================= c c========================= Hg Chemistry Begin ========================== c real cl2c ! Cl2 concentration (ppm) real chgiip, chgiipc ! Adsorbed Hg concentrations (ug/m3) real brc, broc ! Br and BrO concentrations (ppm) real henso2,tfso2 ! Henry's Law for SO2 (M/atm) real heno3,tfo3 ! Henry's Law for O3 (M/atm) real pm10 ! Total PM10 concentrations (ug/m3) integer iocean ! Ocean flag: 1 = ocean, 0 = other c real thg2 ! Total Hg(2) real convhgfac ! conv factor for hg (ppm to ug/m3) c c-----variables for Hg adsorption on primary particles real, parameter :: MWNA = 23.0 ! Mol. wt. of sodium real, parameter :: MWCL = 35.5 ! Mol. wt. of chlorine atom real, parameter :: MWNACL = MWNA + MWCL ! Mol. wt. of seasalt real, parameter :: RHOSEAS = 2.0E6 ! Density of seasalt (g/m3) real dfine, dcoarse ! diameter of fine and coarse particles (um) c c-----surface area concentration (m2/m3) and volume concentrations (cm3/m3) of c-----fine and coarse particles real safine,sacoarse,vfine,vcoarse c real cna,cpcl,vseas,fseas ! Na,Cl concs; vol. conc of seasalt, fraction ! of fine PM that is seasalt c c========================== Hg Chemistry End =========================== c c-----variables for SOACHEM/VBSCHEM c real avgox(NOXID) ! Average oxidant conc. (ppm) real con_prec(NPREC) ! Precursor conc. (ppm) real con_CG(NCG) ! CG conc. (ppm) real con_NV(NNV) ! Non-volatile product conc. (ppm) real yfac(3,NPREC_S) ! 3 CG products per each precursor c c======================== DDM Begin ======================= c real sddm(MXTRSP+NHETRXN*MXFDDM) real avgsen(MXFDDM*NOXID) ! Average sensitivities of oxidants (ppm) real sen_prec(MXFDDM*NPREC) ! Sensitivities of precursors (ppm) real sen_CG(MXFDDM*NCG) ! Sensitivities of CGs (ppm) real sen_NV(MXFDDM*NNV) ! Sensitivities of NVs (ppm) c c======================== DDM End ======================= c c-----subroutine names to be called by EBISOLV and LSCALL c external ebirxn2, ebirate2, hr_hox2, hr_nox2, hr_nxy2, hr_pan2 external ebirxn3, ebirate3, hr_hox3, hr_nox3, hr_nxy3, hr_pan3 external ebirxn4, ebirate4, hr_hox4, hr_nox4, hr_nxy4, hr_pan4 external ebirxn5, ebirate5, hr_hox5, hr_nox5, hr_nxy5, hr_pan5 external ebirxn6, ebirate6, hr_hox6, hr_nox6, hr_nxy6, hr_pan6 external hr_xo3, hr_xo4, hr_xo_dum external lsrate2 external lsrate3 external lsrate4 external lsrate5 external lsrate6 external ddmrate2, hddmjac2 external ddmrate3, hddmjac3 external ddmrate4, hddmjac4 external ddmrate5, hddmjac5 external ddmrate6, hddmjac6 c c-----Local variables c integer i,j,k,is,l,ij,irxn,ihg,iht,nj,nz,nn integer*8 idx integer ierr integer aero_flag integer iozon,irt_cel integer ipa_idx integer lsm integer nitsum(MXLAYER) integer icig c real dtchem,dtfact real tcell,pcell,convfac,cldph,aerph,ctrns,zenith,cliq real o2,ch4,atm,h2 real con(MXSPEC+1) real avgcnc(MXSPEC+1) real hght,alb,trn real desto3,prdo3n,prdo3v real depth,press0,ucomp,vcomp,wind real avit(MXCELLS,MXCELLS) real sumit(MXLAYER) c c --- the delcon array stores the delta concentrations at c various stages of the chemistry process. c c delcon(1,*) -- change in concentration after just gas phase c chemistry (no aerosol chemistry applied) c delcon(2,*) -- change in concentration after gas phase c and aqueous phase chemistry c delcon(3,*) -- change in concentration after all chemistry c is complete c delcon(4,*) -- special instance to capture the change in c concentration due to gas phase and aqueous c phase but also includes the transfer of c some nitrate species during the aqueous phase c delcon(5,*) -- change in concentration due to SOA polymerization c delcon(6,*) -- change in concentration due to particle-phase SOA photolysis c delcon(7,*) -- change in concentration due to in-cloud SOA formation c integer ispc real o3old, o3new c real delcon(7,MXSPEC+1) real cold(MXSPEC) real cnew(MXSPEC) real rtcon(MXTRSP) real soil(MXSMSPC) real vegy(MXSMSPC) real reemit(MXSMSPC) real flu(NLU) c c-----Entry point c call flush(6) call flush(iout) c dtchem = dt/3600. con(nspec+1) = 0. ltrapirr = .FALSE. c c-----Aerosol timing parameters: c time_aero : time to call aerosol routine (HHMM) c date_aero : date to call aerosol routine (YYJJJ) c dt_aero : time interval between calls to aerosol routines (min) c aero_dt : actual time interval for each grid (hr) (accumulated dtchem) c aero_flag = 0 if( naero.GT.0 .AND. (aeropt.EQ.'CF' .OR. aeropt.EQ.'CMU')) then aero_flag = 1 aero_dt(igrd) = aero_dt(igrd) + dtchem if ( (timec(igrd)-time_aero(igrd)) .ge. -0.01 .and. & (datec(igrd) .eq. date_aero(igrd)) ) then aero_flag = 2 endif endif c igrdchm = igrd do 91 k = 1,nlay c c---Update list of variables private to each thread for updated c---Hg treatment, PK, Ramboll 06/11 c nitsum(k)=1 sumit(k)=1. c c$omp parallel default(shared) c$omp& private(i,j,l,is,ispc,con, c$omp& ij,iozon,trn,alb,hght,ctrns, c$omp& zenith,delcon,tcell,pcell,cliq,cldph,aerph,atm,O2,H2,CH4, c$omp& sddm,ldoipr,ldoirr,ipa_idx,titrt,irxn,rrxn_irr, c$omp& ierr,convfac,convhgfac,patmp,idx,nn,cold,cnew,avgcnc,o3old, c$omp& o3new,irt_cel,iocean,brc,broc,cl2c,henso2,tfso2,heno3,tfo3, c$omp& pm10,iht,ihg,ltrapirr,chg0,chg2,chgiip,chgiipc,thg2, c$omp& thg20,cna,cpcl,vseas,fseas,dfine,dcoarse,vfine,vcoarse, c$omp& safine,sacoarse,avgox,con_prec,con_CG,con_NV,yfac, c$omp& avgsen,sen_prec,sen_CG,sen_NV, c$omp& nj,nz,rtcon,desto3,prdo3n,prdo3v, c$omp& depth,press0,vcomp,ucomp,wind,lsm, c$omp& flu,soil,vegy,reemit) c$omp& copyin(/ijkgrd/) c c$omp do schedule(guided) collapse(2) c do 90 j = ja,jz !2,nrow-1 do 89 i = ia,iz !2,ncol-1 ltrapirr = .FALSE. ichm = i jchm = j kchm = k avit(i,j) = 0. c c-----skip chemistry if fine grid exists in this cell c if (idfin(i,j).gt.igrd) goto 89 c c-----For the gas phase species (numbered 1 to ngas) c Pass concentration to CON, converting from umol/m3 to ppm c tcell = tempk(i,j,k) pcell = press(i,j,k) convfac = densfac*(273./tcell)*(pcell/1013.) cldph = cph(i,j,k) aerph = aph(i,j,k) tchm = tcell wchm = water(i,j,k) do is = 1,nrad con(is) = conc(i,j,k,is) if (con(is).lt.0.) then write(iout,'(//,a)') 'ERROR in CHEMDRIV:' write(iout,*) 'Negative concentration before chem' write(iout,*) 'igrd, i, j, k = ', igrd,i,j,k do l = 1,nspec write(iout,'(i3,2x,a7,e10.3)') l,spname(l),con(l) enddo call camxerr() endif con(is) = amax1(bdnl(is),con(is)) enddo do is = nrad+1,ngas con(is) = conc(i,j,k,is)/convfac if (con(is).lt.0.) then write(iout,'(//,a)') 'ERROR in CHEMDRIV:' write(iout,*) 'Negative concentration before chem' write(iout,*) 'igrd, i, j, k = ', igrd,i,j,k do l = 1,nspec write(iout,'(i3,2x,a7,e10.3)') l,spname(l),con(l) enddo call camxerr() endif con(is) = amax1(bdnl(is),con(is)) enddo c c-----Load any aerosols c if (ngas.lt.nspec) then do is=ngas+1,nspec con(is) = conc(i,j,k,is) if (con(is).lt.0.) then write(iout,'(//,a)') 'ERROR in CHEMDRIV:' write(iout,*) 'Negative concentration before chem' write(iout,*) 'igrd, i, j, k = ', igrd,i,j,k do l = 1,nspec write(iout,'(i3,2x,a7,e10.3)') l,spname(l),con(l) enddo call camxerr() endif con(is) = amax1(bdnl(is),con(is)) enddo endif c c-----Mechanism 10 is done here c if (idmech.eq.10) then call chem10(dtchem,tcell,pcell,water(i,j,k),con) goto 88 endif c c======================== DDM Begin ==================================== c c ---- load the sensitivities for this cell into 2-D array --- c if ( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then call loaddm(.FALSE.,m1,m2,m3,ntotsp,ptconc(iptrsa), & i,j,k,nddmsp,nspec,ngas,nrad,NHETRXN,sddm, & convfac) if ( lhddm ) & call hddmadj(.FALSE.,nddmsp,con,sddm) endif c c======================== DDM End ======================================= c c-----Determine thermal rate constants c call ktherm(tcell,pcell) c c-----Calculate the rate constant for heterogeneous hydrolysis of N2O5 c call khetero(tcell,pcell,water(i,j,k),con, & (lddm.AND.lddmcalc(igrd)),nddmsp,sddm) c c-----Load local values of ozone, terrain ht, albedo, zenith angle, altitude c call getznth(cellat(i,j),cellon(i,j),timec(igrd), & datec(igrd),itzon,zenith,ldark(i,j)) ldchm = ldark(i,j) c alb = albedo(i,j) ij = i + (j-1)*m1 iozon = icdozn(iptr2d-1+ij) iocean = icdocn(iptr2d-1+ij) c trn = topo(i,j)/1000. hght = height(i,j,k)/2000. if (k.gt.1) hght = (height(i,j,k) + height(i,j,k-1))/2000. ctrns = 0. if (.not.ldark(i,j)) ctrns = cldtrns(i,j,k) c c-----Determine photolysis rates through interpolation of look-up table c call kphoto(iozon,alb,trn,hght,zenith,ctrns,ldark(i,j), & tcell,pcell) c c-----Accumulate JNO2 ([1/hr]*[hr]) c ajno2(i,j,k) = ajno2(i,j,k) + rk(jno2rxn) * dtchem c c======================== Source Apportion Begin ======================= c c --- store the current value of the species needed for tracer "chemistry" ---- c if( ltrace ) then o3old = con(ko3) do ispc=1,nrad delcon(1,ispc) = -con(ispc)*convfac delcon(2,ispc) = -con(ispc)*convfac delcon(3,ispc) = -con(ispc)*convfac delcon(4,ispc) = -con(ispc)*convfac cold(ispc) = con(ispc) enddo do ispc=nrad+1,ngas delcon(1,ispc) = -con(ispc)*convfac delcon(2,ispc) = -con(ispc)*convfac delcon(3,ispc) = -con(ispc)*convfac delcon(4,ispc) = -con(ispc)*convfac cold(ispc) = con(ispc)*convfac enddo do ispc=ngas+1,nspec delcon(1,ispc) = -con(ispc) delcon(2,ispc) = -con(ispc) delcon(3,ispc) = -con(ispc) delcon(4,ispc) = -con(ispc) cold(ispc) = con(ispc) enddo delcon(5:7,:) = 0.0 ! reset to zero delcon(1:7,nspec+1) = 0.0 ! species not in mechanism c c --- call routine to recalibrate so that tracer species stays c on track with model species: needed because the numerics c cause tracer species to drift c if( tectyp .NE. RTRAC .AND. tectyp .NE. RTCMC ) & call recalib(m1,m2,m3,ntotsp, & ptconc(iptrsa),i,j,k,cold) c c --- if ozone or nitrate tracers are requested, c set the flag for process analysis -- this will save needed c information from the solver routines --- c c --- set variables for OSAT3 ---- c if( lozone .OR. lnitrate ) then ltrapirr = .TRUE. do irxn=1,nreact rrxn_irr(irxn) = 0.0 enddo endif endif c c========================= Source Apportion End ======================== c c-----Assign SOA precursor concentrations before chemistry integration c if( aero_flag .GT. 0 ) then if ( lvbs ) then if (idmech.eq.2 .or. idmech.eq.3 & .or. idmech.eq.4) then ! cb6 con_prec(1) = con(kBENZ) con_prec(2) = con(kTOL ) con_prec(3) = con(kXYL ) con_prec(4) = con(kISOP) con_prec(5) = con(kTERP) con_prec(6) = con(kSQT ) elseif (idmech.eq.5) then ! saprc07tc con_prec(1) = con(kBENZ) con_prec(2) = con(kTOLU) + con(kARO1) con_prec(3) = con(kOXYL) + con(kMXYL) + con(kPXYL) & + con(kB124) + con(kARO2) con_prec(4) = con(kISOP) con_prec(5) = con(kAPIN) + con(kTERP) con_prec(6) = con(kSESQ) elseif (idmech.eq.6) then ! cb05 con_prec(1) = 0.0 con_prec(2) = con(kTOL ) con_prec(3) = con(kXYL ) con_prec(4) = con(kISOP) con_prec(5) = con(kTERP) con_prec(6) = con(kSQT ) endif con_prec( 7) = con(kIVOG) con_prec( 8) = con(kIVOD) con_prec( 9) = con(kIVOA) con_prec(10) = con(kIVOB) else ! soap if (idmech.eq.2 .or. idmech.eq.3 & .or. idmech.eq.4) then ! cb6 con_prec(1) = con(kBENZ) con_prec(2) = con(kTOL ) con_prec(3) = con(kXYL ) con_prec(5) = con(kISOP) con_prec(6) = con(kTERP) con_prec(7) = con(kSQT ) elseif (idmech.eq.5) then ! saprc07tc con_prec(1) = con(kBENZ) con_prec(2) = con(kTOLU) + con(kARO1) con_prec(3) = con(kOXYL) + con(kMXYL) + con(kPXYL) & + con(kB124) + con(kARO2) con_prec(5) = con(kISOP) con_prec(6) = con(kAPIN) + con(kTERP) con_prec(7) = con(kSESQ) elseif (idmech.eq.6) then ! cb05 con_prec(1) = 0.0 con_prec(2) = con(kTOL ) con_prec(3) = con(kXYL ) con_prec(5) = con(kISOP) con_prec(6) = con(kTERP) con_prec(7) = con(kSQT ) endif con_prec(4) = con(kIVOA) c c======================== DDM Begin ======================= c if ( lddm ) then if (idmech.eq.2 .or. idmech.eq.3 & .or. idmech.eq.4) then ! cb6 do is = 1, nddmsp sen_prec( is) = sddm(iptddm(kBENZ)-1+is) sen_prec( nddmsp+is) = sddm(iptddm(kTOL )-1+is) sen_prec(2*nddmsp+is) = sddm(iptddm(kXYL )-1+is) sen_prec(4*nddmsp+is) = sddm(iptddm(kISOP)-1+is) sen_prec(5*nddmsp+is) = sddm(iptddm(kTERP)-1+is) sen_prec(6*nddmsp+is) = sddm(iptddm(kSQT )-1+is) enddo elseif (idmech.eq.5) then ! saprc07tc do is = 1, nddmsp sen_prec( is) = sddm(iptddm(kBENZ)-1+is) sen_prec( nddmsp+is) = sddm(iptddm(kTOLU)-1+is) & + sddm(iptddm(kARO1)-1+is) sen_prec(2*nddmsp+is) = sddm(iptddm(kOXYL)-1+is) & + sddm(iptddm(kMXYL)-1+is) & + sddm(iptddm(kPXYL)-1+is) & + sddm(iptddm(kB124)-1+is) & + sddm(iptddm(kARO2)-1+is) sen_prec(4*nddmsp+is) = sddm(iptddm(kISOP)-1+is) sen_prec(5*nddmsp+is) = sddm(iptddm(kAPIN)-1+is) & + sddm(iptddm(kTERP)-1+is) sen_prec(6*nddmsp+is) = sddm(iptddm(kSESQ)-1+is) enddo elseif (idmech.eq.6) then ! cb05 do is = 1, nddmsp sen_prec( is) = 0.0 sen_prec( nddmsp+is) = sddm(iptddm(kTOL )-1+is) sen_prec(2*nddmsp+is) = sddm(iptddm(kXYL )-1+is) sen_prec(4*nddmsp+is) = sddm(iptddm(kISOP)-1+is) sen_prec(5*nddmsp+is) = sddm(iptddm(kTERP)-1+is) sen_prec(6*nddmsp+is) = sddm(iptddm(kSQT )-1+is) enddo endif do is = 1, nddmsp sen_prec(3*nddmsp+is) = sddm(iptddm(kIVOA)-1+is) enddo endif c c======================== DDM End ======================= c endif endif c c-----Chemistry integration, pass subroutines for mechanism used c atm = 1.e6 O2 = 2.095e5 CH4 = 1.75 H2 = 0.60 c c======================== Process Analysis Begin ======================= c ldoipr = .FALSE. ldoirr = .FALSE. if( lproca ) then if( lipr ) then if( ipa_cel(i+i0,j+j0,k) .GT. 0 ) then ipa_idx = ipa_cel(i+i0,j+j0,k) ldoipr = .TRUE. endif endif if( lirr ) then ltrapirr = .TRUE. if( ipa_cel(i+i0,j+j0,k) .GT. 0 ) then ipa_idx = ipa_cel(i+i0,j+j0,k) ldoirr = .TRUE. endif titrt = 0.0 do irxn=1,nreact rrxn_irr(irxn) = 0.0 enddo endif endif c c======================== Process Analysis End ========================= c if ( idsolv .EQ. IDEBI ) then if (idmech.eq.2) then call ebisolv(ebirxn2,ebirate2,hr_hox2,hr_nox2, & hr_nxy2,hr_pan2,hr_xo_dum, & dtchem,water(i,j,k),atm,O2, & CH4,H2,tcell,ldark(i,j),con,avgcnc,avit(i,j), & hddmjac2,nddmsp,sddm,NOXID,avgsen, & (lddm.AND.lddmcalc(igrd)), & nirrrxn,rrxn_irr,ltrapirr) elseif (idmech.eq.3) then call ebisolv(ebirxn3,ebirate3,hr_hox3,hr_nox3, & hr_nxy3,hr_pan3,hr_xo3, & dtchem,water(i,j,k),atm,O2, & CH4,H2,tcell,ldark(i,j),con,avgcnc,avit(i,j), & hddmjac3,nddmsp,sddm,NOXID,avgsen, & (lddm.AND.lddmcalc(igrd)), & nirrrxn,rrxn_irr,ltrapirr) elseif (idmech.eq.4) then call ebisolv(ebirxn4,ebirate4,hr_hox4,hr_nox4, & hr_nxy4,hr_pan4,hr_xo4, & dtchem,water(i,j,k),atm,O2, & CH4,H2,tcell,ldark(i,j),con,avgcnc,avit(i,j), & hddmjac4,nddmsp,sddm,NOXID,avgsen, & (lddm.AND.lddmcalc(igrd)), & nirrrxn,rrxn_irr,ltrapirr) elseif (idmech.eq.5) then call ebisolv(ebirxn5,ebirate5,hr_hox5,hr_nox5, & hr_nxy5,hr_pan5,hr_xo_dum, & dtchem,water(i,j,k),atm,O2, & CH4,H2,tcell,ldark(i,j),con,avgcnc,avit(i,j), & hddmjac5,nddmsp,sddm,NOXID,avgsen, & (lddm.AND.lddmcalc(igrd)), & nirrrxn,rrxn_irr,ltrapirr) elseif (idmech.eq.6) then call ebisolv(ebirxn6,ebirate6,hr_hox6,hr_nox6, & hr_nxy6,hr_pan6,hr_xo_dum, & dtchem,water(i,j,k),atm,O2, & CH4,H2,tcell,ldark(i,j),con,avgcnc,avit(i,j), & hddmjac6,nddmsp,sddm,NOXID,avgsen, & (lddm.AND.lddmcalc(igrd)), & nirrrxn,rrxn_irr,ltrapirr) endif c elseif ( idsolv .EQ. IDLSOD ) then if (idmech.eq.2) then call lscall(lsrate2,dtchem,water(i,j,k),atm, & O2,CH4,H2,con,avgcnc,nirrrxn, & rrxn_irr,ierr,ltrapirr,0) elseif (idmech.eq.3) then call lscall(lsrate3,dtchem,water(i,j,k),atm, & O2,CH4,H2,con,avgcnc,nirrrxn, & rrxn_irr,ierr,ltrapirr,0) elseif (idmech.eq.4) then call lscall(lsrate4,dtchem,water(i,j,k),atm, & O2,CH4,H2,con,avgcnc,nirrrxn, & rrxn_irr,ierr,ltrapirr,0) elseif (idmech.eq.5) then call lscall(lsrate5,dtchem,water(i,j,k),atm, & O2,CH4,H2,con,avgcnc,nirrrxn, & rrxn_irr,ierr,ltrapirr,0) elseif (idmech.eq.6) then call lscall(lsrate6,dtchem,water(i,j,k),atm, & O2,CH4,H2,con,avgcnc,nirrrxn, & rrxn_irr,ierr,ltrapirr,0) endif endif c c======================== HDDM Begin ==================================== c if ( lhddm .AND. lddmcalc(igrd) ) then if (idmech.eq.2) then call hddmchem(hddmjac2,ebirxn2,ebirate2,ddmrate2, & nreact,ngas,ngas, & nddmsp,nrateddm,ntermddm,nhddm, & iprate,ipterm,wfterm,iphddm, & dtchem,water(i,j,k),atm,O2,CH4,H2, & rk,avgcnc,sddm,ierr) elseif (idmech.eq.3) then call hddmchem(hddmjac3,ebirxn3,ebirate3,ddmrate3, & nreact,ngas,ngas, & nddmsp,nrateddm,ntermddm,nhddm, & iprate,ipterm,wfterm,iphddm, & dtchem,water(i,j,k),atm,O2,CH4,H2, & rk,avgcnc,sddm,ierr) elseif (idmech.eq.4) then call hddmchem(hddmjac4,ebirxn4,ebirate4,ddmrate4, & nreact,ngas,ngas, & nddmsp,nrateddm,ntermddm,nhddm, & iprate,ipterm,wfterm,iphddm, & dtchem,water(i,j,k),atm,O2,CH4,H2, & rk,avgcnc,sddm,ierr) elseif (idmech.eq.5) then call hddmchem(hddmjac5,ebirxn5,ebirate5,ddmrate5, & nreact,ngas,ngas, & nddmsp,nrateddm,ntermddm,nhddm, & iprate,ipterm,wfterm,iphddm, & dtchem,water(i,j,k),atm,O2,CH4,H2, & rk,avgcnc,sddm,ierr) elseif (idmech.eq.6) then call hddmchem(hddmjac6,ebirxn6,ebirate6,ddmrate6, & nreact,ngas,ngas, & nddmsp,nrateddm,ntermddm,nhddm, & iprate,ipterm,wfterm,iphddm, & dtchem,water(i,j,k),atm,O2,CH4,H2, & rk,avgcnc,sddm,ierr) endif if (ierr.ne.0) then write(iout,'(//,a)') 'ERROR in CHEMDRIV:' write(iout,*) & 'Zero determinant in SGEFA in HDDMCHEM at ',ierr write(iout,*) 'igrd, i, j, k = ', igrd,i,j,k call camxerr() endif endif c c======================== HDDM End ======================================= c c-----Perform gas-phase chemistry for SOA precursors c if( aero_flag .GT. 0 ) then avgox(1) = amax1(avgcnc(kO) , bdnl(kO) ) avgox(2) = amax1(avgcnc(kOH) , bdnl(kOH) ) avgox(3) = amax1(avgcnc(kO3) , bdnl(kO3) ) avgox(4) = amax1(avgcnc(kNO3), bdnl(kNO3)) avgox(5) = amax1(avgcnc(kNO) , bdnl(kNO) ) avgox(6) = amax1(avgcnc(kHO2), bdnl(kHO2)) if ( lvbs ) then con_CG( 1) = con(kVAS1) con_CG( 2) = con(kVAS2) con_CG( 3) = con(kVAS3) con_CG( 4) = con(kVAS4) con_CG( 5) = con(kVBS1) con_CG( 6) = con(kVBS2) con_CG( 7) = con(kVBS3) con_CG( 8) = con(kVBS4) con_CG( 9) = con(kVAP1) con_CG(10) = con(kVAP2) con_CG(11) = con(kVAP3) con_CG(12) = con(kVAP4) con_CG(13) = con(kVCP1) con_CG(14) = con(kVCP2) con_CG(15) = con(kVCP3) con_CG(16) = con(kVCP4) con_CG(17) = con(kVFP1) con_CG(18) = con(kVFP2) con_CG(19) = con(kVFP3) con_CG(20) = con(kVFP4) call vbschem(iout,tcell,pcell,dtchem, & avgox,con_prec(:NPREC_V), & con_CG(:NVBS*NSET),con_NV(:NSET)) if (idmech.ne.5 .and. kSQT.lt.nspec+1) & con(kSQT) = amax1( con_prec( 6), bdnl(kSQT) ) con(kIVOG) = amax1( con_prec( 7), bdnl(kIVOG) ) con(kIVOD) = amax1( con_prec( 8), bdnl(kIVOD) ) con(kIVOA) = amax1( con_prec( 9), bdnl(kIVOA) ) con(kIVOB) = amax1( con_prec(10), bdnl(kIVOB) ) con(kVAS1) = con_CG( 1) con(kVAS2) = con_CG( 2) con(kVAS3) = con_CG( 3) con(kVAS4) = con_CG( 4) con(kVBS1) = con_CG( 5) con(kVBS2) = con_CG( 6) con(kVBS3) = con_CG( 7) con(kVBS4) = con_CG( 8) con(kVAP1) = con_CG( 9) con(kVAP2) = con_CG(10) con(kVAP3) = con_CG(11) con(kVAP4) = con_CG(12) con(kVCP1) = con_CG(13) con(kVCP2) = con_CG(14) con(kVCP3) = con_CG(15) con(kVCP4) = con_CG(16) con(kVFP1) = con_CG(17) con(kVFP2) = con_CG(18) con(kVFP3) = con_CG(19) con(kVFP4) = con_CG(20) else con_CG(1) = con(kCG1) con_CG(2) = con(kCG2) con_CG(3) = con(kCG3) con_CG(4) = con(kCG4) c c======================== DDM Begin ======================= c if ( lddm ) then do is = 1, nddmsp sen_CG( is) = sddm(iptddm(kCG1)-1+is) sen_CG( nddmsp+is) = sddm(iptddm(kCG2)-1+is) sen_CG(2*nddmsp+is) = sddm(iptddm(kCG3)-1+is) sen_CG(3*nddmsp+is) = sddm(iptddm(kCG4)-1+is) enddo endif c c======================== DDM End ======================= c call soachem(iout,tcell,pcell,dtchem, & avgox,con_prec(:NPREC_S), & con_CG(:NSOAP),con_NV(:NNVOL),yfac, & nddmsp,avgsen,sen_prec,sen_CG,sen_NV,lddm) c c======================== DDM Begin ======================= c if ( lddm ) then if (idmech.ne.5 .and. kSQT.lt.nspec+1) then do is = 1, nddmsp sddm(iptddm(kSQT)-1+is) = sen_prec(6*nddmsp+is) enddo endif do is = 1, nddmsp sddm(iptddm(kIVOA)-1+is) = sen_prec(3*nddmsp+is) sddm(iptddm(kCG1)-1+is) = sen_CG( is) sddm(iptddm(kCG2)-1+is) = sen_CG( nddmsp+is) sddm(iptddm(kCG3)-1+is) = sen_CG(2*nddmsp+is) sddm(iptddm(kCG4)-1+is) = sen_CG(3*nddmsp+is) enddo endif c c======================== DDM End ======================= c if (idmech.ne.5 .and. kSQT.lt.nspec+1) & con(kSQT) = amax1( con_prec(7), bdnl(kSQT) ) con(kIVOA) = amax1( con_prec(4), bdnl(kIVOA) ) con(kCG1) = con_CG(1) con(kCG2) = con_CG(2) con(kCG3) = con_CG(3) con(kCG4) = con_CG(4) endif endif c c======================== Source Apportion Begin ======================= c c --- subtract the current values of the species from the before values c to get the delta values --- c if( ltrace ) then do ispc=1,nrad delcon(1,ispc) = delcon(1,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac delcon(4,ispc) = delcon(4,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac if( aeropt .NE. 'CF' .OR. aero_flag .EQ. 0 ) then delcon(2,ispc) = delcon(2,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac endif enddo do ispc=nrad+1,ngas delcon(1,ispc) = delcon(1,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac delcon(4,ispc) = delcon(4,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac if( aeropt .NE. 'CF' .OR. aero_flag .EQ. 0 ) then delcon(2,ispc) = delcon(2,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac endif enddo do ispc=ngas+1,nspec delcon(1,ispc) = delcon(1,ispc) + & AMAX1(bdnl(ispc),con(ispc)) delcon(4,ispc) = delcon(4,ispc) + & AMAX1(bdnl(ispc),con(ispc)) if( aeropt .NE. 'CF' .OR. aero_flag .EQ. 0 ) then delcon(2,ispc) = delcon(2,ispc) + & AMAX1(bdnl(ispc),con(ispc)) endif enddo endif c c======================== Source Apportion End ======================= c c======================= Process Analysis Begin ======================== c if( ldoipr ) then c c --- Chemistry change in umol/m3 units (except for HG chemistry) ---- c do is=1,nrad cipr(IPR_CHEM, ipa_idx, is) = & cipr(IPR_CHEM, ipa_idx, is) + & con(is)-conc(i,j,k,is) enddo do is=nrad+1,ngas cipr(IPR_CHEM, ipa_idx, is) = & cipr(IPR_CHEM, ipa_idx, is) + & con(is)*convfac-conc(i,j,k,is) enddo endif c c======================== Process Analysis End ========================= c c c========================= DDM Begin =================================== c c --- put the changed values back into the gridded array --- c if( (lddm .OR. lhddm) .AND. lddmcalc(igrd) ) then if( lhddm ) & call hddmadj(.TRUE.,nddmsp,con,sddm) call loaddm(.TRUE.,m1,m2,m3,ntotsp,ptconc(iptrsa), & i,j,k,nddmsp,nspec,ngas,nrad,NHETRXN,sddm, & convfac) endif c c========================= DDM End ===================================== c c c-----Perform aerosol chemistry if it's time c if ( aero_flag.gt.0 ) then cliq = cwc(i,j,k) if (lcig(igrd)) cliq = amax1(cliq,cigwtr(i,j,k)) if (tcell.lt.273.) & cliq = amax1(0.,cliq*(tcell - tamin)/(273. - tamin)) if( aeropt.eq.'CF' ) then c c======================== DDM Begin ==================================== c c ---- load the sensitivities for this cell into 2-D array --- c if ( lddm ) then call loaddmpm(.FALSE., & m1,m2,m3,ntotsp,ptconc(iptrsa), & i,j,k,nddmsp,NSDDMPM,sddm,convfac) endif c c======================== DDM End ======================================= c icig = 0 if (lcig(igrd)) icig = 2 call aerochem_cf(water(i,j,k),tcell,pcell,cliq,cldph, & aerph,ajno2(i,j,k),con_NV, & con,convfac,dtchem,aero_dt(igrd), & delcon,sen_NV,sddm,nddmsp,NSDDMPM, & ldoipr,ipa_idx,aero_flag,.FALSE.,icig, & ldark(i,j)) c c========================= DDM Begin =================================== c c --- put the changed values back into the gridded array --- c if ( lddm ) then call loaddmpm(.TRUE., & m1,m2,m3,ntotsp,ptconc(iptrsa), & i,j,k,nddmsp,NSDDMPM,sddm,convfac) endif c c========================= DDM End ===================================== c elseif( aeropt.eq.'CMU' ) then call aerochem_cmu(water(i,j,k),tcell,pcell,cliq,cldph, & con_NV, & MXSPEC,con,bdnl,convfac, & timec(igrd),dtchem,aero_dt(igrd), & aero_flag,ldark(i,j)) endif endif c c-----Prepare for surface model c if (k.eq.1 .and. lsrfmod) then ij = i + (j-1)*m1 depth = height(i,j,1) press0 = press(i,j,1) - depth*(press(i,j,2) - & press(i,j,1))/height(i,j,2) ucomp = (uwind(i,j,1) + uwind(i-1,j,1))/2. vcomp = (vwind(i,j,1) + vwind(i,j-1,1))/2. wind = sqrt(ucomp**2 + vcomp**2) wind = amax1(0.1,wind) do l = 1,NLU flu(l) = fsurf(i,j,l) enddo do lsm = 1,nsmspc soil(lsm) = solmas(i,j,lsm) vegy(lsm) = vegmas(i,j,lsm) enddo call srfmod(datec(igrd),dtchem,dx(j+j0),dy,cellat(i,j), & cellon(i,j),snow(i,j),snowrat(i,j), & snowfrc(i,j),snowalb(i,j),zenith,ctrns, & ldark(i,j),pwr(i,j,1),tempk(i,j,1),tsurf(i,j), & press(i,j,1),press0,wind,depth,nsmspc, & flu,lrdlai,lai(i,j),soil,vegy,reemit) do lsm = 1,nsmspc solmas(i,j,lsm) = soil(lsm) vegmas(i,j,lsm) = vegy(lsm) reemis(i,j,lsm) = reemit(lsm) enddo endif c c========================= Hg Chemistry Begin ========================== c c-----Call Hg chemistry modules c if (khg0.lt.nspec+1) then c c-----Assign Cl2, Br and BrO concentrations based on height(km), ocean c-----and day/night c do ihg = 1,NHTHAL if (hght .le. hthal(ihg)) then iht = ihg go to 115 endif enddo iht = NHTHAL 115 continue brc = 0. broc = 0. cl2c = 0.0 if (.not. ldark(i,j)) then if (iocean .eq. 1) then brc = brwprof(iht) broc = browprof(iht) cl2c = cl2day(iht) else brc = brlprof(iht) broc = brolprof(iht) endif else cl2c = cl2nite(iht) endif c c-----Gas-phase chemistry c c c-----Put adsorbed Hg(II) into gas-phase Hg(2) for chemistry calculations c-----convert adsorbed Hg(II) to ppm before adding c convhgfac = MWHG*convfac thg2 = con(khg2) + (con(khgiip)+con(khgiipc))/convhgfac thg20 = thg2 ! For process analysis call hggaschem(con(khg0),thg2,con(ko3), & con(kh2o2),con(koh),brc,broc, & tcell,pcell,dt) con(khg0) = max(bdnl(khg0),con(khg0)) thg2 = max(bdnl(khg2),thg2) c c======================= Process Analysis Begin ======================== c if( ldoipr ) then chg0 = con(khg0)*convfac ! store HG0 in umol/m3 chg2 = thg2*convfac ! store total HG2 in umol/m3 cipr(IPR_CHEM, ipa_idx, khg0) = & cipr(IPR_CHEM, ipa_idx, khg0) + & chg0 - conc(i,j,k,khg0) cipr(IPR_CHEM, ipa_idx, khg2) = & cipr(IPR_CHEM, ipa_idx, khg2) + & chg2 - thg20*convfac endif c c======================== Process Analysis End ========================= c c-----Aqueous-phase chemistry (conditional) c cliq = cwc(i,j,k) if (tcell.lt.273.) then cliq = amax1(0.,cliq*(tcell - tamin)/(273. - tamin)) endif if (cliq.ge.cwmin .and. tcell.ge.tamin ) then c c-----Calculate total PM concentrations c henso2 = henry0(kso2) tfso2 = tfact(kso2) heno3 = henry0(ko3) tfo3 = tfact(ko3) pm10 = 0. do is = ngas+1,nspec if (is.ne.kph2o) pm10 = pm10 + con(is) enddo call hgaqschem(con(khg0),thg2,con(ko3), & con(kso2),con(koh),con(kho2), & con(khcl),cl2c,pm10,cldph,cliq, & tcell,pcell,dt,henso2,tfso2, & heno3,tfo3) con(khg0) = max(bdnl(khg0),con(khg0)) thg2 = max(bdnl(khg2),thg2) c c======================= Process Analysis Begin ======================== c if( ldoipr ) then cipr(IPR_AQCHEM, ipa_idx, khg0) = & cipr(IPR_AQCHEM, ipa_idx, khg0) + & con(khg0)*convfac - chg0 cipr(IPR_AQCHEM, ipa_idx, khg2) = & cipr(IPR_AQCHEM, ipa_idx, khg2) + & thg2*convfac - chg2 endif c c======================== Process Analysis End ========================= c endif c c-----Adsorption of Hg(2) to particles c-----Calculate volume concentrations of fine & coarse primary particles vfine = 0. if(kpec < nspec+1)vfine = vfine + con(kpec)/roprt(kpec) if(kpoa < nspec+1)vfine = vfine + con(kpoa)/roprt(kpoa) if(kfprm < nspec+1)vfine = vfine + con(kfprm)/roprt(kfprm) c-----Add sea-salt to vfine if(kna == nspec+1)then vseas = nacl/RHOSEAS else ! If Na is a modeled species, then PCl is too cna = con(kna)/MWNA cpcl = con(kpcl)/MWCL if(cna < cpcl) then vseas = cna*MWNACL/roprt(kna) else vseas = cpcl*MWNACL/roprt(kpcl) end if end if vfine = vfine + vseas dfine = sqrt(dcut(kph2o,1)*dcut(kph2o,2)) !H2O is a reqd. species c vcoarse = 0. if (kccrs < nspec+1) then vcoarse = vcoarse + con(kccrs)/roprt(kccrs) dcoarse = sqrt(dcut(kccrs,1)*dcut(kccrs,2)) end if if (kcprm < nspec+1) then vcoarse = vcoarse + con(kcprm)/roprt(kcprm) dcoarse = sqrt(dcut(kcprm,1)*dcut(kcprm,2)) end if c c-----calculate number conc, surface area conc of fine and coarse primary PM if (vfine > 0.)then fseas = vseas/vfine fseas = MAX(0.0,MIN(1.0,fseas)) safine = 6.0 * vfine / dfine ! Units: m2/m3 else fseas = 1.0 safine = 0.0 end if if (vcoarse > 0.)then sacoarse = 6.0 * vcoarse / dcoarse ! Units: m2/m3 else sacoarse = 0.0 end if c if (safine > 0. .or. sacoarse > 0.) then chgiip = 0. chgiipc = 0. chg2 = thg2 * convhgfac ! in ug/m3 call hgadsorb(tcell,chg2,chgiip,chgiipc,fseas,safine, & sacoarse) c c-----assign total hg2 back to individual species con(khg2) = chg2/convhgfac ! back to ppm con(khg2) = max(bdnl(khg2),con(khg2)) con(khgiip) = max(bdnl(khgiip),chgiip) con(khgiipc) = max(bdnl(khgiipc),chgiipc) else con(khg2) = thg2 con(khgiip) = bdnl(khgiip) con(khgiipc) = bdnl(khgiipc) endif endif c c========================== Hg Chemistry End =========================== c c c======================== Source Apportion Begin ======================= c c --- subtract the current values of the species from the before values c to get the delta values --- c if( ltrace ) then do ispc=1,nrad cnew(ispc) = con(ispc)*convfac delcon(3,ispc) = delcon(3,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac enddo do ispc=nrad+1,ngas cnew(ispc) = con(ispc)*convfac delcon(3,ispc) = delcon(3,ispc) + & AMAX1(bdnl(ispc),con(ispc)) * convfac enddo do ispc=ngas+1,nspec cnew(ispc) = con(ispc) delcon(3,ispc) = delcon(3,ispc) + & AMAX1(bdnl(ispc),con(ispc)) enddo o3new = AMAX1(bdnl(ko3),con(ko3)) c if( ldark(i,j) ) then desto3 = delcon(2,ko3) prdo3n = 0.0 prdo3v = 0.0 else if( lozone ) then call o3prdsa(rrxn_irr,convfac,delcon(2,ko3), & prdo3n,prdo3v,desto3) if( delcon(2,ko3) .lt. desto3 ) then desto3 = delcon(2,ko3) prdo3n = 0.0 prdo3v = 0.0 endif endif endif c c --- rrxn_irr: ppm -> umol/m3 c if( ltrapirr ) then do irxn=1,nreact rrxn_irr(irxn) = rrxn_irr(irxn) * convfac enddo endif c c --- SA: call routine to do the tracer species "chemistry", just makes c adjustments for production or decay of the regular model species ---- c if( lozone .OR. lnitrate ) then call osatsa(m1,m2,m3,igrd,i,j,k,prdo3n,prdo3v,desto3, & nspec,delcon,convfac,cold,rrxn_irr,dtchem) endif c if( lsulfate .OR. lnitrate .OR. lsoa .OR. lprimary .OR. & lmercury ) then call psatsa(m1,m2,m3,igrd,i,j,k,nspec,delcon, & NPREC_S,yfac) call repartsa(m1,m2,m3,igrd,i,j,k,cnew,delcon) endif endif c c --- RTRAC chemistry c if( tectyp .EQ. RTRAC ) then if (nrtherm .GT. 0 .OR. nrtphot .GT. 0 ) then if( ipa_cel(i+i0,j+j0,k) .GT. 0 ) then irt_cel = ipa_cel(i+i0,j+j0,k) else irt_cel = -9 endif call chemrt(m1,m2,m3,igrd,i,j,k,i0,j0,pcell,tcell,cold,cnew, & avgcnc(koh),avgcnc(ko3),avgcnc(kno3),dtchem,irt_cel) endif if( lparttn ) & call partitionrt(m1,m2,m3,nspec,igrd,i,j,k,cnew) if( k.eq.1 .and. lsrfmodrt ) & call srfmodrt(m1,m2,igrd,i,j,dtchem,zenith,ctrns, & ldark(i,j)) c c --- RTCMC chemistry c else if( tectyp .EQ. RTCMC ) then c c --- set rate constants c if( ktype .EQ. 1 ) then call krtc(tcell,pcell) else call ksci(tcell,pcell) endif if( njschm .GT. 0 ) then do nj = 1,njschm srkrtc(ijschm(nj)) = 0. if (ijpntr(nj).gt.0) then srkrtc(ijschm(nj)) = rk(ijpntr(nj)) else do nz = nzschm,1,-1 if( zenith .LE. zenschm(nz) ) & srkrtc(ijschm(nj)) = rjschm(nz,nj)*ctrns enddo endif rkrtc(ijschm(nj)) = DBLE( srkrtc(ijschm(nj)) ) enddo endif c c --- do chemistry c do l = 1,ntotsp idx = DBLE(iptrsa-1) + DBLE(i) + & DBLE(m1)*DBLE(j-1) + DBLE(m1)*DBLE(m2)*DBLE(k-1) + & DBLE(m1)*DBLE(m2)*DBLE(m3)*DBLE(l-1) rtcon(l) = ptconc(idx)/convfac enddo if( isolv .EQ. 1) then call dlsdrv(dtchem,water(i,j,k),atm,O2,CH4,H2, & avgcnc,rtcon,ierr,.false.) elseif( isolv .EQ. 2) then call slsdrv(dtchem,water(i,j,k),atm,O2,CH4,H2, & avgcnc,rtcon,ierr,.false.) elseif( isolv .EQ. 3) then call rbkdrv(dtchem,water(i,j,k),atm,O2,CH4,H2, & avgcnc,rtcon,ierr,.false.) endif do l = 1,ntotsp idx = DBLE(iptrsa-1) + DBLE(i) + & DBLE(m1)*DBLE(j-1) + DBLE(m1)*DBLE(m2)*DBLE(k-1) + & DBLE(m1)*DBLE(m2)*DBLE(m3)*DBLE(l-1) ptconc(idx) = rtcon(l)*convfac enddo endif c c========================= Source Apportion End ======================== c c c======================= Process Analysis Begin ======================== c c --- Chemical reaction tracking in ppm units c Account for titration reaction in TRAP at night --- c if( ldoirr ) then do irxn=1,nirrrxn cirr(ipa_idx, irxn) = & cirr(ipa_idx, irxn) + rrxn_irr(irxn) enddo if ( titrt .GT. 0.0 ) then if ( idmech.EQ.5 ) then cirr(ipa_idx, 7) = & cirr(ipa_idx, 7) + titrt else cirr(ipa_idx, 3) = & cirr(ipa_idx, 3) + titrt endif endif endif c c --- calculate the chemical process analysis (CPA) variables c time weight CPA variables that are not cumulative by dtfact --- c if( lirr .AND. (l3dflag .OR. k .EQ. 1) ) then nn = 0 dtfact = dt/dtout/60. if( idmech.EQ.2 ) then call cpamech2(rrxn_irr,rk,dtfact,MXREACT,patmp,ntotsp, & nn,ldark(i,j)) elseif( idmech.EQ.3 ) then call cpamech3(rrxn_irr,rk,dtfact,MXREACT,patmp,ntotsp, & nn,ldark(i,j)) elseif( idmech.EQ.4 ) then call cpamech4(rrxn_irr,rk,dtfact,MXREACT,patmp,ntotsp, & nn,ldark(i,j)) elseif( idmech.EQ.5 ) then call cpamech5(rrxn_irr,rk,dtfact,MXREACT,patmp,ntotsp, & nn,ldark(i,j)) elseif( idmech.EQ.6 ) then call cpamech6(rrxn_irr,rk,dtfact,MXREACT,patmp,ntotsp, & nn,ldark(i,j)) else write(iout,'(//,a)') 'ERROR in CHEMDRIV:' write(iout,'(2A,I2)')'Chemical Process Analysis (CPA)', & ' does not support chemical mechanism: ',idmech write(iout,*) & 'Please choose another mechanism and try again.' call camxerr() endif c c --- save the cloud adjustment to photolysis rates --- c nn = nn + 1 ptname(nn) = 'J_CLDADJ' patmp(nn) = ctrns * dtfact c c --- save the aerosol pH --- c nn = nn + 1 ptname(nn) = 'AER_PH' patmp(nn) = aerph * dtfact c c --- save radical concentrations unless accumulating throught run --- c if( .NOT. lcpacum) & call cparad(con, ngas, patmp, ntotsp, nn, dtfact) c c --- add CPA values to gridded array in native units --- c do l=1,ntotsp idx = DBLE(iptrsa-1) + DBLE(i) + & DBLE(m1)*DBLE(j-1) + DBLE(m1)*DBLE(m2)*DBLE(k-1) + & DBLE(m1)*DBLE(m2)*DBLE(nlay)*DBLE(l-1) ptconc(idx) = ptconc(idx) + patmp(l) enddo endif c c======================== Process Analysis End ========================= c c-----Pass CON back to concentration, convert gases from ppm to umol/m3 c 88 continue cph(i,j,k) = cldph aph(i,j,k) = aerph do is=1,nrad conc(i,j,k,is) = con(is) enddo do is=nrad+1,ngas conc(i,j,k,is) = con(is)*convfac enddo if (ngas.lt.nspec) then do is=ngas+1,nspec conc(i,j,k,is) = amax1(con(is),bdnl(is)) enddo endif c 89 continue 90 continue c c$omp end parallel c sumit(k) = 1. c nitsum(k) = 1 c do i = ia,iz c do j = ja,jz c sumit(k) = sumit(k) + avit(i,j) c if (idfin(i,j).le.igrd ) nitsum(k) = nitsum(k) + 1 c enddo c enddo c 91 continue c write(idiag,'(a,2i5,2f8.1)') c & ('EBI ',igrd,k,timec(igrd),sumit(k)/nitsum(k),k=1,nlay) c c-----If aerosol chemistry was called c - reset aero_dt & ajno2 c - increase time_aero (& date_aero) by dt_aero (or multiple of dt_aero) c if (aero_flag.EQ.2) then aero_dt(igrd) = 0.0 ajno2 = 0.0 92 continue call uptime(time_aero(igrd),date_aero(igrd),60.*dt_aero) if ( time_aero(igrd) .le. timec(igrd) .and. & date_aero(igrd) .eq. datec(igrd) ) goto 92 endif c return end