program advect ! box model version of regional model, hank implicit none include "chem_params.inc" !number of species, indices of sp. include 'cchem.com' !character list of species include 'ichem.com' !mapping indices of reactions include 'rchem.com' !MW and stoich. coef include 'irate.com' !type of reaction (hv, Arrhen, etc) include 'rrate.com' !variables for calc. Arrhen and TROE rxns include 'amap.com' !maps species to/from gas and aqueous arrays include 'hvchem_jvalchgs.com' c **** constants *** common /dltprms/dtchem, dt_dbg real dtchem ! the chemistry timestep real dt_dbg ! how often to do a debug print common /constant/pi,r,gr,cp real pi ! pi = PI real r ! gas constant for dry air = RD real gr ! acceleration of gravity = GZERO real cp ! specific heat of dry air at constant pressure = CPD real na ! Avogadro's number (molec/mol) real mwair ! MW of air (g/mol) real rate(nrmx) ! reaction rate constants of gas rxns real ratea(nrmx) ! reaction rate constants of cw rxns real qca, qva ! cloud water kg/kg or g/g and water vapor g/g real lww ! cloud water cm3/cm3 real ph ! pH of drops integer ncpts ! = 1 if cloud exists integer npts ! = 1 for clear sky common/cldpar/qca, lww, ph, ncpts, npts real press, rho, ttemp, rhoair, zalt c ktc is the rate constant of diffusion from gas to liquid c aheff is henry's law coefficients for aqueous species c phrat is phase ratio = Ca/Cg real ktc(naq), aheff(naq) real phrat(naq) real rhow ! water density c cloud water and rain data integer nctim, incw, np integer icld parameter(np=2000) real cwtim(np), tair_in(np), rho_in(np), pres_in(np), z(np), w(np) real qv_in(np), qc_in(np), qi_in(np), qr_in(np), qs_in(np), _ qg_in(np), secofday(np) c--mcb c **** meteorological variables *** real q ! specific humidity (kg/kg) real t ! temperature (K) real spec(nspec) ! gas species concentration _ ,aspec(naq) ! cloud drop species _ ,aqspc(naq) ! cloud drop species in Molar units c *** time: real starttime ! start time of integration -- read in Chemset real endtime ! end time of integration real tim ! the model tim starting from 0 to timf (in seconds) real timf ! the final time of the model from beginning (in seconds) real timt ! time value used for photolysis rate calculations real startcld ! time to start cloud c indexes integer i,ii,j,k,l,ir,isp,i2,n,ip c unit numbers for input/output: integer & inrt, ! inrt: react unit number & inph ! inph: photolysis unit number real cinit(nspec), cinit1(nspec) real dens, tsteps, writtim integer icount, nt logical iprnt c photolysis rate variables (+ common block included above) integer mr1, jtim, jtimp, it, irphot(30) real tday, tj, pslope, tdif, dtj, otim real delttj character*4 ci(0:nspec) character*7 aname character*20 finame external & chemset, xinit, getph, rates, glrates, chem c----------------------------------------------------------------------- c Initialization c c map gas and aqueous phase indices c do i=1,nspec gasmap(i) = 0 end do c++mcb -- map aqueous species to gas species aqmap(0) = 0 aqmap(ko3a) = ko3 aqmap(kh2o2a) = kh2o2 aqmap(kho2a) = kho2 aqmap(koha) = koh aqmap(kso2a) = kso2 aqmap(kso4a) = kso4 c++mcb -- map gas species to aqueous species gasmap(ko3) = ko3a gasmap(kh2o2) = kh2o2a gasmap(kho2) = kho2a gasmap(koh) = koha gasmap(kso2) = kso2a gasmap(kso4) = kso4a c----------------------------------------------------------------------- c initialize the species mixing ratios to zero call xinit(spec,nspec,0.) call xinit(aspec,naq,0.) call xinit(aqspc,naq,0.) c define density of water (kg/m3) rhow = 1000. c assign unit numbers for files (remnant of HANK model) inrt = 15 inph = 11 incw = 14 icount = 0 ! for write output time of mixing ratios c----------------------------------------------------------------------- c read the input file and set up each species P and L terms call chemset(inrt, cinit, nspec, starttime, endtime, t, rhoair) c----------------------------------------------------------------------- write(*,*) write(*,*) '-----------------------------------------------------' write(*,*) 'after chemset: nreact, naqr ', nreact, naqr write(*,*) ' aqmap gasmap ' do isp=1,naq write(*,'(a,2i6)') cx(aqmap(isp)), aqmap(isp), gasmap(aqmap(isp)) end do c initialize constants na = 6.0221e23 mwair = 28.97 pi=2.*asin(1.) gr=9.81 r=287. cp=1004. c c initialize time parameters c dtchem = 10. dt_dbg = 60. tim = starttime timf = endtime writtim = dtchem tsteps = (timf-tim)/writtim nt = tsteps + 1 c Calculate air density, water vapor mixing ratio dens = rhoair * na/mwair *1.e-3 !air density [molecules/cm3] !! q = cinit(kh2o) *18./28.97 !water vapor mixing ratio (cinit is mol/mol; q is g/g) q = 12.3e-3 *18./28.97 !! print*,'dens= ',dens, rhoair, 'H2O = ',cinit(kh2o), q print*,'dens= ',dens, rhoair, 'H2O = ', q c Assign initial vmr (cinit read in chemset) to gas species array (spec) do i=1,nspec spec(i) = cinit(i) enddo c----------------------------------------------------------------------- c Open output file and print initial concentration in file open(unit=17,file='thermo.out',form='formatted') open(unit=18,file='gas.out',form='formatted') open(unit=19,file='cw.out',form='formatted') open(unit=20,file='total.out',form='formatted') open(unit=21,file='aqcw.out',form='formatted') ! aqcw.out has Molar units instead of vmr 1020 FORMAT(100(E13.5,' ')) c write headers in output files write(17,'(a,4x,10(3x,a11,2x))') ' Time', ' z_(km) ', _ ' T_(K) ',' P_(hPa) ', 'rho_(kg/m3)', 'qv_(kg/kg) ', _ 'qc_(kg/kg) ',' pH ' write(18,'(a,4x,100(3x,a10,2x))') ' Time', (cx(j), j=1,nspec) write(19,'(a,4x,100(3x,a10,2x))') ' Time', _ (cx(aqmap(j)), j=1,naq) write(20,'(a,4x,100(3x,a10,2x))') ' Time', _ (cx(aqmap(j)), j=1,naq) write(21,'(a,4x,100(3x,a10,2x))') ' Time', _ (cx(aqmap(j)), j=1,naq) c write initial concentrations into output files write(18,1020) tim, (spec(k) ,k=1,nspec) write(19,1020) tim, (aspec(k) ,k=1,naq) write(20,1020) tim, (spec(aqmap(k)) + aspec(k) ,k=1,naq) ! aqcw.out initial time is written after cloud water data read in c----------------------------------------------------------------------- c count time steps icount=icount+1 write(*,*) ' icount after initialization ', icount c++mcb -- read in cloud water data from external file c startcld = 0.*60. ! start cloud after 30 minutes of running gas chemistry finame='cwrn' ! input file (cw=cloud water, rn=rain but rain is not addressed in this model version) open(file=finame,unit=incw) read(incw,*) ! skip header i = 0 59 continue i = i+1 read(incw,*,err=2002,end=23) secofday(i), qc_in(i) c use secofday to get altitude and cwtim cwtim(i) = secofday(i) z(1) = 1.5 ! altitude of cloud point (prescribed, used for j-values) z(i) = z(1) + cwtim(i) * 0.0 ! 0.0 could be set to an updraft speed to represent a moving parcel cwtim(i) = cwtim(i) + startcld go to 59 ! keep reading file until each reaches its end 23 continue nctim = i-1 write(*,*) starttime + startcld cwtim(1) = starttime ! set first time in CWRN file to starttime to be sure LWC=0 from start do i = 1,5 write(*,*) "firstcld five times ", i, z(i), cwtim(i), qc_in(i) end do do i = nctim-2,nctim write(*,*) " lastcld three times ", i, z(i), cwtim(i), qc_in(i) end do c--mcb c Initial guess at pH ph = 4.5 c Calculate air density, water vapor mixing ratio dens = rhoair * na/mwair *1.e-3 !air density !! spec(kh2o) = q * mwair/18. !! print*,'dens= ',dens, rhoair, 'H2O = ',spec(1),' mol/mol ', !! _ spec(1)*dens, ' molec/cm3 ', q,' kg/kg ', t c icld is a flag = 0 (no cloud) or 1 (cloud exists) c icld allows aqueous concentrations (aqspc) to be 0 if there's no cloud icld = nint(qc_in(1)/max(qc_in(1), 1.e-20)) write(*,*) 'ICLD = ', icld, qc_in(1) do k=1,naq aqspc(k) = aspec(k) * 1.e6/(mwair*max(qc_in(1), 1.e-20)) * icld ! aqspc [mol X/L H2O] and aspec [vmr] end do write(21,1020) tim, (aqspc(k) ,k=1,naq) c**********************************************************************c c Time loop: c**********************************************************************c Time_loop: do write(*,*) 'Time: ', tim, rhoair, q, t, dens, aspec(kso2a) c set cloud water mixing ratios timt = tim do i=1,nctim-1 if(timt .ge. cwtim(i) .and. timt .lt. cwtim(i+1)) then qca = qc_in(i) + (timt-cwtim(i)) * (qc_in(i+1)-qc_in(i)) / _ (cwtim(i+1)-cwtim(i)) dens = rhoair * na/mwair *1.e-3 !air density, molecules/cm3 press= rhoair*287.*t/100. !use ideal gas law instead of interpolation to be sure rho, T, p are consistent zalt = z(i) + (timt-cwtim(i)) * (z(i+1)-z(i)) / ! used for choosing j-value altitude in subr. rates _ (cwtim(i+1)-cwtim(i)) !! q = cinit(kh2o) *18./mwair !water vapor mixing ratio go to 24 else qca = 0. endif end do 24 continue !! spec(kh2o) = q * mwair/18. ! mol H2O/mol air !! write(*,*) tim, zalt, qca, t, rhoair, q, ' kg/kg ', !! _ spec(kh2o), ' mol/mol set qca, q H2O ' if(tim .eq. starttime) then write(17,1030) tim, zalt, t, press, rhoair, q, qca, ph end if c----------------------------------------------------------------------- c adjust spec for new air density and put into molecules/cm3 units do isp=1,nspec spec(isp) = spec(isp)*dens end do do isp=1,naq aspec(isp) = aspec(isp)*dens end do c++mcb -- Find the cloudy grid points ncpts = 0 npts = 1 if(qca .gt. 1.e-12) then lww = rhoair * qca/(rhow*1000.) ncpts = 1 npts = 0 write(*,'(2e12.5)') qca, q write(*,'(1p,3e12.5,2i5)') lww, rhoair, rhow, ncpts, npts c======================================================================= c for species that are considered only in the aqueous phase ! aspec(kco3m) = aspec(kco3m) + spec(kco3mg) ! spec(kco3mg) = 0. aspec(kso4a) = aspec(kso4a) + spec(kso4) spec(kso4) = 0. else ncpts = 0 npts = 1 do isp=1,naq spec(aqmap(isp)) = spec(aqmap(isp)) + aspec(isp) aspec(isp) = 0. end do endif c----------------------------------------------------------------------- c get the reaction rates c----------------------------------------------------------------------- ! get pH of cloud water (qca) ! call getph(t, rhoair, qca, aspec, ph) ph = 4.5 ! get rate coefficients for all reactions call rates(t, rhoair, q, qca, ph, rate, ratea, tim, zalt) ! get rate coefficients for gas-liquid transfer call glrates(t, rhoair, ph, ktc, lww, aheff, phrat) c----------------------------------------------------------------------- c solve the chemistry call chem(rate,t,q,spec, ratea, ktc,aheff,phrat,aspec, tim) c----------------------------------------------------------------------- c adjust spec back to vmr do isp=1,nspec spec(isp) = spec(isp)/dens end do do isp=1,naq aspec(isp) = aspec(isp)/dens end do c increment the time tim = tim+dtchem print*,'elapsed time: ',tim c c write to the output file otim=real(icount)*writtim iprnt=.false. print*,'elapsed time: ',tim, otim, icount+1, nt if(tim.ge.otim.and.(icount+1).le.nt) then iprnt=.true. icount=icount+1 write(17,1030) tim, zalt, t, press, rhoair, q, qca, ph write(18,1020) tim,((spec(k)),k=1,nspec) write(19,1020) tim, (aspec(k) ,k=1,naq) write(20,1020) tim, ((spec(aqmap(k)) + aspec(k)) ,k=1,naq) icld = nint(qca/max(qca, 1.e-20)) do k=1,naq aqspc(k) = aspec(k) * 1.e6/(mwair*max(qca, 1.e-20)) * icld end do write(21,1020) tim, (aqspc(k) ,k=1,naq) endif if(tim >= timf) then exit Time_loop endif enddo Time_loop 1030 FORMAT(5(F13.5,' '),2(E13.5,' '),F13.5 ) 1001 continue c write(17,1030) tim, zalt, t, press, rhoair, qva, qca, ph print*,'times written ',nt,icount stop 'AQCHEM Finished' 2000 continue write(6,1140) finame stop 2010 continue write(6,1150) finame stop 2002 continue write(6,1140) finame stop 1140 FORMAT ('Error in Reading Input File:',A20) 1150 FORMAT ('Unexpected End-Of-File in Input File:',A20) end program advect subroutine chemset(inrt, cinit, mspec, starttime, endtime, _ tair, rhoair) implicit none include 'chem_params.inc' ! ! this routine sets the reactions and writes a map of the ! reactions into icmap. icmap is labeled by the species 1....n. for each species ! the ordering is as follows: ! ! # of sources ! # of sinks ! then for each source reaction the following five numbers are read ! until all the source reactions are accounted for ! reaction number ! label of stoich coefficient in stc ! 1 reactant ! 2 reactant ! 3 reactant ! then for each sink reaction the following five numbers are read ! until all the sink reactions are accounted for ! reaction number ! label of stoich coefficient in stc ! 1 reactant ! 2 reactant ! 3 reactant ! 4 reactant ! ! variables: ! ! c common block storage c c cw c cx character list of species c stc matrix of standard stoichometry coefficients c icmap reaction matrix for each species c nreact total number of reactions c a temperature independent coeff for arrehenius reactions c s temperature dependent coeff for arrehenius reactions c stc stoichmetry coefficients c nspec number of chemically interactive species c c c temporary storage c c isrc temporary storage of map of sinks c isnk temporary storage of map of sources c iblank blank character c c2 index for special reaction c molec reaction molecules c sc stoichemetry coefficients c mr1 number of photolysis reactions c mr2 number of arrhenius reactions c mr3 number of m reactions c mr4 number of troe reactions c mr5 number of special reactions c iwhich index indentifier for special reactions c inrt input file numbers c finame input file name c itype1 - photolysis - already assigned at read c itype2 - simple arrhenius c itype3 - a + m or a + b + m, possibly also arrhenius c itype4 - troe - already assigned at read c itype5 - special functions for rate constants c*** troe reaction data, molecule, cm3, sec units c ak0300 - zero pressure 300k rate constant, third order c an - temperature exponent for zero pressure rate constant c aki300 - high pressure 300k rate constant, second order c am - temperature exponent for high pressure rate constant c bas - base of exponentiation, 0.6 for most c aequil - pre-exponential of equilibrium constant c tequil - activation temperature of equilibrium constant c** c input/output integer inrt, mspec real cinit(mspec) real starttime, endtime real tair, rhoair c temporary storage integer mr1,mr2,mr3,mr4,mr5,mr6,mr7 integer ntx,ncx,nx c integer c & isrc(5*nrps), c & isnk(5*nrps) character*20 finame character*2 c2(nrmx) character*4 molec(nrmx,7) character*4 iblank real sc(nrmx,7) c indexes integer ix,i,ir,j,k,ik,n,jr,it,iy,iz integer nsnks,nsrcs, nsnksa,nsrcsa integer ip,ir1,is1,ip1,ip2,ip3,ip4 real csnk character*80 astrg integer irphot(30), lr1 c chemical storage block include 'cchem.com' include 'ichem.com' include 'rchem.com' include 'irate.com' include 'amap.com' include 'rrate.com' include 'hvchem_jvalchgs.com' external opndat, maprxns DATA IBLANK /' '/ do i=0,nspec cx(i)=iblank enddo do i=0,naq cxa(i)=iblank enddo c weight the molecular weights c molecules/mole c cw=grams/mole c 1.e3 is conversion to molecules/kg FINAME = 'react' CALL OPNDAT( INRT, FINAME ) read(inrt,*) nreact,mr1,mr2,mr3,mr4,mr5,mr6,mr7,mzphot print*,'Number of reactions ',nreact print*,'Number of photolysis reactions ',mr1 print*,'Number of arrhenius reactions ',mr2 print*,'Number of m reactions ',mr3 print*,'Number of troe reactions ',mr4 print*,'Number of special reactions ',mr5 print*,'Number of aqueous reactions ',mr6 print*,'Number of aqueous photo reactions ',mr7 print*,'Number of new.jv files for altitude dependence',mzphot c read in initial time of simulation, final time, c temperature (K), and air density (kg/m3) read(inrt,*) starttime, endtime, tair, rhoair print*,'Start/Stop times ', starttime, endtime print*,'Temperature, Air density ', tair, rhoair c read in species names, MWs and initial concentration READ(inrt,*) nx ! total number of chemical species print*,'Number of total species ',nx if(nx.ne.nspec) then print*,'MISMATCH BETWEEN PARAMETER AND REACT FILES' print*,'REACT, TOTAL SPECIES ',nx print*,'PARAM, TOTAL SPECIES ',nspec stop endif print*,' ' print*,'SPECIES ' do i=1,nspec read(inrt,*) cx(i),cw(i), cinit(i) ! mol X / mol air for all species write(6,602) i,cx(i),cw(i), cinit(i) enddo 602 format(i3,1x,a4,1x,f6.2,3x,1p,e13.5) c stoich coefficients: set the standard ones call xinit(stc,(2*nspec),0.) call xinit(sc,(nrmx*7),0.) stc(0)=99. stc(1)=1. stc(2)=2. c read in reaction file print*,' ' do 15 j = 1, nreact read(inrt,1090,err=2000,end=2010) c2(j), | (molec(j,i),i=1,3),(sc(j,i),molec(j,i),i=4,7),a(j),s(j) write(6,1091) j,c2(j), | (molec(j,i),i=1,3),(sc(j,i),molec(j,i),i=4,7),a(j),s(j) 15 continue write(*,*) write(*,*) ' rxns read in ' do 17 i = 1, mr4 read(inrt,1080,err=2000,end=2010) itype4(i), | ak0300(i),an(i),aki300(i),am(i),bas(i), | aequil(i),tequil(i) 17 continue write(*,*)'done reading troe reaction rates ' close(inrt) c----------------------------------------------------------------------c c read in photolysis frequencies if(mr1 .gt. 0) then do iz=1,mzphot if(iz.le.10) then write(cz,'(i1)') iz-1 else write(cz,'(i2)') iz-1 endif jvfile = "./clear.jv_" // cz open(unit=31,file=jvfile) read(31,*) read(31,*,err=2000,end=2010) lr1, nphot read(31,*) read(31,'(12i6)',err=2000,end=2010) (irphot(i), i=1,lr1) read(31,*) read(31,'(7e11.3)',err=2000,end=2010) (tjval(j), j=1,nphot) do i=1,lr1 read(31,'(a)') astrg write(*,'(a,3x,a)') jvfile, astrg read(31,'(7e11.3)',err=2000,end=2010) _ (xjval(iz,j,i), j=1,nphot) if(astrg(1:4) .eq. "terp" .or. astrg(1:4) .eq. "MACR") then do j=1,nphot xjval(iz,j,i) = 0.5 * xjval(iz,j,i) enddo write(*,*) astrg(1:4), xjval(iz,48,i) endif enddo close(31) enddo endif write(*,*)'done photolysis frequencies ' c----------------------------------------------------------------------c write(*,*) write(*,*) ' cw and rain names ' do i=1,naq cxa(i) = cx(aqmap(i)) // 'a' do j=1,5 if(cxa(i)(j:j) .eq. ' ') then cxa(i)(j:j) = '_' endif end do write(*,'(a,3x,a)') cxa(i) end do C check reactions and verify assignments of types c itype1(i) - photolysis - already assigned at read c itype2(i) - simple arrhenius c itype3(i) - a + m or a + b + m, possibly also arrhenius c itype4(i) - troe - already assigned at read c itype5(i) - special functions for rate constants c itype6(i) - aqueous reactions need rate constants modified c itype7(i) - aqueous photolysis reactions nr1 = 0 nr2 = 0 nr3 = 0 nr4 = 0 nr5 = 0 nr6 = 0 nr7 = 0 do 30 j = 1, nreact C--------------------------------------------------------------------------- C Count photolysis reactions: if (molec(j,2) .eq. 'HV ' .and. c2(j)(1:1) .ne. 'a') then nr1 = nr1 + 1 itype1(nr1)=j endif C--------------------------------------------------------------------------- C Identify and count simple Arrhenius type: if (s(j) .ne. 0.) then nr2 = nr2 + 1 itype2(nr2) = j endif do 22 k = 1, 3 c--------------------------------------------------------------------------- c identify and count reactions with m as a reagent (not troe): if (molec(j,k) .eq. 'M ') then nr3 = nr3 + 1 itype3(nr3) = j endif 22 continue c--------------------------------------------------------------------------- c count troe reactions if ((molec(j,2).eq.'(M) ').or.(molec(j,3).eq.'(M) ')) then nr4 = nr4 + 1 endif c--------------------------------------------------------------------------- c identify and count reactions requiring special treatment: c these have a lower case 's' in column 2 if (c2(j)(2:2) .eq. 's') then nr5 = nr5 + 1 itype5(nr5) = j endif c--------------------------------------------------------------------------- c identify and count aqueous reactions c these have a lower case 'a' in column 1 if (c2(j)(1:1) .eq. 'a') then nr6 = nr6 + 1 itype6(nr6) = j c write(8,'(30x,a,2i4)') 'aqueous: ', j, nr6 endif naqr = nr6 C--------------------------------------------------------------------------- C Count aqueous photolysis reactions: if (molec(j,2) .eq. 'HV ' .and. c2(j)(1:1) .eq. 'a') then nr7 = nr7 + 1 itype7(nr7)=j endif C--------------------------------------------------------------------------- C Convert blank stoichiometry coefficients to unity (except if species is C also blank. Note: specific products cannot be "turned off" by C simply giving them zero stoichiometry coefficients. do 32 k = 1, 7 if (molec(j,k) .eq. 'HV ') molec(j,k) = iblank if (molec(j,k) .eq. 'M ') molec(j,k) = iblank c if (molec(j,k) .eq. 'N2 ') molec(j,k) = iblank c if (molec(j,k) .eq. 'O2 ') molec(j,k) = iblank if (molec(j,k) .eq. '(M) ') molec(j,k) = iblank if (sc(j,k) .eq. 0.) sc(j,k)=1. if (molec(j,k) .eq. iblank) sc(j,k) = 0. 32 continue 30 continue c map aqueous photolysis reactions to gas analogies do j=1,nr7 jr = itype7(j) do it=1,nr1 ir = itype1(it) if(molec(ir,1) .eq. molec(jr,1)) then hvmap(j) = ir endif end do write(*,*) 'hvmap ', j, hvmap(j) end do if (nr1 .ne. mr1) write(6,*) nr1,mr1, | 'warning: incorrect number of photolysis reactions' if (nr2 .ne. mr2) write(6,*) nr2, mr2, | 'warning: incorrect number of arrhenius reactions' if (nr3 .ne. mr3) write(6,*) nr3, mr3, | 'warning: incorrect number of m reactions' if (nr4 .ne. mr4) write(6,*) nr4, mr4, | 'warning: incorrect number of troe reactions' if (nr5 .ne. mr5) write(6,*) nr5, mr5, | 'warning: incorrect number of special reactions' if (nr6 .ne. mr6) write(6,*) nr6, mr6, | 'warning: incorrect number of aqueous reactions' if (nr7 .ne. mr7) write(6,*) nr7, mr7, | 'warning: incorrect number of aqueous photo reactions' print*,' ' print*,'REACTION MECHANISM CHECKED CORRECTLY ' print*,' ' C--------------------------------------------------------------------------- C find correct index for special reactions: C 1 - HO2 + HO2 -> H2O2 + O2 C 2 - HNO3 + HO -> H2O + NO3 C 3 - CO + HO -> CO2 + H do 50 i = 1, nr5 iwhich(i) = 0 if((molec(itype5(i),1) .eq. 'HO2 ') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 1 if((molec(itype5(i),1) .eq. 'HNO3') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 2 if((molec(itype5(i),1) .eq. 'CO ') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 3 if((molec(itype5(i),1) .eq. 'O3 ') .and. | (molec(itype5(i),2) .eq. ' ')) iwhich(i) = 4 if((molec(itype5(i),1) .eq. 'k031') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 5 if((molec(itype5(i),1) .eq. 'XOOH') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 6 if((molec(itype5(i),1) .eq. 'DMS') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 7 c account for special reactions that take place in the aqueous phase if(c2(itype5(i))(1:1) .eq. 'a') then if((molec(itype5(i),1) .eq. 'HO2 ') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 11 if((molec(itype5(i),1) .eq. 'HO ') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 12 if((molec(itype5(i),1) .eq. 'O3 ') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 13 if((molec(itype5(i),1) .eq. '2011') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 14 if((molec(itype5(i),1) .eq. 'a011') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 15 if((molec(itype5(i),1) .eq. 'SO2 ') .and. | (molec(itype5(i),2) .eq. 'H2O2')) iwhich(i) = 16 if((molec(itype5(i),1) .eq. 'SO2 ') .and. | (molec(itype5(i),2) .eq. 'O3 ')) iwhich(i) = 17 if((molec(itype5(i),1) .eq. 'CO2 ') .and. | (molec(itype5(i),2) .eq. 'HO ')) iwhich(i) = 18 if((molec(itype5(i),1) .eq. 'CO2 ') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 19 if((molec(itype5(i),1) .eq. 'CO3M') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 20 if((molec(itype5(i),1) .eq. 'NO3 ') .and. | (molec(itype5(i),2) .eq. 'HO2 ')) iwhich(i) = 21 if((molec(itype5(i),1) .eq. 'HO ') .and. | (molec(itype5(i),2) .eq. 'ad21')) iwhich(i) = 22 if((molec(itype5(i),1) .eq. 'HO ') .and. | (molec(itype5(i),2) .eq. 'aa21')) iwhich(i) = 23 if((molec(itype5(i),1) .eq. 'HO ') .and. | (molec(itype5(i),2) .eq. 'ak33')) iwhich(i) = 24 if((molec(itype5(i),1) .eq. 'HO ') .and. | (molec(itype5(i),2) .eq. 'a021')) iwhich(i) = 25 if((molec(itype5(i),1) .eq. 'HO ') .and. | (molec(itype5(i),2) .eq. 'ao23')) iwhich(i) = 26 if((molec(itype5(i),1) .eq. 'H2O2') .and. | (molec(itype5(i),2) .eq. 'aa21')) iwhich(i) = 27 if((molec(itype5(i),1) .eq. 'H2O2') .and. | (molec(itype5(i),2) .eq. 'ad21')) iwhich(i) = 28 if((molec(itype5(i),1) .eq. 'H2O2') .and. | (molec(itype5(i),2) .eq. 'ak33')) iwhich(i) = 29 c if((molec(itype5(i),1) .eq. 'HO2 ') .and. c | (molec(itype5(i),2) .eq. 'CL2M')) iwhich(i) = 14 endif !aq.phase special reactions if(iwhich(i) .eq. 0) then write(6,*) | 'warning: could not identify special reaction' endif write(*,*) 'specrxn ', itype5(i), iwhich(i), i 50 continue c determine the reaction matrix icmap call maprxns(molec, sc, iblank, itype6, nspec) c find rates for all reactions do 3000 ix=1,nspec ! print*,ix,' ' write(*,*) ' ' print*,' reactions for species: ',cx(ix) nsnks=icmap(ix,1) nsrcs=icmap(ix,2) print*, ' SINKS: ' do 3001 ir=1,nsnks ip=(ir-1)*5 ir1=icmap(ix,2+ip+1) is1=icmap(ix,2+ip+2) ip1=icmap(ix,2+ip+3) ip2=icmap(ix,2+ip+4) ip3=icmap(ix,2+ip+5) write(6,601) ir,ir1,stc(is1),cx(ip1),cx(ip2),cx(ip3),a(ir1),s(ir1) 3001 continue iy = gasmap(ix) if(iy .ge. 1 .and. iy .le. naq) then nsnksa=icmapa(iy,1) print*,' ' print*, ' AQUEOUS SINKS: ' do ir=1,nsnksa ip=(ir-1)*5 ir1=icmapa(iy,2+ip+1) is1=icmapa(iy,2+ip+2) ip1=icmapa(iy,2+ip+3) ip2=icmapa(iy,2+ip+4) ip3=icmapa(iy,2+ip+5) write(6,601) ir,ir1,stc(is1),cx(aqmap(ip1)),cx(aqmap(ip2)), _ cx(aqmap(ip3)),a(ir1),s(ir1) end do endif nsrcs=icmap(ix,2) print*,' ' print*,' SOURCES: ' do 3002 ir=1,nsrcs ip=(ir-1)*5 ir1=icmap(ix,2+5*nsnks+ip+1) is1=icmap(ix,2+5*nsnks+ip+2) ip1=icmap(ix,2+5*nsnks+ip+3) ip2=icmap(ix,2+5*nsnks+ip+4) ip3=icmap(ix,2+5*nsnks+ip+5) write(6,601) ir,ir1,stc(is1),cx(ip1),cx(ip2),cx(ip3), & a(ir1),s(ir1) 3002 continue iy = gasmap(ix) if(iy .ge. 1 .and. iy .le. naq) then nsrcsa=icmapa(iy,2) print*,' ' print*,' AQUEOUS SOURCES: ' do ir=1,nsrcsa ip=(ir-1)*5 ir1=icmapa(iy,2+5*nsnksa+ip+1) is1=icmapa(iy,2+5*nsnksa+ip+2) ip1=icmapa(iy,2+5*nsnksa+ip+3) ip2=icmapa(iy,2+5*nsnksa+ip+4) ip3=icmapa(iy,2+5*nsnksa+ip+5) write(6,601) ir,ir1,stc(is1), cx(aqmap(ip1)), cx(aqmap(ip2)), _ cx(aqmap(ip3)),a(ir1),s(ir1) end do endif 3000 continue return 2000 continue write(6,1140) finame stop 2010 continue write(6,1150) finame stop 601 format( i3,1x,i3,1x,f4.1,1x,a4,1x,a4,1x,a4,1x,e10.3,1x,e10.3) 1040 FORMAT (12I10) 1050 FORMAT (7E10.3) 1080 FORMAT (I5,1X,2(E10.3,1X,F6.2,1X),F5.2,2(1X,E10.3)) 1090 FORMAT (A2,3(A4,1X),2X, _ 3(F5.2,1X,A4,1X),F5.2,1X,A4,E8.1,1X,E8.1) 1091 FORMAT (I3,1x,A2,3(A4,1X),3(F5.2,1X,A4,1X),F5.2,1X,A4,2E10.3) 1140 FORMAT ('Error in Reading Input File:',A20) 1150 FORMAT ('Unexpected End-Of-File in Input File:',A20) end c----------------------------------------------------------------------- subroutine maprxns(molec, sc, iblank, iaq, nsp) c determine the reaction matrix icmap implicit none include 'chem_params.inc' c Input arguments integer nsp integer iaq(nrmx) ! index for aqueous reactions character*4 molec(nrmx,7) ! reaction molecules character*4 iblank ! blank character real sc(nrmx,7) ! stoichiometric coefficients c common blocks include 'cchem.com' include 'rchem.com' include 'ichem.com' include 'amap.com' c indices integer ix, ir, j, is, jp, jspec, ip, irct, i, na, n, iy integer nsrcs, nsnks integer nsrcsa, nsnksa integer nrct integer isrc(5*nrps+2) , ! temporary storage of map of sinks & isnk(5*nrps+2) ! temporary storage of map of sources integer isrca(5*nrps+2), ! temporary storage of map of sinks & isnka(5*nrps+2) ! temporary storage of map of sources real csnk c----------------------------------------------------------------------c c initialize do ir=1,5*nrps+2 do ix=1,nsp icmap(ix,ir) = 0 end do isnk(ir) = 0 isrc(ir) = 0 end do do ir=1,5*nrps+2 do ix=1,naq icmapa(ix,ir) = 0 end do isnka(ir) = 0 isrca(ir) = 0 end do c loop through all the species do 1000 ix=1,nsp call iinit(isnk, (5*nrps+2),0) call iinit(isrc, (5*nrps+2),0) call iinit(isnka,(5*nrps+2),0) call iinit(isrca,(5*nrps+2),0) c initialize number of sources and sinks for each species to 0 nsnks = 0 nsrcs = 0 nsnksa = 0 nsrcsa = 0 c loop through all the rxns and look for sinks and sources for species ix c first sinks, then sources do 900 ir=1,nreact csnk=0 ! stoich coefficient for sinks c if species appears as reactant sink is found for species ix for reaction ir do 110 j=1,3 if(molec(ir,j).eq.cx(ix)) then csnk=csnk+sc(ir,j) endif 110 continue if(csnk.ne.0) then ! sink reaction is found c set index of sink reaction for species ix do na=1,naqr if(ir .eq. iaq(na)) then nsnksa=nsnksa+1 ! total number of sink reactions isnka((nsnksa-1)*5+1)=ir c set index of stoich coefficient for species ix for reaction ir: c csnk is the stoich coeff. do is=1,2*nspec+1 if(abs((stc(is)-csnk)).lt..01) then isnka((nsnksa-1)*5+2)=is goto 113 elseif(stc(is).eq.0) then stc(is)=csnk isnka((nsnksa-1)*5+2)=is goto 113 elseif(is.gt.2*nspec) then print*,'not enough space for stoich. coefficients' print*,'increase dimensions of rstoich' stop endif end do 113 continue c fill in the index of the reactants nrct=0 do j=1,3 if(molec(ir,j).eq.iblank) goto 114 do irct=1,nspec if(molec(ir,j).eq.cx(irct)) then if((cx(irct) .eq. cx(ix)) .and. & (isnka(((nsnksa-1)*5)+2+3) .eq. 0)) then isnka(((nsnksa-1)*5)+2+3)=gasmap(irct) goto 114 else nrct=nrct+1 isnka(((nsnksa-1)*5)+2+nrct)=gasmap(irct) goto 114 endif endif end do 114 continue end do go to 145 !skip gas reactions endif end do c Sink of gas reaction nsnks=nsnks+1 ! total number of sink reactions isnk((nsnks-1)*5+1)=ir c set index of stoich coefficient for species ix for reaction ir: c csnk is the stoich coeff. do 120 is=1,2*nspec+1 if(abs((stc(is)-csnk)).lt..01) then isnk((nsnks-1)*5+2)=is goto 121 elseif(stc(is).eq.0) then stc(is)=csnk isnk((nsnks-1)*5+2)=is goto 121 elseif(is.gt.2*nspec) then print*,'not enough space for stoich. coefficients' print*,'increase dimensions of rstoich' stop endif 120 continue 121 continue c fill in the index of the reactants nrct=0 do j = 1,3 if(molec(ir,j) /= iblank) then do irct = 1,nspec if(molec(ir,j).eq.cx(irct)) then if((cx(irct) .eq. cx(ix)) .and. & (isnk(((nsnks-1)*5)+2+3) .eq. 0)) then isnk(((nsnks-1)*5)+2+3)=irct cycle else nrct=nrct+1 isnk(((nsnks-1)*5)+2+nrct)=irct cycle endif endif enddo endif enddo endif ! end of sink indexing for reaction ir 145 continue c c if species appears as product source is found for species ix for reaction ir do 150 jp=4,7 if(molec(ir,jp).eq.cx(ix)) then ! source reaction is found jspec=jp do na=1,naqr if(ir .eq. iaq(na)) then nsrcsa=nsrcsa+1 ! total number of sink reactions isrca((nsrcsa-1)*5+1)=ir c set index of stoich coefficient for species ix for reaction ir: c stoich(ir,jspec) do is=1,2*nspec+1 if(abs((stc(is)-sc(ir,jspec))).lt..01) then isrca((nsrcsa-1)*5+2)=is goto 157 elseif(stc(is).eq.0) then stc(is)=sc(ir,jspec) isrca((nsrcsa-1)*5+2)=is goto 157 elseif(is.gt.2*nspec) then print*,'not enough space for stoich. coefficients' print*,'increase dimensions of rstoich' stop endif end do 157 continue c fill in the index of the reactants nrct=0 do j=1,3 if(molec(ir,j).eq.iblank) goto 158 do irct=1,nspec if(molec(ir,j).eq.cx(irct)) then nrct=nrct+1 isrca(((nsrcsa-1)*5)+2+nrct)=gasmap(irct) goto 158 endif end do 158 continue end do go to 151 !source reaction found; get out of loop end if end do c set index of source reaction for species ix nsrcs=nsrcs+1 ! number of source reactions isrc(((nsrcs-1)*5)+1)=ir c set index of stoich coefficient for species ix for reaction ir: c stoich(ir,jspec) do 160 is=1,2*nspec+1 if(abs((stc(is)-sc(ir,jspec))).lt..01) then isrc((nsrcs-1)*5+2)=is goto 161 elseif(stc(is).eq.0) then stc(is)=sc(ir,jspec) isrc((nsrcs-1)*5+2)=is goto 161 elseif(is.gt.2*nspec) then print*,'not enough space for stoich. coefficients' print*,'increase dimensions of rstoich' stop endif 160 continue 161 continue c fill in the index of the reactants nrct=0 do 180 j=1,3 if(molec(ir,j).eq.iblank) goto 180 do 170 irct=1,nspec if(molec(ir,j).eq.cx(irct)) then nrct=nrct+1 isrc(((nsrcs-1)*5)+2+nrct)=irct goto 180 endif 170 continue 180 continue goto 151 ! source reaction found: skip the loop endif ! end of source indexing for reaction ir 150 continue 151 continue 900 continue c looped through all reactions: c total number of sources and sinks known for species ix: c now set icmap for species ix. c gas reactions: icmap(ix,1)=nsnks icmap(ix,2)=nsrcs ip=2 do 400 i=1,nsnks*5 ip=2+i icmap(ix,ip)=isnk(i) 400 continue do 500 i=1,nsrcs*5 icmap(ix,ip+i)=isrc(i) 500 continue is = jp c aqueous reactions: iy = gasmap(ix) if(iy .ge. 1 .and. iy .le. naq) then icmapa(iy,1)=nsnksa icmapa(iy,2)=nsrcsa ip=2 do i=1,nsnksa*5 ip=2+i icmapa(iy,ip)=isnka(i) end do do i=1,nsrcsa*5 icmapa(iy,ip+i)=isrca(i) end do endif ! 1 <= iy <= naq 1000 continue ! icmap set for all species return end c======================================================================c subroutine getph(t, rho, qc, aspec, ph) c this subroutine computes the pH based on trial and error implicit none include 'chem_params.inc' c c Input parameters real t ! temperature (K) real rho ! air density (kg/m3) real qc ! CW mixing ratio (g/kg) real aspec(naq) ! (molec/cm3) c Output parameter real ph c Local variables: real lwc ! liquid water content (cm3 H2O / cm3 air) real na ! Avogadro's number real k1s ! dissociation constant of H2SO3 <--> HSO3- real k2s ! dissociation constant of HSO3 <--> SO3= real k1c ! dissociation constant of H2CO3 <--> HCO3- real k2c ! dissociation constant of HCO3 <--> CO3= real k1o2m ! dissociation constant of HO2 <--> O2- real k1fa ! dissociation constant of HCOOH <--> HCOO- real k1aca ! dissociation constant of CH3COOH <--> CH3COO- + H+ real k1aom ! dissociation constant of HOCH2COOH <--> HOCH2COO- + H+ real k1agm ! dissociation constant of CHOCOOH <--> CHOCOO- + H+ real k1apm ! dissociation constant of CH3COCOOH <--> CH3COCOO- + H+ real k1axm ! dissociation constant of HOOCCOOH <--> HOOCCOO- + H+ real k2axm ! dissociation constant of HOOCCOO- <--> -OOCCOO- + H+ real fact ! temperature factor for dissoc. constants real hion ! [H+] in mol/L or M real ahion ! [H+] in mol/L or M (temporary) real f_o2m ! fraction of HO2(a) + O2- that is O2- real f_fo ! fraction of HCOOH(a) + HCOO- that is HCOO- real f_hso3 ! fraction of S(IV) that is HSO3- real f_hco3 ! fraction of H2CO3 + HCO3- + CO3= that is HCO3- real f_aca ! fraction of CH3COOH(a) + CH3COO- that is CH3COO- real f_aom ! fraction of HOCH2COOH(a) + HOCH2COO- that is HOCH2COO- real f_agm ! fraction of CHOCOOH(a) + CHOCOO- that is CHOCOO- real f_apm ! fraction of CH3COCOOH(a) + CH3COCOO- that is CH3COCOO- real f_axm ! fraction of HOOCCOOH(a) + HOOCCOO- + -OOCCOO- that is HOOCCOO- real so4a ! aqueous conc of SO4 [M] real nh4a ! aqueous conc of NH4 [M] real no3a ! aqueous conc of HNO3 [M] real so2a ! aqueous conc of S(IV) [M] real co2a ! aqueous conc of H2CO3 + HCO3- + CO3= [M] real ho2a ! aqueous conc of HO2 + O2- [M] real faa ! aqueous conc of HCOOH + HCOO- [M] real aca ! aqueous conc of CH3COOH + CH3COO- [M] real aoa ! aqueous conc of HOCH2COOH + HOCH2COO- [M] real aga ! aqueous conc of CHOCOOH + CHOCOO- [M] real apa ! aqueous conc of CH3COCOOH + CH3COCOO- [M] real axa ! aqueous conc of HOOCCOOH + HOOCCOO- + -OOCCOO- [M] !---------------------------------------------------------------------- ! Henry's Law Coefficients are written as Heff = KH + KH K1/[H+] + KH K1 K2/[H+]^2 ! where KH=kh_298*exp(dhr*(1/T-1/298)) ! where K1=k1_298*exp(k1hr*(1/T-1/298)) ! where K2=k2_298*exp(k2hr*(1/T-1/298)) include 'henryslawconsts.f' ! this has both declarations and code and must be placed at the transition !---------------------------------------------------------------------- c Electroneutrality equation: c Na+ + NH4+ + H+ = 2SO4= + NO3- + HSO3- + HCO3- + Cl- + O2- + HCOO- + CH3COO- + HOCH2COO- + CHOCOO- + CH3COCOO- + HOOCCOO- + -OOCCOO- c assume Na+ = Cl- (i.e. NaCl aerosol) c assume HSO3- dominates S(IV) (which it does for pH of most cloudwater) c to get: c H+ + NH4+ = 2*SO4= + NO3- + HSO3- + HCO3- + O2- + HCOO- c assume aspec(kso4) = SO4= entirely (all SO4 is in drops) c aspec(khno3a) = NO3- (all HNO3 in drops is NO3- and not HNO3(a)) c other species found by dissociation equilibrium: c aspec(kso2) = concentration of S(IV) family, so HSO3- = f_hso3 * aspec(kso2) c where f_hso3 = HSO3/S(IV) = K1*H+ / (H+*H+ + K1*H+ + K1*K2) c aspec(kfa) = concen of HCOOH(a) + HCOO-, so HCOO- = f_fo * aspec(kfa) c where f_fo = K1/(H+ + K1) ! 7/26/2017 add organic acids c H+ + NH4+ = 2*SO4= + NO3- + HSO3- + HCO3- + O2- + HCOO- + CH3COO- + HOCH2COO- + CHOCOO- + CH3COCOO- + HOOCCOO- c ignoring -OOCCOO- That is, we are assuming most of the oxalate in cloud water is the single ion (same assumption as for H2SO3) c organic anions are found in same way formate is found, see aspec(kfa) comment above c NH4+ is assigned. c because these fractions depend on H+, the calculation will be iterated c to determine H+ and pH if(qc .lt. 1.e-12) return !! na = 6.0221e23 c!! Convert concentrations to mol/L !! lwc = rho*qc * 1.e-6 !1.e-6 converts g H2O to m3 H2O !! nh4a = 7.16e-6 !From conditions of the simulation where NH4=7.16uM, SO4=21.066uM !!! nh4a = 40.46e-6 !From ALSC measurements where NH4=40.46uM, SO4=15.73uM (but still use SO4=21.066uM) !! so4a = aspec(kso4a) *1000./ (na * lwc) !! no3a = aspec(khno3a)*1000./ (na * lwc) !! so2a = aspec(kso2a) *1000./ (na * lwc) !! co2a = aspec(kco2a) *1000./ (na * lwc) !! ho2a = aspec(kho2a) *1000./ (na * lwc) !! faa = aspec(kfaa) *1000./ (na * lwc) !! aca = aspec(ka021a)*1000./ (na * lwc) !! aoa = aspec(kao23a)*1000./ (na * lwc) !! aga = aspec(kad21a)*1000./ (na * lwc) !! apa = aspec(kak33a)*1000./ (na * lwc) !! axa = aspec(kaa21a)*1000./ (na * lwc) c----------------------------------------------------------------------- !! write(*,*) 'SO4, HNO3, SO2, CO2, HO2, HCOOH ' !! write(*,'(1p,(8e12.4))') aspec(kso4a), aspec(khno3a), !! _ aspec(kso2a), aspec(kco2a), !! _ aspec(kho2a), aspec(kfaa) !! write(*,'(1p,(8e12.4))') so4a, no3a, so2a, co2a, ho2a, faa c!! use initial guess of pH for this case !! if(so4a.eq.0. .and. no3a.eq.0. .and. so2a.eq.0. .and. !! _ co2a.eq.0. .and. ho2a.eq.0. .and. faa .eq.0.) return !! fact = (1./t) - 1./298. !! k1fa = k1_298(kfaa) *exp(k1hr(kfaa) *fact) !1.77e-4 !! k1s = k1_298(kso2a)*exp(k1hr(kso2a)*fact) !1.74e-2*exp(-1940.*fact) !! k2s = k2_298(kso2a)*exp(k2hr(kso2a)*fact) !6.22e-8*exp(-1960.*fact) !! k1c = k1_298(kco2a)*exp(k1hr(kso2a)*fact) !7.7e-7*exp(1000.*fact) !! k2c = k2_298(kco2a)*exp(k2hr(kso2a)*fact) !4.84e-11*exp(1760.*fact) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! k1aca= k1_298(ka021a)*exp(k1hr(ka021a)*fact) ! CH3COOH --> CH3COO- + H+ !! k1aom= k1_298(kao23a)*exp(k1hr(kao23a)*fact) ! HOCH2COOH --> HOCH2COO- + H+ !! k1agm= k1_298(kad21a)*exp(k1hr(kad21a)*fact) ! CHOCOOH --> CHOCOO- + H+ !! k1apm= k1_298(kak33a)*exp(k1hr(kak33a)*fact) ! CH3COCOOH --> CH3COCOO- + H+ !! k1axm= k1_298(kaa21a)*exp(k1hr(kaa21a)*fact) ! HOOCCOOH --> HOOCCOO- + H+ !! k2axm= k2_298(kaa21a)*exp(k2hr(kaa21a)*fact) ! HOOCCOO- --> -OOCCOO- + H+ !! write(*,*) 'k1ho2, k1fa, k1s, k2s, k1c, k2c ' !! write(*,'(1p,(15e12.4))') k1o2m, k1fa, k1s, k2s, !! _ k1c, k2c, k1aca, k1aom, k1agm, k1apm, k1axm, k2axm !! write(*,*) 'initial pH = ',ph !! 5 continue !! hion = 10.**(-ph) !!! write(*,*) ' hion = ',hion !! f_o2m = k1o2m /(k1o2m + hion) !! f_fo = k1fa/(k1fa + hion) !! f_hso3 = k1s*hion/(hion*hion + k1s*hion + k1s*k2s) !! f_hco3 = k1c*hion/(hion*hion + k1c*hion + k1c*k2c) !! f_aca = k1aca/(k1aca + hion) !! f_aom = k1aom/(k1aom + hion) !! f_agm = k1agm/(k1agm + hion) !! f_apm = k1apm/(k1apm + hion) !! f_axm = k1axm*hion/(hion*hion + k1axm*hion + k1axm*k2axm) !!! write(*,*) 'f_o2m, f_fo, f_hso3, f_hco3, f_aca, f_aom, f_agm ' !!! write(*,'(1p,(15e12.4))') f_o2m, f_fo, f_hso3, f_hco3, !!! _ f_aca, f_aom, f_agm, f_apm, f_axm !! ahion = 2*so4a + no3a + f_hso3*so2a + f_hco3*co2a + f_o2m*ho2a + !! _ f_fo*faa + f_aca*aca + f_aom*aoa + f_agm*aga + !! _ f_apm*apa + f_axm*axa - nh4a !! ph = -log10(ahion) !! if(ph .lt. 0. .or. ph .gt. 10. .or. ahion .lt. 0.) then !! write(*,*) 'pH = ', ph, ' stopping program' !! write(*,*) 'ahion, hion, 2*so4a, no3a, f_hso3*so2a, f_hco3*co2a, !! _f_o2m*ho2a, f_fo*faa, f_aca*aca, f_aom*aoa, f_agm*aga, f_apm*apa, !! _ f_axm*axa, nh4a ' !! write(*,'(1p,(14e12.4))') ahion, hion, 2*so4a, no3a, !! _ f_hso3*so2a, f_hco3*co2a, f_o2m*ho2a, f_fo*faa, !! _ f_aca*aca, f_aom*aoa, f_agm*aga, f_apm*apa, f_axm*axa, !! _ nh4a !! stop !! endif !! if(abs(hion-ahion) .gt. 0.001*hion) go to 5 !! write(*,*) 'pH = ', ph !! write(*,*) 'ahion, hion, 2*so4a, no3a, f_hso3*so2a, f_hco3*co2a, !! _f_o2m*ho2a, f_fo*faa, f_aca*aca, f_aom*aoa, f_agm*aga, f_apm*apa, !! _ f_axm*axa, nh4a ' !! write(*,'(1p,(14e12.4))') ahion, hion, 2*so4a, no3a, !! _ f_hso3*so2a, f_hco3*co2a, f_o2m*ho2a, f_fo*faa, !! _ f_aca*aca, f_aom*aoa, f_agm*aga, f_apm*apa, f_axm*axa, !! _ nh4a end subroutine rates(t, rho, qv, qc, ph, rate, ratea, time, zalt) c this subroutine computes the rates for the various reactions implicit none include 'chem_params.inc' c c input parameters c t temp c rho air density (kg/m3) c qv water vapor mixing ratio (g/g) c qc cloud water mixing ratio c ph ph of liquid water c c output parameters c rate gas-phase reaction rates c ratea cloud reaction rates c c common block storage c c dx grid length x direction (m) c dy grid length y direction (m) c p0 reference pressure c dtchem time step c c cw molecular weights c cx character list of species c stc matrix of standard stoichometry coefficients c icmap reaction matrix for each species c nreact total number of reactions c a temperature independent coeff for arrehenius reactions c s temperature dependent coeff for arrehenius reactions c stc stoichmetry coefficients c c*** chemical reaction types c nr1,nr2,nr3,nr4,nr5 number of each reaction type c itype1,...itype5 reaction identifier for each reaction type c itype1 - photolysis c itype2 - simple arrhenius c itype3 - a + m or a + b + m, possibly also arrhenius c itype4 - troe - already assigned at read c itype5 - special functions for rate constants c c*** troe reaction data, molecule, cm3, sec units c ak0300 - zero pressure 300k rate constant, third order c an - temperature exponent for zero pressure rate constant c aki300 - high pressure 300k rate constant, second order c am - temperature exponent for high pressure rate constant c bas - base of exponentiation, 0.6 for most c aequil - pre-exponential of equilibrium constant c tequil - activation temperature of equilibrium constant c c temporary storage c c rate reaction rate c iwhich index indentifier for special reactions c press press c conv conversion factor c m concentration of air c h2o concentration of water c **variables used in calculating rates c fact c trel c ak0t c akit c f1 c f2 c f3 c f4 c forw c fh2o c rk0 c rk2 c rk3 c input variables real qv, qc, t, rho real ph real time real zalt ! altitude (m) of the parcel, needed for jvalues c output variable real rate(nrmx) real ratea(nrmx) c common /dltprms/dtchem, dt_dbg real dtchem ! the chemistry timestep real dt_dbg ! how often to do a debug print c chemical storage block include 'cchem.com' include 'ichem.com' include 'rchem.com' include 'irate.com' include 'rrate.com' include 'hvchem_jvalchgs.com' c temporary storage integer mr1, jtim, jtimp, irphot(30), ip real tday, tj, pslope, tdif, dtj, otim, delttj, utctime real zdif, zval, dzj, rk1 ! to include altitude interpolation integer jz, jzp1 real m, h2o, press real :: fact,trel,ak0t,akit,f1,f2,f3,f4,forw,fh2o, & rk0,rk2,rk3,conv,rcld real dmsnum, dmsden integer icld integer i,j,k,it,ir, kern real rstr(nrmx) real lww !liquid water content of cloud water real hion !H+ concentration in cloud/rain drops real na !Avogadro's number real f_ho2 !fraction of HO2 + O2- that is HO2 real f_o2m !fraction of HO2 + O2- that is O2- real f_fa !fraction of HCOOH + HCOO- that is HCOOH real f_fo !fraction of HCOOH + HCOO- that is HCOO- real f_hco3 !fraction of CO2 + HCO3- that is HCO3- real k1fa !dissociation const for formic acid real k1c !dissociation const for CO2 real k1s !dissociation const for H2SO3-->HSO3- real k2s !dissociation const for HSO3-->SO3= real k1o2m !dissociation const for CO2 real k1aca !dissociation const for CO2 real k1aom !dissociation const for CO2 real k1agm !dissociation const for CO2 real k1apm !dissociation const for CO2 real k1axm !dissociation const for CO2 real k2axm !dissociation const for CO2 real f_hso3 !fraction of SO2 + HSO3- + SO3= that is HSO3- real f_so3 !fraction of SO2 + HSO3- + SO3= that is SO3= real rk !temporary rate constant for a01A+HOA real rkh2o !rate constant for O1D + H2O real rkm !rate constant for O1D + M c+b+ylbj_newaq real f_so2 !fraction of SO2 + HSO3- + SO3= that is SO2 real f_aa !fraction of CH3COOH + CH3COO- that is CH3COOH (acetic acid) real f_ao !fraction of CH3COOH + CH3COO- that is CH3COO- real f_ga !fraction that is CHOCOOH (glyoxylic acid) real f_go !fraction that is CHOCOO- real f_oa !fraction that is COOHCOOH (oxalic acid) real f_oo !fraction that is COOHCOO- real f_ooo !fraction that is COO-COO- real f_pa !fraction that is CH3COCOOH (pyruvic acid) real f_po !fraction that is CH3COCOO- real f_gca !fraction that is HOCH2COOH (glycolic acid) real f_gco !fraction that is HOCH2COO- c+e+ylbj_newaq !---------------------------------------------------------------------- ! Henry's Law Coefficients are written as Heff = KH + KH K1/[H+] + KH K1 K2/[H+]^2 ! where KH=kh_298*exp(dhr*(1/T-1/298)) ! where K1=k1_298*exp(k1hr*(1/T-1/298)) ! where K2=k2_298*exp(k2hr*(1/T-1/298)) include 'henryslawconsts.f' ! this has both declarations and code and must be placed at the transition !---------------------------------------------------------------------- rate(1:nreact) = a(1:nreact) ratea(1:nreact) = 0. write(*,*) 'beg rates ', rate(4) c find concentrations of air and water c conv converts to from kg/m3 to grams/cm3 press = rho*287.*t ! N/m2 conv = 1.e-3*rho ! g/cm3 m = conv* (6.0221e23/(28.97 )) ! molec/cm3 h2o = conv*qv*(6.0221e23/(18.0152)) ! g air/cm3 * g H2O/g air * mol H2O/g H2O * molecules/mol --> molecules/cm3 write(*,*) 'M, H2O ', m, h2o, press, rho, qv, t, time C C PHOTOLYSIS REACTIONS C if (nr1 .gt. 0) then tday = 24.*3600. !!! Time is local time, but data in JVALUES/clear.jv files are in UTC time !! e.g. max jvalue at 66600 s = 18:30 UTC --> 1:30 pm local solar noon (verified on NOAA solar calculator) !! --> add 5 hours to tj !! utctime = time + 5.*3600. !mcb 08/23/2016 Changed TUV v5.3 calculation to include time zone and have max jvalue at solar noon utctime = time tj = utctime - tday*real(int(utctime/tday)) delttj = tjval(2)-tjval(1) jtim=int(tj/delttj) +1 jtimp=jtim+1 if(jtim.eq.nphot) then jtimp=1 elseif(jtim.gt.nphot) then print*,'something is wrong ' stop endif tdif = tjval(jtimp) - tjval(jtim) dtj = tj - tjval(jtim) print*,tjval(jtim),tj,tjval(jtimp),jtim print*,xjval(1,jtim,1),xjval(1,jtimp,1),tdif, dtj zdif = 1000. ! assume 1000 m altitude difference in input files -- it's how it is set up jz = int(zalt/zdif) + 1 jzp1 = jz + 1 zval = (jz-1)*1000. dzj = zalt - zval print*, zval, zalt, jz, xjval(jz,jtim,1), xjval(jzp1,jtim,1), dzj if(mod(time,300.).eq.0.) print*,'photolysis rates' do 30 it=1,nr1 ir=itype1(it) pslope=(xjval(jz,jtimp,it) - xjval(jz,jtim,it)) / tdif rk1 = abs(xjval(jz,jtim,it) + pslope*dtj) !jvalue at lower altitude pslope=(xjval(jzp1,jtimp,it) - xjval(jzp1,jtim,it)) / tdif rk2 = abs(xjval(jzp1,jtim,it) + pslope*dtj) !jvalue at higher altitude rate(ir) = rk1 + (rk2 - rk1) * dzj/zdif !jvalue at desired altitude if(mod(time,300.).eq.0.) print*,it,ir,rate(ir), rk1,rk2 30 continue 341 format(f10.1,1x,i2,1x,i2,e12.5) endif ! end of photolysis block c aqueous photolysis rates: if(nr7 .gt. 0) then do it=1,nr7 ir=itype7(it) ip = hvmap(it) ratea(ir) = 1.5 * rate(ip) end do endif C C THIS CALCULATES RATES CONSTANTS FOR SIMPLE ARRHENIUS C EXPRESSIONS FOR SELECTED REACTIONS if (nr2 .gt. 0) then do it = 1, nr2 ir=itype2(it) rate(ir) =a(ir)*exp(s(ir)/t) rstr(ir) = rate(ir) enddo endif do it=4,6 write(*,*) 'after arr ', it, rate(it) end do C THIS CALCULATES RATES CONSTANTS FOR SIMPLE ARRHENIUS C EXPRESSIONS WHICH INCLUDE THE COLLISION PARTNER M C NOTE THAT TROE AND REVERSE TROE ARE EXCLUDED. IF (NR3 .GT. 0) then do it = 1, nr3 ir=itype3(it) fact = 3.3556e-03 - 1./t rate(ir) = m*a(ir)*exp(fact*s(ir)) end do endif ! write(*,*) 'after M ', rate(4) do it=4,6 write(*,*) 'after M ', it, rate(it) end do IF (NR4 .GT. 0) then C THIS CALCULATES RATE CONSTANTS FOR REACTIONS OF THE TYPE C A + B + (M) -> C + (M) and C C + (M) -> A + B + (M) C USING TROE EXPRESSIONS AND EQUILIBRIUM CONSTANTS FOR THE REVERSE C TREL - temperature / 300. C AK0T - low press lim k at current temp C AKIT - high press. lim k at current temp C F1,F2,F3,F4 - temporary values C FORW - forward rate const. C______________________________________________________________________ do it = 1, nr4 ir=itype4(it) fact = (298. - t)/(298.*t) trel = t/300. ak0t = ak0300(it) * trel**(-an(it)) akit = aki300(it) * trel**(-am(it)) f1 = ak0t*m/akit f2 = ak0t*m/(1. + f1) f3 = log(max(f1,1.e-30))/2.303 f4 = 1./(1. + f3*f3) forw = f2 * bas(it)**f4 c print*,'troe reaction for reaction ',ir c print*,forw,aequil(it),tequil(it),fact rate(ir) = forw / ( aequil(it) * | exp(tequil(it)*(fact))) end do endif ! write(*,*) 'after (M) ', rate(4) do it=4,6 write(*,*) 'after (M) ', it, rate(it) end do if (nr5 .gt. 0) then c SPECIAL REACTIONS do 70 it = 1, nr5 ir=itype5(it) write(*,*) 'start special ', ir, rate(ir) C HO2+HO2 if (iwhich(it) .eq. 1) then fact = (1./t)- 3.3556e-03 rk2 = 1.7e-12*exp( 600.*(fact)) rk3 = 4.9e-32*exp(1000.*(fact)) fh2o = 1.4e-21*exp(2200./t) rate(ir) = (rk2 + rk3*m)*(1. + fh2o*h2o) endif C HNO3+HO if (iwhich(it) .eq. 2) then rk0 = 7.2e-15*exp(785./t) rk2 = 4.1e-16*exp(1440./t) rk3 = 1.9e-33*exp(725./t) rate(ir) = rk0 + rk3*m/(1. + rk3*m/rk2) endif C CO+HO if (iwhich(it) .eq. 3) then c print*,'pressure: ',press,press/101325. rate(ir)=1.5e-13*(1. + 0.6*press/(101325.)) endif C O3+hv if (iwhich(it) .eq. 4) then write(*,*) 'O3+hv ', rate(ir), h2o, m rkh2o = 2.2e-10 rkm = 2.9e-11 rate(ir)= rate(ir) * rkh2o*h2o / (rkh2o*h2o + rkm*m) endif C CH3COCH3+HO or k031+HO if (iwhich(it) .eq. 5) then rate(ir) = 8.8e-12*exp(-1320/t) + 1.7e-14*exp(423/t) endif C XOOH+HO if (iwhich(it) .eq. 6) then rate(ir) = (t*t)*7.69e-17*exp(253/t) endif C DMS+HO if (iwhich(it) .eq. 7) then dmsnum = t*exp(-234./t) + 8.46e-10*exp(7230/t) + _ 2.68e-10*exp(7810./t) dmsden = 1.04e11*t + 88.1*exp(7460./t) rate(ir) = dmsnum/dmsden endif c Special Reactions in Aqueous Reaction Mechanism c First determine [H+] concentration hion = 10.**(-ph) C HO2aq + O2- !! if (iwhich(it) .eq. 11) then !! fact = (1./t)- (1./298.) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! f_ho2 = hion/(k1o2m + hion) !! rate(ir)=9.7e7*exp(-1060.*fact)*f_o2m*f_ho2 + ! This is HO2 + O2M rxn and HO2 + HO2 rxn !! + 8.3e5 * exp(-2720.*fact)*f_ho2*f_ho2 !! endif C!! HOaq + O2- !! if (iwhich(it) .eq. 12) then !! fact = (1./t)- (1./298.) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! f_ho2 = hion/(k1o2m + hion) !! rate(ir)=(f_o2m*1.e10) * exp(-1500.*fact) !! endif C!! O3aq + O2- !! if (iwhich(it) .eq. 13) then !! fact = (1./t)- (1./298.) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! rate(ir)=rstr(ir) * f_o2m !! endif C!! CH3OOaq + O2- !! if (iwhich(it) .eq. 14) then !! fact = (1./t)- (1./298.) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! rate(ir)=rstr(ir) * f_o2m !! endif C!! HCOOHaq + OHaq and HCOO- + OHaq !! if (iwhich(it) .eq. 15) then !! fact = (1./t)- (1./298.) !! k1fa = k1_298(kfaa) *exp(k1hr(kfaa) *fact) !1.77e-4 !! f_fa = k1fa/(k1fa + hion) !! f_fo = hion/(k1fa + hion) c!! rk = (1.6e8*f_fo + 2.5e9*f_fa) * exp(-1510.*fact) ! coefficients from Lelieveld and Crutzen (1991) !! rk = 1.3e8*f_fo*exp(-1000.*fact) + !! _ 3.2e9*f_fa*exp(-1000.*fact) ! coefficients from Ervens et al (2008) !! rate(ir)= rk !! write(*,*) 'fa-oh ', !! _ 1.3e8*f_fo, 3.2e9*f_fa, rate(ir), ir !! endif c SO2(aq) + H2O2(aq) if(iwhich(it) .eq. 16) then !SO2 aqchem fact = (1./t)- (1./298.) k1s = k1_298(kso2a)*exp(k1hr(kso2a)*fact) !1.74e-2*exp(-1940.*fact) k2s = k2_298(kso2a)*exp(k2hr(kso2a)*fact) !6.22e-8*exp(-1960.*fact) f_hso3 = k1s*hion/(hion*hion + k1s*hion + k1s*k2s) !fraction of S(IV) that's HSO3 rk = 7.45e7*exp(-4430.*fact) ! Hoffmann and Calvert (1985) rate(ir) = f_hso3 * rk*hion/(1. + 13.*hion) write(*,*) 'so2per', rk, hion, hion/(1. + 13.*hion), + f_hso3, rate(ir) endif c SO2(aq) + O3(aq) if(iwhich(it) .eq. 17) then !SO2 aqchem fact = (1./t)- (1./298.) k1s = k1_298(kso2a)*exp(k1hr(kso2a)*fact) !1.74e-2*exp(-1940.*fact) k2s = k2_298(kso2a)*exp(k2hr(kso2a)*fact) !6.22e-8*exp(-1960.*fact) f_hso3 = k1s*hion/(hion*hion + k1s*hion + k1s*k2s) !fraction of S(IV) that's HSO3 f_so3 = k1s*k2s /(hion*hion + k1s*hion + k1s*k2s) !fraction of S(IV) that's HSO3 rk = 3.7e5*exp(-5530.*fact)*f_hso3 + ! Hoffmann and Calvert (1985) + 1.5e9*exp(-5280.*fact)*f_so3 rate(ir) = rk write(*,*) 'so2-o3', 3.7e5*exp(-5530.*fact)*f_hso3, + 1.5e9*exp(-5280.*fact)*f_so3, f_hso3, f_so3, + rate(ir) endif C!! HCO3- + OHaq !! if (iwhich(it) .eq. 18) then !! fact = (1./t)- (1./298.) !! k1c = k1_298(kco2a)*exp(k1hr(kso2a)*fact) !7.7e-7*exp(1000.*fact) !! f_hco3 = k1c/(k1c + hion) !! rk = 1.0e7*f_hco3 * exp(-1510.*fact) ! for test rk = 0. !! rate(ir)= rk !! endif C!! HCO3- + O2- !! if (iwhich(it) .eq. 19) then !! fact = (1./t)- (1./298.) !! k1c = k1_298(kco2a)*exp(k1hr(kso2a)*fact) !7.7e-7*exp(1000.*fact) !! f_hco3 = k1c/(k1c + hion) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! rk = 1.5e6*f_hco3 * f_o2m * exp(-1510.*fact) ! for test rk = 0. !! rate(ir)= rk !! endif C!! CO3- + O2- !! if (iwhich(it) .eq. 20) then !! fact = (1./t)- (1./298.) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! rk = 4.0e8 * f_o2m * exp(-1510.*fact) ! for test rk = 0. !! rate(ir)= rk !! endif C!! NO3 + O2- !! if (iwhich(it) .eq. 21) then !! fact = (1./t)- (1./298.) !! k1o2m= k1_298(kho2a)*exp(k1hr(kho2a)*fact) ! HO2 --> O2- + H+ !! f_o2m = k1o2m/(k1o2m + hion) !! rate(ir)=rstr(ir) * f_o2m !! endif c+b+ylbj_newaq C!! HO + ad21 +ylbj_newaq ! ad21 = CHOCOOH = glyoxylic acid !! if (iwhich(it) .eq. 22) then !! fact = (1./t)- (1./298.) !! k1agm= k1_298(kad21a)*exp(k1hr(kad21a)*fact) ! CHOCOOH --> CHOCOO- + H+ !! f_go = k1agm/(k1agm + hion) !! f_ga = hion/(k1agm + hion) !! rk = 3.6e8*f_ga * exp(-1000.*fact) + !! _ 2.9e9*f_go * exp(-4300.*fact) !! rate(ir)= rk !! endif C!! HO + aa21 +ylbj_newaq ! aa21 = HOOCCOOH = oxalic acid !! if (iwhich(it) .eq. 23) then !! fact = (1./t)- (1./298.) !! k1axm= k1_298(kaa21a)*exp(k1hr(kaa21a)*fact) ! HOOCCOOH --> HOOCCOO- + H+ !! k2axm= k2_298(kaa21a)*exp(k2hr(kaa21a)*fact) ! HOOCCOO- --> -OOCCOO- + H+ !! f_oa = hion*hion/(hion*hion + k1axm*hion + k1axm*k2axm) !! f_oo = k1axm*hion/(hion*hion + k1axm*hion + k1axm*k2axm) !! f_ooo= k1axm*k2axm/(hion*hion + k1axm*hion + k1axm*k2axm) !! rk = 1.4e6*f_oa + !! _ 1.9e8*exp(-2800.*fact)*f_oo + !! _ 1.6e8*exp(-4300.*fact)*f_ooo !! rate(ir) = rk !! endif C!! HO + ak33 +ylbj_newaq ! ak33 = CH3COCOOH = pyruvic acid !! if (iwhich(it) .eq. 24) then !! fact = (1./t)- (1./298.) !! k1apm= k1_298(kak33a)*exp(k1hr(kak33a)*fact) ! CH3COCOOH --> CH3COCOO- + H+ !! f_po = k1apm/(k1apm + hion) !! f_pa = hion/(k1apm + hion) !! rk = 3.2e8*exp(-1800.*fact)*f_pa + ! Herrmann et al (2015) !! _ 7.1e8*exp(-3000.*fact)*f_po !! rate(ir)= rk !! endif C!! HO + a021 +ylbj_newaq ! a021 = CH3OOH = acetic acid !! if (iwhich(it) .eq. 25) then !! fact = (1./t)- (1./298.) !! k1aca= k1_298(ka021a)*exp(k1hr(ka021a)*fact) ! CH3COOH --> CH3COO- + H+ !! f_ao = k1aca/(k1aca + hion) !! f_aa = hion/(k1aca + hion) !! rk = 1.5e7*f_aa * exp(-1330.*fact) + !! _ 1.0e8*f_ao * exp(-1800.*fact) !! rate(ir)= rk !! endif C!! HO + ao23 ! ao23 = HOCH2COOH = glycolic acid !! if (iwhich(it) .eq. 26) then !! fact = (1./t)- (1./298.) !! k1aom= k1_298(kao23a)*exp(k1hr(kao23a)*fact) ! HOCH2COOH --> HOCH2COO- + H+ !! f_gca = hion/(k1aom + hion) !! f_gco = k1aom/(k1aom + hion) !! rk = 1.2e9*f_gco + 5.4e8*f_gca !! rate(ir)= rk !! endif C!! H2O2 + aa21 +ylbj_newaq ! aa21 = HOOCCOOH = oxalic acid !! if(iwhich(it) .eq. 27) then !! fact = (1./t)- (1./298.) !! k1axm= k1_298(kaa21a)*exp(k1hr(kaa21a)*fact) ! HOOCCOOH --> HOOCCOO- + H+ !! k2axm= k2_298(kaa21a)*exp(k2hr(kaa21a)*fact) ! HOOCCOO- --> -OOCCOO- + H+ !! f_oa = hion*hion/(hion*hion + k1axm*hion + k1axm*k2axm) !! f_oo = k1axm*hion/(hion*hion + k1axm*hion + k1axm*k2axm) !! rk = 1.1e-1*f_oa + 1.5e-4*f_oo !! rk = 0. !! endif C!! H2O2 + ad21 +ylbj_newaq ! ad21 = CHOCOOH = glyoxylic acid !! if (iwhich(it) .eq. 28) then !! fact = (1./t)- (1./298.) !! k1agm= k1_298(kad21a)*exp(k1hr(kad21a)*fact) ! CHOCOOH --> CHOCOO- + H+ !! f_go = k1agm/(k1agm + hion) !! f_ga = hion/(k1agm + hion) !! rk = 0.9*f_ga + 16.5*f_go !! rk = 0. !! rate(ir)= rk !! endif C!! H2O2 + ak33 +ylbj_newaq ! ak33 = CH3COCOOH = pyruvic acid !! if (iwhich(it) .eq. 29) then !! fact = (1./t)- (1./298.) !! k1apm= k1_298(kak33a)*exp(k1hr(kak33a)*fact) ! CH3COCOOH --> CH3COCOO- + H+ !! f_po = k1apm/(k1apm + hion) !! f_pa = hion/(k1apm + hion) !! rk = 0.12*f_pa + 0.75*f_po ! Herrmann et al (2015) !! rk = 0. !! rate(ir)= rk !! endif c+e+ylbj_newaq C CL2- + HO2 and CL2- + O2- c if (iwhich(it) .eq. 14) then c fact = (1./t)- 3.3556e-03 c f_o2m = 3.5e-5/(3.5e-5 + hion) c f_ho2 = hion/(3.5e-5 + hion) cc f_o2m = 1.6e-5/(1.6e-5 + hion) ! K1 from Ervens et al. 2008 and her earlier refs cc f_ho2 = hion/(1.6e-5 + hion) c rk = (4.5e9*f_ho2 + 1.0e9*f_o2m) * exp(-1510.*fact) c rate(ir)= rk c endif C CL2- reactions c if (iwhich(it) .eq. 10) then c fact = (1./t)- 3.3556e-03 c k1cl = 5.3e-6 c f_cl2m= k1fa/(k1fa + hion) c f_cla = hion/(k1fa + hion) c f_clm = hion/(k1fa + hion) c rk = (1.6e8*f_fo + 2.5e9*f_fa) * exp(-1510.*fact) c rate(ir)= rk c endif write(*,*) 'end special ', ir, rate(ir) 70 continue endif ! write(*,*) 'after special ', rate(4) do it=4,6 write(*,*) 'after special ', it, rate(it) end do na = 6.0221e20 if(nr6 > 0) then c modify aqueous reaction rate constants (molar units to molec/cm3 units) lww = rho * qc/1.e6 rcld = real(nint(lww / max(lww,1.e-20))) if( nr7 > 0 ) then do it=1,nr7 ir = itype6(it) ratea(ir) = ratea(ir) *rcld rate(ir) = 0. end do endif write(*,*) 'in rates lww: ', lww, rcld do it=nr7+1,nr6 ir = itype6(it) if(ir < 7) then write(*,'(i5,1p3e12.3)') ir, rate(ir), na, lww endif ratea(ir) = rate(ir)/(na*max(lww, 1.e-20)) * rcld rate(ir) = 0. end do endif ! write(*,*) 'after aq ', rate(4) write(*,*) do it=4,6 write(*,*) 'after aqueous ', it, ratea(it) end do if(mod(time,600.) .eq. 0.) then print*,'reaction rates at time ', time do i=1,nreact write(*,'(i5,1pe11.2,e11.2,0p,f8.1)') i,rate(i), ratea(i),t enddo endif end subroutine rates c----------------------------------------------------------------------c c heff= Kh*Ra*T and is dimensionless c c----------------------------------------------------------------------c subroutine glrates(tair, rho, ph, ktc, lww, heff, phrat) implicit none include 'chem_params.inc' c Input variables real tair, rho, ph, lww c Output variables real ktc(naq) ! first-order rate for transfer of species c into/out of drop real heff(naq) ! dimensionless Henry's law coefficient real phrat(naq) ! phase ratio = Ca/Cg c Local variables integer ik, it real na, ahion data na /6.0221e23/ c Externals external gasflux, !gets ktc _ solub !gets heff, phrat ahion = 10.**(-ph) call gasflux(lww, tair, ktc) call solub(tair, ahion, lww, rho, heff, phrat) end subroutine glrates subroutine gasflux(lww, tair, ktc) c Parameterization of gas to drop transfer; follows Schwartz (1986) c implicit none include 'chem_params.inc' c Input variables real lww, tair c Output variables real ktc(naq) c Local variables integer ik, n real third parameter(third=1./3.) real taudg, taui real nd, amean, na, acc(naq) real ru, dg, thvel parameter(ru=8.314e7, dg=0.10, na=6.0221e23) real nor, rhol, pi !Nor=0.08 (cm-4), rhol=1g/cm3, pi data nor, rhol/0.08, 1./ real mw(naq) c----------------------------------------------------------------------- c c Initialize c acc = 0.05 is default value from Lelieveld & Crutzen 1991 c acc for VOCs are from McNeill et al (2012) supplementary info and c her references: (some refs may be for KH) c 5 = Betterton and Hoffmann (1988) c 6 = Lim et al (2005) c 7 = Kroll et al (2005) c 8 = Herrmann et al (2005) CAPRAM 3.0 c 11 = Khan et al (1995) c 12 = Davidovits et al (1995) c 17 = Sander and Crutzen (1996) c 16 = Sander (2012) NIST chemistry webbook c 18 = McNeill et al Atmos Chem (2012) c 19 = Lelieveld & Crutzen 1991 c 20 = Kolb et al (2010) pi = 4.*atan(1.) acc(ko3a) = 5.3e-4 acc(kh2o2a) = 0.02 acc(koha) = 0.05 acc(kho2a) = 0.2 acc(kso2a) = 0.05 acc(kso4a) = 0.2 ! assumed same as HNO3 mw(ko3a) = 48. mw(kh2o2a) = 34. mw(koha) = 17. mw(kho2a) = 33. mw(kso2a) = 64. mw(kso4a) = 96. ktc(1:naq) = 0. if(lww >= 1.e-15) then ! cloud water amean = 0.001 ! cm taudg = amean**2/(3.*dg) do n = 1,naq thvel = sqrt(8.*ru*tair/(pi*mw(n))) !cm/s taui = 4.*amean/(3.*thvel*acc(n)) if(taudg > 0. .or. taui > 0.) then ktc(n) = 1./(taudg + taui) endif end do endif end subroutine gasflux subroutine solub(temp, ahp, lww, rho, heff, phrat) ! ! calculate the effective henry's law constant and the phase ratio ! for each species ! implicit none include 'chem_params.inc' include 'cchem.com' include 'amap.com' c Input variables real temp, lww, rho real ahp c Output variables real heff(naq), phrat(naq) c Local variables integer isp, n real khen, k1, k2, ra c----------------------------------------------------------------------- c functions real ek, h298, hhr, tt, t ek(h298,hhr,tt) = h298*exp( hhr*((298. - tt)/(tt*298.)) ) !---------------------------------------------------------------------- ! Henry's Law Coefficients are written as Heff = KH + KH K1/[H+] + KH K1 K2/[H+]^2 ! where KH=kh_298*exp(dhr*(1/T-1/298)) ! where K1=k1_298*exp(k1hr*(1/T-1/298)) ! where K2=k2_298*exp(k2hr*(1/T-1/298)) !---------------------------------------------------------------------- include 'henryslawconsts.f' ! this has both declarations and code and must be placed at the transition !---------------------------------------------------------------------- ! calculate effective henry's law constants ! organic compounds KH from Sander (2015) ACP, but chosen value has specific reference !---------------------------------------------------------------------- do isp=1,naq khen = ek(kh_298(isp), dhr(isp), temp) k1 = ek(k1_298(isp), k1hr(isp), temp) k2 = ek(k2_298(isp), k2hr(isp), temp) heff(isp) = khen*(1. + k1/ahp *(1. + k2/ahp)) end do write(*,*) 'in solub', temp, lww ra = 8314./101325. do n=1,naq heff(n) = heff(n)*ra*temp !dimensionless form phrat(n) = heff(n)*lww !Phase ratio write(*,'(i4,1p,2e12.5,3x,a4)') n, heff(n), phrat(n), _ cx(aqmap(n)) end do end subroutine solub subroutine chem(rate,t,q,spec,ratea,ktc,heff,phrat,aspec,time) !---------------------------------------------------------------------- c this subroutine computes the rates for the various reactions c then solves the rate equations. they are solved with the Euler- c Backward iterative method c c Aug. 14/95 c For steady-state species, the fluxes and losses in the boundary c Layer are included in the chemistry scheme c c c Mar '96: adjusted for concentration mixing ratios c took away deposition and flux terms from lower layer c allowed chemistry in the first layer c c variables: see chemset for chemical storage block and parameters c c temporary storage c nsrcs number of sources for each species c nsnks number of sinks for each species c srca(b) sum of total sources for each species a=time n, b=time n+1 c xlsa(b) sum of total sinks for each species a=time n, b=time n+1 c xa(b) species concentration a=time n, b=time n+1 c xspec species concentration !---------------------------------------------------------------------- implicit none include 'chem_params.inc' !c++mcb !---------------------------------------------------------------------- ! dummy arguments !---------------------------------------------------------------------- real, intent(in) :: t, q real, intent(in) :: time real, intent(in) :: rate(nrmx) ! gas phase rxnt constants real, intent(in) :: ratea(nrmx) ! aqueous phase rxnt constants real, intent(in) :: ktc(naq) ! diffusion transfer rate coefs. real, intent(in) :: heff(naq) ! effective henry's law coef real, intent(in) :: phrat(naq) ! phase ratio = Ca/Cg real, intent(inout) :: spec(nspec) ! gas phase chemical species concentrations real, intent(inout) :: aspec(naq) ! aqueous phase chemical species concentrations c commons c chemical storage block real qca, lww, ph integer ncpts, npts common/cldpar/qca, lww, ph, ncpts, npts c grid storage block common /dltprms/dtchem, dt_dbg real dtchem ! the chemistry timestep real dt_dbg ! how often to do a debug print c chemical storage block include 'cchem.com' include 'ichem.com' include 'rchem.com' include 'amap.com' common /constant/pi,r,gr,cp real pi ! pi = PI real r ! gas constant for dry air = RD real gr ! acceleration of gravity = GZERO real cp ! specific heat of dry air at constant pressure = CPD c temporary storage real :: gasSrc, gasSnk real :: press, convu real :: xa(nspec), xb(0:nspec) integer :: ik, ir1, iy, iu integer :: ixa, nsnka, nsrca real xac(naq) real xbc(0:naq) real :: aqSrc, aqSnk real :: xold, xbn, xbcn real :: xtmp real :: ctot1t(naq), ctot1(naq) integer :: icyc, ncyc real :: xls, src, p1,p2, mair, mh2o real :: ccrit, xnew logical :: more,pconv logical :: debugging logical :: conv(nspec) c indexes integer :: ii,i3,i,j,k,ki,iter,ix,ir,kern,iii integer :: ip,nr1,is1,ip1,ip2,ip3 integer :: iintt integer :: icount integer :: imxc, istr external soln !c++mcb !---------------------------------------------------------------------- ccrit = 1.e-4 ncyc = 1 Cycle_loop: do icyc=1,ncyc xb(0) = 1. xbc(0) = 1. ! conv: density !** The solution is somewhat sensitive to the !** Order of the steady-state species in the input file xb(1:nspec) = spec(1:nspec) ! molecules/cm3 where( xb(1:nspec) < 1.e-10 ) xb(1:nspec) = 0. endwhere xbc(1:naq) = aspec(1:naq) write(*,*) 'in chem, species concentration at time =', time write(*,'(10e11.4)') (spec(i), i=1,8) ! iterate the solution !++mcb Do Gauss-Seidel method over two groups in this order, ! (1) nspec and clear air points, and (2) all species on cloudy points has_cldwtr: if(ncpts == 0) then xa(1:nspec) = xb(1:nspec) Iter_loop: do iter = 1,50 if(iter > 2) then more = any( .not. conv(1:nspec) ) else more = .true. endif if(more) then ! continue iterating over some gridpoints ! solve for gas phase species Gas_spc_loop: do ix=1,nspec ! loop over gas phase species gasSnk = 0. debugging = mod(time,dt_dbg) == 0. .and. ix==kso2 if( debugging ) write(56,*) 'Gas sinks ', time, iter !---------------------------------------------------------------------- ! set the total sinks for each gas phase species !---------------------------------------------------------------------- do ir = 1,icmap(ix,1) ! loop over number of sink reactions ip=(ir-1)*5 nr1=icmap(ix,2+ip+1) ! index of reaction rate is1=icmap(ix,2+ip+2) ! index of stoich. coeff. ip1=icmap(ix,2+ip+3) ! index of species #1 ip2=icmap(ix,2+ip+4) ! index of species #2 ip3=icmap(ix,2+ip+5) ! index of species #3 !---------------------------------------------------------------------- ! computation of the loss, where the loss is divided by the ! species concentration xb(i,ip3) !---------------------------------------------------------------------- gasSnk = gasSnk + stc(is1)*rate(nr1)*xb(ip1)*xb(ip2) if( debugging ) then write(56,'(2i4,4e12.5,4x,e12.5)') ir, nr1, rate(nr1), $ xb(ip1), xb(ip2), xb(ip3), $ stc(is1)*rate(nr1)*xb(ip1)*xb(ip2) endif end do if( debugging ) then write(56,*) 'Gas sources ', time, iter endif !---------------------------------------------------------------------- ! set the total sources for each gas phase species !---------------------------------------------------------------------- gasSrc = 0. do ir = 1,icmap(ix,2) ! loop over species (number of sources) ip=((ir-1)+icmap(ix,1))*5 nr1=icmap(ix,2+ip+1) ! index of reaction rate is1=icmap(ix,2+ip+2) ! index of stoich. coeff. ip1=icmap(ix,2+ip+3) ! index of species #1 ip2=icmap(ix,2+ip+4) ! index of species #2 ip3=icmap(ix,2+ip+5) ! index of species #3 !---------------------------------------------------------------------- ! compute the source !---------------------------------------------------------------------- gasSrc = gasSrc $ + stc(is1)*rate(nr1)*xb(ip1)*xb(ip2)*xb(ip3) if( debugging ) then write(56,'(2i4,4e12.5,4x,e12.5)') ir, nr1, rate(nr1), $ xb(ip1), xb(ip2), xb(ip3), $ stc(is1)*rate(nr1)*xb(ip1)*xb(ip2)*xb(ip3) endif end do !---------------------------------------------------------------------- ! computation new gas phase concentrations !---------------------------------------------------------------------- xnew = (xa(ix) + dtchem*gasSrc)/(1. + gasSnk*dtchem) conv(ix) = abs(xnew-xb(ix)) <= abs(xb(ix))*ccrit if(iter < 30) then xb(ix) = xnew else xtmp = xb(ix) xb(ix) = (xnew + xtmp) * 0.5 endif end do Gas_spc_loop else exit Iter_loop endif ! no more points end do Iter_loop if( iter > 50 ) then write(*,*) $ 'MORE THAN 50 ITERATIONS; ', $ count(.not.conv(1:nspec)),' not converged' endif elseif(ncpts > 0) then has_cldwtr !---------------------------------------------------------------------- ! gas and aqueous phase !---------------------------------------------------------------------- do ix = 1,nspec iy = gasmap(ix) xa(ix) = xb(ix) if(iy > 0) then xac(iy) = xbc(iy) ! define total concentration of species in aqueous phase ctot1(iy) = xa(ix) + xac(iy) ctot1t(iy) = ctot1(iy) endif end do Iter_loopa: do iter = 1,50 if(iter > 2) then more = any( .not. conv(1:nspec) ) else more = .true. endif if(more) then All_spc_loop: do ix = 1,nspec ! loop over species debugging = mod(time,dt_dbg) == 0. .and. ix==kso2 !---------------------------------------------------------------------- ! set the gas phase sinks for each species !---------------------------------------------------------------------- if( debugging ) write(56,*) 'Gas sinks ', time, iter gasSnk = 0. do ir = 1,icmap(ix,1) ! loop over number of sink reactions ip=(ir-1)*5 nr1=icmap(ix,2+ip+1) ! index of reaction rate is1=icmap(ix,2+ip+2) ! index of stoich. coeff. ip1=icmap(ix,2+ip+3) ! index of species #1 ip2=icmap(ix,2+ip+4) ! index of species #2 ip3=icmap(ix,2+ip+5) ! index of species #3 !---------------------------------------------------------------------- ! computation of the gas phase loss: ! the loss is divided by the species concentration xb(i,ip3) !---------------------------------------------------------------------- gasSnk = gasSnk + stc(is1)*rate(nr1)*xb(ip1)*xb(ip2) if( debugging ) then write(56,'(2i4,4e12.5,4x,e12.5)') ir, nr1, rate(nr1), $ xb(ip1), xb(ip2), xb(ip3), $ stc(is1)*rate(nr1)*xb(ip1)*xb(ip2) endif end do !---------------------------------------------------------------------- ! find the gas sources for each species !---------------------------------------------------------------------- if( debugging ) write(56,*) 'Gas sources ', time, iter gasSrc = 0. do ir = 1,icmap(ix,2) ! loop over species (number of sources) ip=((ir-1)+icmap(ix,1))*5 nr1=icmap(ix,2+ip+1) ! index of reaction rate is1=icmap(ix,2+ip+2) ! index of stoich. coeff. ip1=icmap(ix,2+ip+3) ! index of species #1 ip2=icmap(ix,2+ip+4) ! index of species #2 ip3=icmap(ix,2+ip+5) ! index of species #3 !---------------------------------------------------------------------- ! computation of the gas phase source !---------------------------------------------------------------------- gasSrc = gasSrc $ + stc(is1)*rate(nr1)*xb(ip1)*xb(ip2)*xb(ip3) if( debugging ) then write(56,'(2i4,4e12.5,4x,e12.5)') ir, nr1, rate(nr1), $ xb(ip1), xb(ip2), xb(ip3), $ stc(is1)*rate(nr1)*xb(ip1)*xb(ip2)*xb(ip3) endif end do aqSnk = 0. aqSrc = 0. iy = gasmap(ix) if(iy > 0) then !---------------------------------------------------------------------- ! set the cloud water sinks for each species !---------------------------------------------------------------------- if( debugging ) write(46,*) 'Aqueous sinks ', time, iter nsnka = icmapa(iy,1) do i = 1,nsnka ip = (i-1)*5 ir1=icmapa(iy,2+ip+1) ! index of reaction rate is1=icmapa(iy,2+ip+2) ! index of stoich. coeff. ip1=icmapa(iy,2+ip+3) ! index of species #1 ip2=icmapa(iy,2+ip+4) ! index of species #2 ip3=icmapa(iy,2+ip+5) ! index of species #3 aqSnk = aqSnk $ + stc(is1) * ratea(ir1) * xbc(ip1) * xbc(ip2) if( debugging ) then write(46,'(2i3,f5.2,5e11.4)') $ i, ir1, stc(is1), ratea(ir1), $ xbc(ip1), xbc(ip2), xbc(ip3), $ stc(is1)*ratea(ir1)*xbc(ip1)*xbc(ip2) endif end do if( debugging ) write(46,*) 'Aqueous sources ', time, iter !---------------------------------------------------------------------- ! set the cloud water sources for each species !---------------------------------------------------------------------- nsrca = icmapa(iy,2) do i = 1,nsrca ip = 5.*nsnka + (i-1)*5 ir1=icmapa(iy,2+ip+1) ! index of reaction rate is1=icmapa(iy,2+ip+2) ! index of stoich. coeff. ip1=icmapa(iy,2+ip+3) ! index of species #1 ip2=icmapa(iy,2+ip+4) ! index of species #2 ip3=icmapa(iy,2+ip+5) ! index of species #3 aqSrc = aqSrc + stc(is1) * ratea(ir1) * $ xbc(ip1) * xbc(ip2) * xbc(ip3) if( debugging ) then write(46,'(2i3,f5.2,5e11.4)') $ i, ir1, stc(is1), ratea(ir1), $ xbc(ip1), xbc(ip2), xbc(ip3), $ stc(is1)*ratea(ir1)*xbc(ip1)*xbc(ip2) endif end do endif !---------------------------------------------------------------------- ! compute new gas and aqueous phase concentrations !---------------------------------------------------------------------- if(iy > 0) then xold = xb(ix) + xbc(iy) else xold = xb(ix) endif call soln(lww, gasSrc, gasSnk, aqSrc, aqSnk, $ ktc, heff, phrat, $ xb, xbc, ctot1t, xa, xac, $ xbn, xbcn, ix, dtchem) if(iy > 0) then xnew = ctot1t(iy) else xnew = xbn endif conv(ix) = abs(xnew-xold) <= abs(xold)*ccrit if(iter < 30) then xb(ix) = xbn if(iy > 0) then xbc(iy) = xbcn endif else xtmp = xb(ix) xb(ix) = (xbn + xtmp) * 0.5 if(iy > 0) then xtmp = xbc(iy) xbc(iy) = (xbcn + xtmp) * 0.5 endif endif end do All_spc_loop else exit Iter_loopa endif ! no .more. points end do Iter_loopa if(iter > 50) then write(*,*) $ 'MORE THAN 50 cldy ITERATIONS; ', $ count(.not.conv(1:nspec)),' not converged' endif endif has_cldwtr spec(1:nspec) = xb(1:nspec) aspec(1:naq) = xbc(1:naq) end do Cycle_loop end subroutine chem subroutine soln(lww, srcg, snkg, srca, snka, $ ktc, heff, phrat, $ csp, aqc, ctot1t, csp1, aqc1, $ xbn, xbcn, n, dt) ! Solve for concentration at t(n+1) via iterated EB implicit none include 'chem_params.inc' include 'amap.com' ! Input variables integer, intent(in) :: n ! species index real, intent(in) :: lww ! LWC (cm3/cm3) real, intent(in) :: srcg, snkg ! gas phase source,sink real, intent(in) :: srca, snka ! aq phase source,sink real, intent(in) :: ktc(naq), heff(naq), phrat(naq) real, intent(in) :: csp(0:nspec), aqc(0:naq) real, intent(in) :: csp1(nspec), aqc1(naq) real, intent(in) :: dt ! time step (s) real, intent(inout) :: ctot1t(naq) ! Output variables real, intent(out) :: xbn, xbcn ! Local variables integer :: ik, ic, it, iy real :: ctotal, cliq, xing real :: snk, srcavg real :: srcgavg, srcaavg, snkgavg, snkaavg real :: todrop, otdrop real :: snkgtmp, srcgtmp, snkatmp, srcatmp real :: cg iy = gasmap(n) cg = csp(n) if(iy > 0) then ! Determine new concentration for species that interact with drops srcgavg = srcg + ktc(iy)*aqc(iy)/heff(iy) srcaavg = srca + ktc(iy)*lww*cg snkgavg = snkg + ktc(iy)*lww snkaavg = snka + ktc(iy)/heff(iy) xbn = (csp1(n) + dt*srcgavg)/(1. + snkgavg*dt) xbcn = (aqc1(iy) + dt*srcaavg)/(1. + snkaavg*dt) ctot1t(iy) = xbcn + xbn else ! Determine new concentration for species that don't interact with drops srcgavg = srcg snkgavg = snkg xbn = (csp1(n) + dt*srcgavg)/(1. + snkgavg*dt) endif end subroutine soln subroutine xinit(buf,iend,x) implicit none integer i,iend real buf(*) real x do i=1,iend buf(i)=x enddo end subroutine xinit subroutine iinit(buf,iend,ival) integer buf(*) do i=1,iend buf(i)=ival enddo end subroutine iinit SUBROUTINE OPNDAT( iucl, clfile ) CHARACTER*20 clfile C----------------------------------------------------------- C Open clfile, on error write error number and terminate. C----------------------------------------------------------- nend=index(clfile,' ') open(unit=iucl, file=clfile(1:nend-1),status='OLD', 1 access='sequential',form='formatted',iostat=ios) if (ios .ne. 0) then print 1000, ios stop end if Return 1000 Format(1x,'error opening data file ',i5) END