subroutine aerochem_cmu(water,tempk,press,lwc_c,cph,cNV, & mxspec_c,con,bdnl,convfac, & t00,dtaq,dtaer,aero_flag,ldark) c c----CAMx v7Beta6 190902 c c AEROCHEM_CMU drives the CMU multi-section aerosol modules for CAMx c c Copyright 1996 - 2018 c Ramboll c c Modifications: c 01/30/2002 (tmg); do not set prevsp to true because each grid cell c will be different c 03/22/2002 (tmg); do not update inorganic gases after soap module c because they are not changed inside it c 04/15/2002 (tmg); changed molecular weights in conversion of c organic gases c 11/25/2002 (tmg); added include file other.inc to track whether c aerosol or aqueous module is called c 12/05/2002 (tmg); added countcell to track which grid cell is c current (for debugging purposes) c 03/09/03 (bkoo) removed code related with lsoap c (SOAP has been merged with inorganic aerosol module) c 12/26/2006 (bkoo) AQ-CHEM now called each chemistry timestep c 12/29/2006 (bkoo) added species mapping for the updated SOA scheme c 01/10/2007 (bkoo) Added species mapping for CB05/SAPRC99 c 05/30/2007 (cae) Skip call to RAQCHEM under high acid conditions. c 08/23/2013 (bkoo) Replaced hard-wired MW's with wtmol & emw c 03/18/2014 (bkoo) Revised for benzene SOA c 08/25/2016 (bkoo) Updated for new SOAP c 09/16/2016 (bkoo) Revised for RAQCHEM updates c 11/22/2017 (bkoo) Updated for revised WTMOL c c Input arguments: c water cell water vapor (ppm) c tempk cell temperature (K) c press cell pressure (mb) c lwc_c cloud water content (g/m3) c cph cloud water pH c cNV NV product concentrations (ppm) c mxspec_c Max number of species c con species concentrations (ppm, ug/m3) c bdnl lower bound concentrations (ppm, ug/m3) c convfac conversion factor: umol/m3 = ppm * convfac c t00 model time (HHMM) c dtaq aqueous chemistry time step (hours) c dtaer general PM chemistry time step (hours) c aero_flag 1 - time for AQ; 2 - time for general PM c ldark darkness flag c c Output arguments: c con species concentrations (ppm, ug/m3) c c Routines called: c GET_PARAM c RAQCHEM c AQDIST c AQCHEM c WDIAMETER c EQPARTO c DDIAMETER c NEWDIST c STEP c AERCHEM c c Called by: c CHEMDRIV c include 'dynamic.inc' include 'camx.prm' include 'aerpar.inc' include 'droppar.inc' include 'camx_aero.inc' include 'camx_aero_cmu.inc' include 'soap.inc' c real*4 cNV(NNV) real*4 con(mxspec_c+1) real*4 bdnl(mxspec_c+1) real*4 lwc_c,tempk,ev,es,water,press,convfac,t00 real*4 dtaq,dtaer,hour,tmin,t1_min,dt_min real*4 prs,qwatr,rhumid real*4 gas(ngas_aq), aerosol(nsect,naers) real*8 t0,t1 real*8 q(ntotal) integer modeaero,aero_flag logical ldark c c-----Variables for NVOL -> SOPA & SOPB c real rlambda,accom,frqtot,wdia real qmas(nsec),qnum(nsec),frq(nsec) c c-----Variables for RADM c real r_gas(15),r_aer(10) real pres_pa,dt_sec,cw_kgm3 real r_sum(7),cut(nsect+1) real dnit,dchlo,dammo,dsulf,dsoac real totpcl,totph2o,totpso4 real s_neg,s_pos real cph real r_1,r_2 integer r_idx(5) c c-----Variables for iron and manganese fractions in dust and flag for using dust fractions c-----instead of standard background values for iron and manganese c-----THESE MUST BE CONSISTENT WITH THOSE DEFINED IN AEROCHEM_CF.F logical :: dust_femn = .FALSE. real, parameter :: frac_Fe = 0.03355, frac_MN = 0.00115 real fe_iii, mn_ii, fe_sol, mn_sol, fe_mmr, mn_mmr c c-----Variables for ISORROPIA (for adjustment after RADM) c real*8 wi(5),wt(5),wt_save(5) real*8 cntrl(2) real*8 gasis(3),aerliq(12),aersld(9),other(6) character*15 scasi c c-----Entry point c modeaero = 0 c c-----Pass some variables to aerosol module common blocks c temp = tempk pres = dble(press/1013.25) prs = press/1013.25 c hour = aint(t00/100.) tmin = t00 - 100.*hour if (tmin.ge.59.99 .and. tmin.le.60.01) tmin = 60. if (tmin.ge.60.0) then tmin = tmin - 60. hour = hour + 1. endif c t1_min = 60.*hour + tmin ! end time (min) for AQCHEM dt_min = dtaq*60. ! delta t (min) for AQCHEM t1 = DBLE(t1_min*60.) ! end time (sec) for AEROCHEM dt = DBLE(dtaer*3600.) ! delta t (sec) for AEROCHEM t0 = t1 - dt ! start time (sec) for AEROCHEM tcom = t0 ! common t (sec) for AEROCHEM c c-----Calculate RH c qwatr = 1.e-6*water*18./28.8 ev = qwatr*press/(qwatr+eps) es = e0*exp((lv/rv)*(1./273.-1./tempk)) rhumid = amin1(0.99,ev/es) rh = rhumid c if (lfrst) then do i = 1,nsecp1 dsec(i) = dsec_c(i) dsecf(i) = dsecf_c(i) enddo endif c c-----Condense NVOL to SOPA/SOPB completely (non-volatile) c rlambda=0.065 ! this should be consistent with that in eqpart.f accom = delta(icg1) ! assuming the same accomodation coefficient for all CGs frqtot = 0.0 do isec = 1, nsec qmas(isec) = -con(kph2o_c+isec-1) ! dry do isp = ngas_c+1, nspec_c, nsec qmas(isec) = qmas(isec) + con(isp+isec-1) enddo qnum(isec) = qmas(isec) / dsec(isec)**3 qmas(isec) = qmas(isec) + con(kph2o_c+isec-1) ! wet wdia = (qmas(isec)/qnum(isec))**(0.33333) frq(isec) = qnum(isec) * wdia / (1. + rlambda / (accom*wdia)) frqtot = frqtot + frq(isec) enddo do isec = 1, nsec frq(isec) = frq(isec) / frqtot con(ksopa_c+isec-1) = con(ksopa_c+isec-1) + & cNV(1) * convfac * emw(ksopa) * frq(isec) ! anthro con(ksopb_c+isec-1) = con(ksopb_c+isec-1) + & cNV(2) * convfac * emw(ksopb) * frq(isec) ! bio enddo c if (lwc_c.ge.aqcwmin .and. tempk.ge.aqtamin) then modeaero = 1 ! flag to indicate aqueous chem module is called c c-----RADM aqueous chemistry c if ( chaq.eq.'RADM' ) then pres_pa = 100.*press dt_sec = dtaq*3600. cw_kgm3 = lwc_c/1000. c c-----Gas species mapping c r_gas(1) = con(kso2_c)*1.e-6 r_gas(2) = con(khno3_c)*1.e-6 r_gas(3) = con(kn2o5_c)*1.e-6 r_gas(4) = co2*1.e-6 r_gas(5) = con(knh3_c)*1.e-6 r_gas(6) = con(khpo_c)*1.e-6 r_gas(7) = con(ko3_c)*1.e-6 if (kfoa_c.lt.nspec_c+1) then r_gas(8) = con(kfoa_c)*1.e-6 else r_gas(8) = foa*1.e-6 endif if (kmhp_c.lt.nspec_c+1 .or. kohp_c.lt.nspec_c+1) then r_1 = con(kmhp_c) + con(kohp_c) r_gas(9) = r_1*1.e-6 r_1 = con(kmhp_c) / r_1 else r_gas(9) = mhp*1.e-6 endif if (kpaa_c.lt.nspec_c+1 .or. kopa_c.lt.nspec_c+1) then r_2 = con(kpaa_c) + con(kopa_c) r_gas(10) = r_2*1.e-6 r_2 = con(kpaa_c) / r_2 else r_gas(10) = paa*1.e-6 endif r_gas(11) = con(kh2so4_c)*1.e-6 if (kgly_c.lt.nspec_c+1) then r_gas(12) = con(kgly_c)*1.e-6 else r_gas(12) = gly*1.e-6 endif if (kmgly_c.lt.nspec_c+1) then r_gas(13) = con(kmgly_c)*1.e-6 else r_gas(13) = mgly*1.e-6 endif if (kglyd_c.lt.nspec_c+1) then r_gas(14) = con(kglyd_c)*1.e-6 else r_gas(14) = glyd*1.e-6 endif r_gas(15) = con(koh_c)*1.e-6 c c-----Aerosol species mapping c r_sum = 0.0 do ksec = 1,nsect r_sum(1) = r_sum(1) + con(kpso4_c+(ksec-1)) r_sum(2) = r_sum(2) + con(kpnh4_c+(ksec-1)) r_sum(3) = r_sum(3) + con(kpno3_c+(ksec-1)) r_sum(4) = r_sum(4) + con(kna_c +(ksec-1)) r_sum(5) = r_sum(5) + con(kpcl_c +(ksec-1)) r_sum(6) = r_sum(6) + con(kcrst_c+(ksec-1)) r_sum(7) = r_sum(7) + con(ksoac_c+(ksec-1)) enddo totpcl = r_sum(5) r_aer(1) = (r_sum(1)/wtmol(1)/convfac)*1.e-6 r_aer(2) = (r_sum(2)/wtmol(2)/convfac)*1.e-6 r_aer(3) = (r_sum(3)/wtmol(3)/convfac)*1.e-6 r_aer(4) = (caco3/(wtmol(8)+wtmol(7))/convfac)*1.e-6 r_aer(5) = (mgco3/(wtmol(9)+wtmol(7))/convfac)*1.e-6 if (r_sum(4).gt.r_sum(5)) then r_aer(6) = (r_sum(5)/wtmol(4)/convfac)*1.e-6 else r_aer(6) = (r_sum(4)/wtmol(5)/convfac)*1.e-6 endif c c-----Mods for Fe, Mn as dust components c if (dust_femn) then c c-----Assign fraction partitioning of FE(III) and MN(II) c if ( ldark ) then fe_iii = 0.9 ! Night time, GS 01july2011 else fe_iii = 0.1 ! Day time, GS 01july2011 endif mn_ii = 1.0 ! Same for day and night, GS 01July2011 c c-----Assign solubility of Fe and Mn c fe_sol = 0.1 ! GS 01July2011 mn_sol = 0.5 ! GS 28July2011 fe_mmr = (frac_Fe*r_sum(6)/wtmol(10)/convfac)*1.e-6 mn_mmr = (frac_MN*r_sum(6)/wtmol(11)/convfac)*1.e-6 r_aer(7) = fe_iii * fe_sol * fe_mmr ! GS 01July2011 r_aer(8) = mn_ii * mn_sol * mn_mmr ! GS 01July2011 else ! default CAMx treatment r_aer(7) = (a3fe/wtmol(10)/convfac)*1.e-6 r_aer(8) = (b2mn/wtmol(11)/convfac)*1.e-6 endif r_aer(9) = (potcl/(wtmol(6)+wtmol(4))/convfac)*1.e-6 c c-----Biogenic SOA tracer c r_aer(10) = (r_sum(7)/wtmol(12)/convfac)*1.e-6 c c-----Call RADM c if (r_aer(1).lt.5.0e-7 .and. r_aer(3).lt.1.5e-6) then ! check high acid conditions call get_param(igrdchm,ichm,jchm,kchm,iout,idiag) call raqchem(tempk,pres_pa,dt_sec,cw_kgm3,r_gas,r_aer, & cph,idiag,iout,igrdchm,ichm,jchm,kchm,wtmol(12)) c c-----Call ISORROPIA for ammonium/nitrate/chloride/water if AQ only c Since RADM-AQ puts all H2SO4 into aerosol phase ISORROPIA won't change c H2SO4/sulfate partitioning c if ( aero_flag.eq.1 ) then cntrl(1) = 0.d0 ! forward problem cntrl(2) = 0.d0 ! solids and liquid allowed wi(1) = r_sum(4)/wtmol(5)*1.e-6 ! total sodium wi(2) = (r_gas(11) + r_aer(1))*convfac ! total sulfate wi(3) = (r_gas(5) + r_aer(2))*convfac ! total ammonium wi(4) = (r_gas(2) + r_aer(3))*convfac ! total nitrate wi(5) = (con(khcl_c)*convfac + r_sum(5)/wtmol(4))*1.e-6 ! total chloride do is = 1, 5 wt_save(is) = wi(is) enddo call isoropia(wi,rh,temp,cntrl,wt,gasis,aerliq,aersld, & scasi,other) r_gas(5) = gasis(1) /convfac ! NH3 (mol/mol) r_aer(2) = (wt_save(3)-gasis(1))/convfac ! ammonium (mol/mol) r_gas(2) = gasis(2) /convfac ! HNO3 (mol/mol) r_aer(3) = (wt_save(4)-gasis(2))/convfac ! nitrate con(khcl_c) = gasis(3)/convfac*1.e6 ! HCL (ppm) con(khcl_c) = amax1(con(khcl_c),bdnl(khcl_c)) totpcl = (wt_save(5)-gasis(3))*wtmol(4)*1.e6 ! chloride (ug/m3) totpcl = amax1(totpcl ,0.0) totph2o = aerliq(8)*18.*1.e6 ! aerosol water (ug/m3) totph2o = amax1(totph2o,0.0) totpso4 = wt_save(2)*wtmol(1)*1.e6 ! sulfate (ug/m3) totpso4 = amax1(totpso4,1.e-12) endif c c-----Map gas back to con c con(kso2_c) = amax1(r_gas(1)*1.e6,bdnl(kso2_c)) ! SO2 (ppm) con(khno3_c) = amax1(r_gas(2)*1.e6,bdnl(khno3_c)) con(kn2o5_c) = amax1(r_gas(3)*1.e6,bdnl(kn2o5_c)) ! N2O5 gas (ppm) con(knh3_c) = amax1(r_gas(5)*1.e6,bdnl(knh3_c)) con(khpo_c) = amax1(r_gas(6)*1.e6,bdnl(khpo_c)) ! H2O2 (ppm) con(ko3_c) = amax1(r_gas(7)*1.e6,bdnl(ko3_c)) ! O3 (ppm) con(kfoa_c) = amax1(r_gas(8)*1.e6,bdnl(kfoa_c)) con(kmhp_c) = amax1(r_gas(9)*1.e6*r_1,bdnl(kmhp_c)) con(kohp_c) = amax1(r_gas(9)*1.e6*(1.-r_1),bdnl(kohp_c)) con(kpaa_c) = amax1(r_gas(10)*1.e6*r_2,bdnl(kpaa_c)) con(kopa_c) = amax1(r_gas(10)*1.e6*(1.-r_2),bdnl(kopa_c)) con(kh2so4_c)= amax1(r_gas(11)*1.e6,bdnl(kh2so4_c)) if (kgly_c.lt.nspec_c+1) then con(kgly_c) = amax1(r_gas(12)*1.e6,bdnl(kgly_c)) ! GLY(ppm) endif if (kmgly_c.lt.nspec_c+1) then con(kmgly_c) = amax1(r_gas(13)*1.e6,bdnl(kmgly_c)) ! MGLY(ppm) endif if (kglyd_c.lt.nspec_c+1) then con(kglyd_c) = amax1(r_gas(14)*1.e6,bdnl(kglyd_c)) ! GLYD(ppm) endif c c-----Calculate the differences c dsulf = amax1(r_aer(1)*1.e6*convfac*wtmol(1),0.0) - r_sum(1) dammo = amax1(r_aer(2)*1.e6*convfac*wtmol(2),0.0) - r_sum(2) dnit = amax1(r_aer(3)*1.e6*convfac*wtmol(3),0.0) - r_sum(3) dchlo = totpcl - r_sum(5) dsoac = amax1(r_aer(10)*1.e6*convfac*wtmol(12),0.0) - r_sum(7) c write(*,*)'RADM-change: ',dsulf,dnit,dammo c c-----Calculate fdist c do ksec = 1,nsect+1 cut(ksec) = SNGL(dsecf(ksec)) enddo call aqdist(nsect,cut,fdist,fdist2) c c-----Add the differences to the corresponding section c do ksec = 1,nsect con(kpso4_c+(ksec-1))=con(kpso4_c+(ksec-1))+dsulf*fdist(ksec) con(kpnh4_c+(ksec-1))=con(kpnh4_c+(ksec-1))+dammo*fdist(ksec) con(kpno3_c+(ksec-1))=con(kpno3_c+(ksec-1))+dnit *fdist(ksec) con(kpcl_c +(ksec-1))=con(kpcl_c +(ksec-1))+dchlo*fdist(ksec) con(ksoac_c+(ksec-1))=con(ksoac_c+(ksec-1))+dsoac*fdist(ksec) enddo c c-----Adjust mass to avoid negative concentration sections c r_idx(1) = kpso4_c r_idx(2) = kpnh4_c r_idx(3) = kpno3_c r_idx(4) = kpcl_c r_idx(5) = ksoac_c do n = 1,5 s_neg = 0.0 s_pos = 0.0 do ksec = 1,nsect if (con(r_idx(n)+(ksec-1)) .lt. 0.0) then s_neg = s_neg - con(r_idx(n)+(ksec-1)) con(r_idx(n)+(ksec-1)) = 0.0 elseif (con(r_idx(n)+(ksec-1)) .gt. 1.e-12) then s_pos = s_pos + con(r_idx(n)+(ksec-1)) endif enddo do ksec = 1,nsect if (con(r_idx(n)+(ksec-1)) .gt. 1.e-12) & con(r_idx(n)+(ksec-1)) = & amax1((1.-s_neg/s_pos)*con(r_idx(n)+(ksec-1)),0.) enddo enddo c c-----Assign new aerosol water to the new total sulfate distribution if AQ only c if ( aero_flag.eq.1 ) then do ksec = 1,nsect con(kph2o_c+(ksec-1)) = totph2o * con(kpso4_c+(ksec-1)) & / totpso4 enddo modeaero = 0 endif endif ! check high acid conditions c c-----Reset non-existing species concentrations c con(nspec_c+1) = 0. c c-----VSRM aqueous chemistry c else do k = 1,ngas_aq gas(k) = 0.0 enddo do k = 1,naers do ks = 1,nsect aerosol(ks,k) = 0.0 enddo enddo if (con(kh2so4_c) .gt. 0.0) then tot_sulf = 0.0 do k = 1,nsect tot_sulf = tot_sulf+con(kpso4_c+(k-1)) enddo if (tot_sulf .gt. 0.0) then do k = 1,nsect sfac = con(kpso4_c+(k-1))/tot_sulf con(kpso4_c+(k-1)) = con(kpso4_c+(k-1)) + & sfac*con(kh2so4_c)*convfac*96. enddo con(kh2so4_c) = 0.0 endif endif c c-----Map con to gas and aerosol gas in ppm; aerosol in ugr/m3 c gas(nga) = con(knh3_c) ! NH3(g) in ppm gas(ngn) = con(khno3_c) ! HNO3(g) in ppm gas(ngc) = con(khcl_c) ! HCl (g) in ppm gas(ngso2) = con(kso2_c) ! SO2(g) in ppm gas(ngh2o2) = con(kh2o2_c) ! H2O2(g) in ppm gas(nghcho) = con(kform_c) ! HCHO(g) in ppm gas(nghcooh) = 0.1*con(kh2o2_c) ! HCOOH(g) in ppm = 0.1*H2O2 gas(nghno2) = con(khono_c) ! HNO2(g) in ppm gas(ngo3) = con(ko3_c) ! O3 in ppm gas(ngoh) = con(koh_c) ! OH in ppm gas(ngho2) = con(kho2_c) ! HO2 in ppm gas(ngno3) = con(kno3_c) ! NO3 in ppm gas(ngno) = con(kno_c) ! NO(g) in ppm gas(ngno2) = con(kno2_c) ! NO2(g) in ppm gas(ngpan) = con(kpan_c) ! PAN(g) in ppm gas(ngch3o2h) = 0.2*con(kh2o2_c) ! CH3OOH(g) in ppm = 0.2*H2O2 gas(ngch3o2) = 1.0e-6 ! CH3O2(g) in ppm gas(ngch3oh) = 1.0e-3 ! CH3OH(g) in ppm = 1 ppb gas(ngch3co3h) = 0.05*con(kh2o2_c) ! CH3C(O)OOH(g) in ppm = 0.05*H2O2 c do knsec=1,nsect aerosol(knsec,naw) = con(kph2o_c+(knsec-1)) ! water aerosol(knsec,naa) = con(kpnh4_c+(knsec-1)) ! ammonium aerosol(knsec,na4) = con(kpso4_c+(knsec-1)) ! sulfate aerosol(knsec,nan) = con(kpno3_c+(knsec-1)) ! nitrate aerosol(knsec,nas) = con(kna_c+(knsec-1)) ! sodium aerosol(knsec,nac) = con(kpcl_c+(knsec-1)) ! chloride aerosol(knsec,nae) = con(kpec_c+(knsec-1)) ! elemental carbon aerosol(knsec,nao) = con(kpoa_c+(knsec-1)) ! primary organics aerosol(knsec,nar) = con(kcrst_c+(knsec-1)) ! crustal aerosol(knsec,nahso5) = 0.0 aerosol(knsec,nahmsa) = 0.0 enddo c call aqchem(gas,aerosol,rhumid,prs,tempk,lwc_c,t1_min, & dt_min,dsecf,ierr) c c-----CAMx doesn't carry HMSA or HSO5 so we put their mass into sulfate (NA4) c do knsec=1,nsect aerosol(knsec,na4) = aerosol(knsec,na4) + & 96./111.*aerosol(knsec,nahmsa) + & 96./113.*aerosol(knsec,nahso5) enddo c c-----Map gas and aerosol back con c con(knh3_c) = amax1(gas(nga) ,bdnl(knh3_c)) con(khno3_c) = amax1(gas(ngn) ,bdnl(khno3_c)) con(khcl_c) = amax1(gas(ngc) ,bdnl(khcl_c)) con(kso2_c) = amax1(gas(ngso2) ,bdnl(kso2_c)) con(kh2o2_c) = amax1(gas(ngh2o2),bdnl(kh2o2_c)) con(kform_c) = amax1(gas(nghcho),bdnl(kform_c)) con(khono_c) = amax1(gas(nghno2),bdnl(khono_c)) con(ko3_c) = amax1(gas(ngo3) ,bdnl(ko3_c)) con(kno3_c) = amax1(gas(ngno3) ,bdnl(kno3_c)) con(kno_c) = amax1(gas(ngno) ,bdnl(kno_c)) con(kno2_c) = amax1(gas(ngno2) ,bdnl(kno2_c)) con(kpan_c) = amax1(gas(ngpan) ,bdnl(kpan_c)) c do knsec=1,nsect con(kph2o_c+(knsec-1)) = aerosol(knsec,naw) con(kpnh4_c+(knsec-1)) = aerosol(knsec,naa) con(kpso4_c+(knsec-1)) = aerosol(knsec,na4) con(kpno3_c+(knsec-1)) = aerosol(knsec,nan) con(kna_c+(knsec-1)) = aerosol(knsec,nas) con(kpcl_c+(knsec-1)) = aerosol(knsec,nac) con(kpoa_c+(knsec-1)) = aerosol(knsec,nao) con(kpec_c+(knsec-1)) = aerosol(knsec,nae) con(kcrst_c+(knsec-1)) = aerosol(knsec,nar) enddo if ( aero_flag.eq.1 ) modeaero = 0 ! if AQ only endif endif c c-----RADM or VSRM -> call AER even if the aqueous module is called c OVSR -> call SOAP alone if the aqueous module is called c if ( aero_flag.eq.2 ) then if ( chaq.eq.'OVSR' .and. modeaero.eq.1 ) then modeaero = 1 else modeaero = 2 endif endif c c-----If neither AQCHEM nor AERCHEM is called we don't call soap c if (modeaero.ne.0) then do kk = 1,ntotal q(kk) = 0.d0 enddo c c-----Map con to q [ug] gas in ppm c do knsec=1,nsec q((knsec-1)*nsp+kh2o)=con(kph2o_c+(knsec-1)) q((knsec-1)*nsp+kna)=con(kna_c+(knsec-1)) q((knsec-1)*nsp+kso4)=con(kpso4_c+(knsec-1)) & / wtmol(1) * emw(kso4) q((knsec-1)*nsp+kno3)=con(kpno3_c+(knsec-1)) & / wtmol(3) * emw(kno3) q((knsec-1)*nsp+knh4)=con(kpnh4_c+(knsec-1)) & / wtmol(2) * emw(knh4) q((knsec-1)*nsp+kcl)=con(kpcl_c+(knsec-1)) & / wtmol(4) * emw(kcl) q((knsec-1)*nsp+ksoa1)=con(ksoa1_c+(knsec-1)) q((knsec-1)*nsp+ksoa2)=con(ksoa2_c+(knsec-1)) q((knsec-1)*nsp+ksoa3)=con(ksoa3_c+(knsec-1)) q((knsec-1)*nsp+ksoa4)=con(ksoa4_c+(knsec-1)) q((knsec-1)*nsp+ksopa)=con(ksopa_c+(knsec-1)) q((knsec-1)*nsp+ksopb)=con(ksopb_c+(knsec-1)) q((knsec-1)*nsp+kpom)=con(kpoa_c+(knsec-1)) q((knsec-1)*nsp+kec)=con(kpec_c+(knsec-1)) q((knsec-1)*nsp+kcrus)=con(kcrst_c+(knsec-1)) enddo c c-----Use MW instead of 100 and actual pressure and temperature for c conversion to ppm for organic gases (tmg, 04/15/02) c now organic gases are given in ppm - bkoo (08/25/03) c q(naer+inh3) = con(knh3_c) q(naer+ihno3) = con(khno3_c) q(naer+ih2so4) = con(kh2so4_c) q(naer+ihcl) = con(khcl_c) q(naer+icg1) = con(kcg1_c) q(naer+icg2) = con(kcg2_c) q(naer+icg3) = con(kcg3_c) q(naer+icg4) = con(kcg4_c) c c-----Call SOAP only c if (modeaero.eq.1) then do i = 1,nsec qt(i) = 0.d0 do ki = 2,nsp ! use dry basis qt(i) = qt(i)+q((i-1)*nsp+ki) ! total mass per section (ug/m3) enddo qn(i) = qt(i)/dsec(i)**3 ! calculate 0th moment (number of particles) enddo ! if density = 1 g/cm3 units are particles/cm3 call wdiameter(q) ! wet diameter ntotalx = 0 ! not used nsecx = 0 ! not used ntotalx2 = ntotal nsecx2 = nsec call eqparto(t1,q) ! equilibrium organic aerosol partitioning call ddiameter(q) ! dry diameter call newdist(t1,q) ! size distribution mapping call step(nsec,q) ! calculate water in each section c c-----Call SOAP + AER c else call aerchem(chaero,q,t0,t1,lfrst,ierr) endif c c-----Map q back to con c do knsec=1,nsec con(kph2o_c+(knsec-1))=q((knsec-1)*nsp+kh2o) con(kna_c +(knsec-1))=q((knsec-1)*nsp+kna) con(kpso4_c+(knsec-1))=q((knsec-1)*nsp+kso4) & * wtmol(1) / emw(kso4) con(kpno3_c+(knsec-1))=q((knsec-1)*nsp+kno3) & * wtmol(3) / emw(kno3) con(kpnh4_c+(knsec-1))=q((knsec-1)*nsp+knh4) & * wtmol(2) / emw(knh4) con(kpcl_c +(knsec-1))=q((knsec-1)*nsp+kcl) & * wtmol(4) / emw(kcl) con(ksoa1_c+(knsec-1))=q((knsec-1)*nsp+ksoa1) con(ksoa2_c+(knsec-1))=q((knsec-1)*nsp+ksoa2) con(ksoa3_c+(knsec-1))=q((knsec-1)*nsp+ksoa3) con(ksoa4_c+(knsec-1))=q((knsec-1)*nsp+ksoa4) con(ksopa_c+(knsec-1))=q((knsec-1)*nsp+ksopa) con(ksopb_c+(knsec-1))=q((knsec-1)*nsp+ksopb) con(kpoa_c +(knsec-1))=q((knsec-1)*nsp+kpom) con(kpec_c +(knsec-1))=q((knsec-1)*nsp+kec) con(kcrst_c+(knsec-1))=q((knsec-1)*nsp+kcrus) enddo c con(knh3_c) =amax1(SNGL(q(naer+inh3)), bdnl(knh3_c)) con(khno3_c) =amax1(SNGL(q(naer+ihno3)), bdnl(khno3_c)) con(kh2so4_c)=amax1(SNGL(q(naer+ih2so4)),bdnl(kh2so4_c)) con(khcl_c) =amax1(SNGL(q(naer+ihcl)), bdnl(khcl_c)) con(kcg1_c) = amax1(SNGL(q(naer+icg1)), bdnl(kcg1_c)) con(kcg2_c) = amax1(SNGL(q(naer+icg2)), bdnl(kcg2_c)) con(kcg3_c) = amax1(SNGL(q(naer+icg3)), bdnl(kcg3_c)) con(kcg4_c) = amax1(SNGL(q(naer+icg4)), bdnl(kcg4_c)) c endif c lfrst = .false. return end