#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: CHEM_RATES ! ! !DESCRIPTION: routines that calculate rates. Also used if there is no ! chemistry (without_chemistry) to calculate Henry coefficients ! for wet deposition. !\ !\ ! !INTERFACE: ! MODULE CHEM_RATES ! ! !USES: ! use GO, only : gol, goPr, goErr IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: rates ! rate constants & Henry law #ifndef without_chemistry public :: calrates #endif ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'chem_rates' ! ! !REVISION HISTORY: ! 26 Mar 2010 - P. Le Sager - took off NO+XO2 three rates ! 15 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ contains !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: RATES ! ! !DESCRIPTION: ! ! calculation of look up tables rate constants/henry coefficients ! ! method ! ------ ! use known array of temperatures (integers) to calculate rate constants ! 3 body reactions are explicitly calculated in chemistry !\ !\ ! !INTERFACE: ! SUBROUTINE RATES( status ) ! ! !USES: ! use toolbox, only : zfarr use chem_param ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 26 Mar 2010 - P. Le Sager - added protex ! 16 Nov 2010 - P. Le Sager - bug fix : kh1 & kh2 are real now. ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ !BOC integer :: itemp, k, i real :: kh1, kh2, br real :: ztrec, zt3rec, temp ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/rates' ! start ! JEW : updated rates JPL(2006), incl. Evaluation Number 15 (March, 2006) ! do k=1,ntemp itemp=k+ntlow temp=float(itemp) ztrec=1./float(itemp) zt3rec=300./float(itemp) !JEW: changed to IUPAC 2013 rates_lut(knoo3,k)=zfarr(2.07e-12,-1400.,ztrec) rates_lut(kho2no,k)=zfarr(3.45e-12,270.,ztrec) !JEW: changed to IUPAC 2012 rates_lut(kno2oha,k)=3.2e-30*(temp/300.)**(-4.5) rates_lut(kno2ohb,k)=3.0e-11 ! JEW : reference IUPAC rates_lut(kohhno3a,k)=zfarr(2.41e-14,460.,ztrec) rates_lut(kohhno3b,k)=zfarr(6.51e-34,1335.,ztrec) rates_lut(kohhno3c,k)=zfarr(2.69e-17,2199.,ztrec) rates_lut(kno2o3,k)=zfarr(1.4e-13,-2470.,ztrec) rates_lut(knono3,k)=zfarr(1.8e-11,110.,ztrec) !JEW: changed to IUPAC 2012 rates_lut(kno2no3a,k)=3.6e-30*(temp/300.)**(-4.1) rates_lut(kno2no3b,k)=1.9e-12*(temp/300.)**0.2 rates_lut(kn2o5a,k)=1.3e-3*(temp/300.)**(-3.5) rates_lut(kn2o5b,k)=9.7e14*(temp/300.)**0.1 rates_lut(kohhono,k)=zfarr(2.5e-12,260.,ztrec) rates_lut(khno4oh,k)=zfarr(3.2e-12,690.,ztrec) rates_lut(kno2ho2a,k)=1.4e-31*(temp/300.)**(-3.1) rates_lut(kno2ho2b,k)=4.0e-12 ! independant of T rates_lut(khno4a,k)=zfarr(4.1e-5,-10650.,ztrec) rates_lut(khno4b,k)=zfarr(6.0e15,-11170.,ztrec) rates_lut(khonoa,k)=7.4e-31*(temp/300.)**(-2.4) rates_lut(khonob,k)=3.3e-11*(temp/300.)**(-0.3) !JEW: changed to JPL2006 rates_lut(kodm,k)=.2095*zfarr(3.3e-11,55.,ztrec)+ & .7808*zfarr(2.15e-11,110.,ztrec) !JEW: changed to JPL2006 rates_lut(kh2ood,k)=zfarr(1.63e-10,60.,ztrec) rates_lut(ko3ho2,k)=zfarr(1.0e-14,-490.,ztrec) rates_lut(ko3oh,k)=zfarr(1.7e-12,-940.,ztrec) !JEW: changed to JPL2006 rates_lut(khpoh,k)=1.8e-12 rates_lut(kfrmoh,k)=zfarr(5.5e-12,125.,ztrec) !JEW: changed to JPL2006 rates_lut(kch4oh,k)=zfarr(2.45e-12,-1775.,ztrec) ! TvN: JPL2015 gives slightly different exponents: ! KCOOHA: 1.4 -> 1.0 ! KCOOHC: -0.6 -> 0.0 rates_lut(kcooha,k)=5.9e-33*zt3rec**1.4 rates_lut(kcoohb,k)=1.1e-12*zt3rec**(-1.3) rates_lut(kcoohc,k)=1.5e-13*zt3rec**(-0.6) rates_lut(kcoohd,k)=2.1e9*zt3rec**(-6.1) !JEW: changed to JPL2006 rates_lut(kohmper,k)=zfarr(3.8e-12,200.,ztrec) rates_lut(kohrooh,k)=zfarr(3.01e-12,190.,ztrec) ! CB05 rates_lut(kho2oh,k)=zfarr(4.8e-11,250.,ztrec) rates_lut(kh2oh,k)=zfarr(2.8e-12,-1800.,ztrec) !JEW: changed to IUPAC 2009 br=1./(1.+498.*exp(-1160./temp)) rates_lut(kmo2ho2a,k)=zfarr(3.8e-13,750.,ztrec) rates_lut(kmo2ho2a,k)=rates_lut(kmo2ho2a,k)*(1.0-br) rates_lut(kmo2ho2b,k)=zfarr(3.8e-13,750.,ztrec) rates_lut(kmo2ho2b,k)=rates_lut(kmo2ho2b,k)*br rates_lut(kmo2no,k)=zfarr(2.3e-12,360.,ztrec) rates_lut(kmo2mo2,k)=zfarr(9.5e-14,390.,ztrec) !JEW: changed to JPL2006 ! TvN: It looks like the first two reactions below ! have updated coefficients in JPL2015: rates_lut(kho2ho2a,k)=zfarr(3.5e-13,430.,ztrec) rates_lut(kho2ho2b,k)=zfarr(1.7e-33,1000.,ztrec) rates_lut(kho2ho2c,k)=zfarr(1.4e-21,2200.,ztrec) rates_lut(kc41,k)=5.8e-16 rates_lut(kc46,k)=zfarr(8.1e-12,270.,ztrec) ! ! Reformulated as defined in IUPAC recommendations ! rates_lut(kc47a,k)=3.28e-28*(temp/300.)**(-6.87) rates_lut(kc47b,k)=1.125e-11*(temp/300.)**(-1.105) rates_lut(kc48a,k)=zfarr(1.1e-5,-10100.,ztrec) rates_lut(kc48b,k)=zfarr(1.9e17,-14100.,ztrec) !JEW: changed to JPL2006 rates_lut(kc49,k)=zfarr(2.9e-12,500.,ztrec) rates_lut(kc50,k)=zfarr(4.3e-13,1040.,ztrec) !------------------------------------------------------ rates_lut(kc52,k)=8.1e-13 rates_lut(kc53,k)=zfarr(1.e15,-8000.,ztrec) rates_lut(kc54,k)=1.6e3 !JEW: changed to JPL2006 rates_lut(kc57,k)=zfarr(5.4e-12,-610.,ztrec) !JEW: changed to JPL2006 rates_lut(kc58,k)=zfarr(8.5e-16,1520.,ztrec) rates_lut(kc59,k)=zfarr(4.6e-14,400.,ztrec) !JEW: changed to JPL2006 ! TvN: JPL2015 recommends different coefficients ! for OH+C2H4 (KC61): ! (1.e-28,-4.5) -> (1.1e-28,-3.5) ! (8.8e-12,-0.85) -> (8.4e-12,-1.75) rates_lut(kc61a,k)=1.e-28*(temp/300.)**(-4.5) rates_lut(kc61b,k)=8.8e-12*(temp/300.)**(-0.85) !JEW: changed to JPL2006 rates_lut(kc62,k)=zfarr(1.2e-14,-2630.,ztrec) !JEW: changed to IUPAC2004 rates_lut(kc73,k)=1.5e-11 ! IUPAC rates_lut(kc76,k)=zfarr(2.7e-11,390.,ztrec) ! IUPAC rates_lut(kc77,k)=zfarr(1.03e-14,-1995.,ztrec) ! IUPAC rates_lut(kc78,k)=zfarr(3.15e-12,-450.,ztrec) ! IUPAC !JEW: changed to JPL2006 rates_lut(kc79,k)=zfarr(2.6e-12,365.,ztrec) rates_lut(kc80,k)=zfarr(1.6e-12,-2200.,ztrec) rates_lut(kc81,k)=zfarr(2.6e-12,365.,ztrec) ! bug rates_lut(kc82,k)=zfarr(7.5e-13,700.,ztrec) ! CB05 rates_lut(kc83,k)=8.e-11 rates_lut(kc84,k)=zfarr(5.9e-13,-360.,ztrec) ! CB05 temp dep rates_lut(kc85,k)=zfarr(7.5e-13,700.,ztrec) ! CB05 rates_lut(KO3PO2,k)=6.0e-34*(temp/300.)**(-2.4) rates_lut(KO3PO3,k)=zfarr(8.0e-12,-2060.,ztrec) rates_lut(KXO2XO2N,k)=6.8e-14 rates_lut(KXO2N,k)=6.8e-14 ! sulfur and ammonia gas phase reactions ! the hynes et al. parameterisation is replaced by chin et al. 1996 !JEW: changed to JPL2006 rates_lut(kdmsoha,k)= 1.11e-11*exp(-240./temp) rates_lut(kdmsohb,k)=1.0e-39*exp(5820./temp) rates_lut(kdmsohc,k)=5.0e-30*exp(6280./temp) rates_lut(kdmsno3,k)=zfarr(1.9e-13,520.,ztrec)!at 298 1.01e-12 !JEW: changed to JPL2006 rates_lut(kso2oha,k)=3.3e-31*(temp/300.)**(-4.3) rates_lut(kso2ohb,k)= 1.6e-12*(temp/300.) ! ! fate of ammonia ! rates_lut(knh3oh,k)= zfarr(1.7e-12,-710.,ztrec)!1.56e-13 at 298k rates_lut(knh2oh,k)= 3.4e-11 rates_lut(knh2no,k)= zfarr(4.0e-12,450.,ztrec)!1.72e-11 rates_lut(knh2no2,k)= zfarr(2.1e-12,650.,ztrec)!1.86e-11 rates_lut(knh2ho2,k)= 3.4e-11 ! TvN: JPL2015 only says value is smaller than 6.0e-21 rates_lut(knh2o2,k)= 6.0e-21 rates_lut(knh2o3,k)= zfarr(4.3e-12,-930.,ztrec)!1.89e-13 at 298k ! ! duplicate rates for NH2O2 ! rates_lut(knh2o2o3,k)= rates_lut(knh2o3,k) rates_lut(knh2o2no,k)= rates_lut(knh2no,k) ! ! for higher organics rates_lut(kohmcho,k) = zfarr(4.4e-12,365.,ztrec) ! IUPAC rates_lut(kohmch2cho,k) = zfarr(4.9e-12,405.,ztrec) rates_lut(kno3mcho,k) = zfarr(1.4e-12,-1860.,ztrec) rates_lut(kno3mch2cho,k) = 6.4e-15 rates_lut(kohole,k) = zfarr(8.2e-12,610.,ztrec) ! IUPAC rates_lut(ko3ole,k) = 1.0e-17 ! Corrected to expression in Williams et al. (ACP,2013) !rates_lut(kno3ole,k) = zfarr(4.0e-14,-400.,ztrec) rates_lut(kno3ole,k) = zfarr(4.6e-14,400.,ztrec) ! ! the rates for additional BVOC's ! ! TvN: JPL2015 only gives two digits: 2.9e-12 rates_lut(kohch3oh,k) = zfarr(2.85e-12,-345.,ztrec) rates_lut(kohhcooh,k) = 4.0e-13 rates_lut(kohmcooh,k) = zfarr(4.2e-14,855.,ztrec) ! TvN: JPL2015 recommends different coefficients: ! (6.9e-12,-1000.) -> (7.66e-12,-1020.) ! (3.0e-12,20.) -> (3.35e-12,0.) ! (7.6e-12,-585.) -> (8.7e-12,-615.) rates_lut(kohc2h6,k) = zfarr(6.9e-12,-1000.,ztrec) rates_lut(kohethoh,k) = zfarr(3.0e-12,20.,ztrec) rates_lut(kohc3h8,k) = zfarr(7.6e-12,-585.,ztrec) ! TvN: JPL2015 recommends different coefficients ! (8.0e-27,-3.5) -> (4.6e-27,-4.) ! (3.0e-11,-1.0) -> (2.6e-11,-1.3) rates_lut(kohc3h6a,k) = 8.0e-27*(temp/300.)**(-3.5) rates_lut(kohc3h6b,k) = 3.0e-11*(temp/300.) **(-1.0) rates_lut(ko3c3h6,k) = zfarr(5.5e-15,-1880.,ztrec) ! IUPAC rates_lut(kno3c3h6,k) = zfarr(4.6e-13,-1155.,ztrec) ! IUPAC ! ! Terpenes ! rates_lut(kohterp,k) = zfarr(1.2e-11,440.,ztrec) ! IUPAC rates_lut(ko3terp,k) = zfarr(6.3e-16,-580.,ztrec) ! IUPAC rates_lut(kno3terp,k) = zfarr(1.2e-12,490.,ztrec) ! IUPAC ! ! Acetone ! rates_lut(kohaceta,k) = zfarr(8.8e-12,-1320.,ztrec) ! IUPAC rates_lut(kohacetb,k) = zfarr(1.7e-14,423.,ztrec) ! IUPAC rates_lut(kaco2ho2,k) = 1.0e-11 rates_lut(kaco2mo2,k) = 3.8e-12 ! IUPAC rates_lut(kaco2no,k) = 8.0e-12 rates_lut(kaco2xo2,k) = 6.8e-14 rates_lut(kxo2xo2n,k) = 6.8e-14 rates_lut(kxo2n,k) = 6.8e-14 ! ! MACR, MVK ! rates_lut(kohmvk,k)=zfarr(2.6e-12,610.,ztrec) rates_lut(kohmacr,k)=zfarr(8.0e-12,380.,ztrec) rates_lut(ko3mvk,k)=zfarr(8.5e-16,-1520.,ztrec) rates_lut(ko3macr,k)=zfarr(1.4e-15,-2100.,ztrec) ! ! Additional peroxy from C3 VOC ! rates_lut(kc3h7o2no,k)=zfarr(7.6e-12,-585.,ztrec) rates_lut(kc3h7o2ho2,k)=zfarr(7.5e-13,700.,ztrec) rates_lut(khypo2no,k)=zfarr(4.2e-12,180.,ztrec) rates_lut(khypo2ho2,k)=zfarr(7.5e-13,700.,ztrec) ! ! methyl peroxy nitrate ! rates_lut(kmo2no2a,k)=2.5e-30*(temp/300.)**(-5.5) rates_lut(kmo2no2b,k)=1.8e-11 rates_lut(kch3o2no2a,k)=9.0E-5*exp(-9690./temp) rates_lut(kch3o2no2b,k)=1.1e16*exp(-10560./temp) ! ! **** solubility Henry equilibrium ! HNO3/so4/nh4 just a very high number to take H and ! dissociation into account ! henry(:,k)=0. henry(ih2o2,k)=1.0e5*exp(6300*ztrec)*6.656e-10 henry(ihno3,k)=1e7 henry(ich3o2h,k)=310.*exp(5200*ztrec)*2.664e-8 henry(ich2o,k)=3000.*exp(7200*ztrec)*3.253e-11 henry(irooh,k)=340.*exp(6000.*ztrec)*1.821e-9 henry(ich2o,k)=henry(ich2o,k)*37 henry(ich3o2no2,k)=2.0*exp(4700.*ztrec)*exp(-4700.*(1/298.15)) henry(iorgntr,k)=zfarr(1.8e-6,6000.,ztrec) henry(iso4,k)=1.e7 henry(inh4,k)=1.e7 henry(imsa,k)=1.e7 henry(iso2,k) =1.2*exp(3120.*ztrec)*2.85e-5 henry(inh3,k) =75.0*exp(3400.*ztrec)*1.10e-5 henry(io3,k)=1.1e-2*exp(2300.*ztrec)*4.45e-4 ! JEW add two new scavenging rates for CH3COCHO and ALD2 ! need KH(eff) due to hydration steps for both species henry(imgly,k) = 3.2e4*48.6 kh1=17.*exp(5000.*ztrec)*exp(-5000.*(1/298.15)) kh2=13.*exp(5700.*ztrec)*exp(-5700.*(1/298.15)) henry(iald2,k) = ((kh1+kh2)/2.)*1.0246 henry(ihcooh,k) = 8900.*exp(6100.*ztrec)*exp(-6100.*(1/298.15)) henry(ich3oh,k) = 220.0*exp(5200.*ztrec)*2.66e-8 henry(imcooh,k)= 4.1e3*exp(6300.*ztrec)*exp(-6300.*(1/298.15)) henry(iethoh,k) = 190.0*exp(6600.*ztrec)*exp(-6600.*(1/298.15)) henry(iacet,k) = 35.*exp(3800.*ztrec)*exp(-3800.*(1/298.15)) end do !k temperature loop ! ! marked tracers: ! henry(io3s,:) = henry(io3,:) !ok status = 0 end subroutine rates !EOC #ifndef without_chemistry !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: calrates ! ! !DESCRIPTION: ! calculate rate constants using lookup table rates_lut ! calculate third bodies ! calculate heterogeneous removal on aerosols ! ! External: CALCHET ! !\ !\ ! !INTERFACE: ! SUBROUTINE CALRATES( region, level, is, js, rjx, rr, reff_cloud, ye, dt_chem, sad_cld_out, sad_ice_out, sad_aer_out, status ) ! ! !USES: ! use global_data, only : region_dat, mass_dat use Dims, only : CheckShape, idate, jm use chem_param use binas, only : Avog, pi use tm5_distgrid, only : dgrid, Get_DistGrid use reaction_data, only : nreac #ifdef with_m7 use m7_data, only : rw_mode, dens_mode use mo_aero_m7, only : nmod, cmr2ras #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region, level, is, js real, intent(in) :: reff_cloud(is:,js:) real, intent(in) :: dt_chem ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: rr(is:,js:,:) ! reaction rates integer, intent(out) :: status real, intent(out) :: sad_cld_out(is:,js:), sad_ice_out(is:,js:), sad_aer_out(is:,js:) ! ! !INPUT/OUTPUT PARAMETERS: ! real, intent(inout) :: rjx(is:,js:,:) ! photolysis rate real, intent(inout) :: ye(is:,js:,:) ! extra 2D fields ! ! !REVISION HISTORY: ! 26 Mar 2010 - P. Le Sager - added protex ! 16 Jun 2011 - P. Le Sager - bug fix from JEW : g_n2o5 values ! 15 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! 12 Jun 2012 - P. Le Sager - fix O3 concentration ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/calc_rates' real, dimension(:,:,:,:), pointer :: rm ! tracer mass ! help variables integer :: itemp,i,j, i1, j1, i2, j2, imode, iaer,n,latindex real :: tr, temp, wv, airp, rx1, rx2, rx3 real :: dum, h2ox, aird, air_vol, o2, o3 real :: x1, x2, xice, xliq ! ! cloud chemistry of n2o5 real :: dg, kt_liq_n2o5, kt_ice_n2o5, kt_liq_ho2, kt_ice_ho2 real :: g_n2o5_l_temp real :: r_ice, r_cloud ! cm real :: sad_ice, sad_cloud, iwc, lwc, o3_p ! Aerosol heterogeneous chemistry real :: kt_aer_n2o5, kt_aer_no3, kt_aer_ho2 real :: sad_aer, sad_aert, aer_conc, aer_dens, aer_rad real :: smr_aer real, parameter :: sad_factor = 4.*pi*1.e-3 ! 4*pi as prefactor for area, 1e-3 to convert air_vol to cm3 ! Standard aerosol properties density and radius for NO3_A, MSA,NH4, and -if not with_m7- SO4 real, parameter :: aer_dens_ref = 1.7 ! assumed particle density [gr/cm3] real, parameter :: aer_rad_ref_so4_LRH = 0.18e-4 ! assumed particle radius [cm] real, parameter :: aer_rad_ref_no3_LRH = 0.2e-4 ! assumed particle radius [cm] real, parameter :: aer_rad_ref_so4_HRH = 0.25e-4 ! assumed particle radius accounting for abs. H2O[cm] real, parameter :: aer_rad_ref_no3_HRH = 0.27e-4 ! assumed particle radius accounting for abs. H2O[cm] ! Parameters to describe subgrid scale mixing real :: zcc,nn real, parameter :: M_HO2 = (1+2*16.)/Avog *1e-3 ! HO2 mass, kg/molec real, parameter :: M_N2O5 = (2*14+5*16.)/Avog *1e-3 ! N2O5 mass, kg/molec real, parameter :: M_NO3 = (14+3*16.)/Avog *1e-3 ! no3 mass, kg/molec real :: C_Thermal_N2O5, C_Thermal_HO2,C_Thermal_NO3, gamma, ttemp, factor, g_n2o5_aer_so4 real :: aer_rad_ref, tot, no3frac, so4frac, r_no3, r_so4 real :: zsgs_mix_n2o5 real,dimension(180) :: h2_surface ! --- begin -------------------------------- CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) call CheckShape( (/i2-i1+1, j2-j1+1, nj /), shape(rjx), status ) IF_ERROR_RETURN(status = 1) call CheckShape( (/i2-i1+1, j2-j1+1, nreac /), shape(rr ), status ) IF_ERROR_RETURN(status = 1) call CheckShape( (/i2-i1+1, j2-j1+1, n_extra/), shape(ye ), status ) IF_ERROR_RETURN(status = 1) rm => mass_dat(region)%rm h2_surface(1:180)=h2_lat(1:180,idate(2)) !$OMP parallel & !$OMP default (none) & !$OMP shared (level, jsr, jer, isr, ier, region, ye, rates_lut, rr, rjx, rm, fscale) & !$OMP private (i, j, jss, jee, rx3, temp, itemp, airp, h2ox, aird, & !$OMP tr, wv, o2, o3, o3_p, rx1, rx2) rx3=0.6 rr(i1:i2,j1:j2,1:nreac)=0.0 do j=j1,j2 do i=i1,i2 temp=ye(i,j,i_temp) itemp=nint(temp-float(ntlow)) itemp=min(max(itemp,1),ntemp) !limit temperatures in loop-up table range airp=ye(i,j,i_pres) ! ! Calculation of relative humidity [%] ! Richardson's approximation for water vapor pressure ! should be calculated in subroutine rates ! h2ox = ye(i,j,ih2on) aird = ye(i,j,iairn) tr=1.-373.15/temp wv=exp((((-.1299*tr-.6445)*tr-1.976)*tr+13.3185)*tr) ye(i,j,irh)=h2ox*temp/(1013.25*wv*7.24291e16) ye(i,j,irh) = max(min(ye(i,j,irh),100.0),0.0) !limit rh between 0-100% o2=0.209476*aird o3 = rm(i,j,level,io3) / ye(i,j,iairm) * aird * fscale(io3) !kg ----> #/cm3 ! !**** calculate three body reaction rates nn=0.75-(1.27*LOG10(0.41)) ! rx1=rates_lut(KNO2OHA,itemp)*aird rx2=rates_lut(KNO2OHB,itemp) rr(i,j,kno2oh)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.41)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(KOHHNO3C,itemp) rx2=rates_lut(KOHHNO3B,itemp)*aird rr(i,j,kohhno3)=rates_lut(KOHHNO3A,itemp)+rx1*rx2/(rx1+rx2) ! ! reformulated following IUPAC 2013 ! nn=0.75-(1.27*LOG10(0.35)) rx1=rates_lut(KNO2NO3A,itemp)*aird rx2=rates_lut(KNO2NO3B,itemp) rr(i,j,kno2no3)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.35)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(KN2O5A,itemp)*aird*exp(-11000/temp) rx2=rates_lut(KN2O5B,itemp)*exp(-11080/temp) rr(i,j,kn2o5)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.35)/(1+(log10(RX1/RX2)/nn)**2)) nn=0.75-(1.27*LOG10(0.4)) rx1=rates_lut(KNO2HO2A,itemp)*aird rx2=rates_lut(KNO2HO2B,itemp) rr(i,j,kno2ho2)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.4)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(KHNO4A,itemp)*aird rx2=rates_lut(KHNO4B,itemp) ! !=((4.1e-5*exp(-10650/T))*M*(6.0e15*exp(-11170/T)))/((4.1e-5*exp(-10650/T))*M+(6.0e15*exp(-11170/T)))*10^(log10(0.4)/(1+(log10((4.1e-5*exp(-10650/T))*M/(6.0e15*exp(-11170/T)))/(0.75-1.27*log10(0.4)))^2)) ! rr(i,j,khno4m)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.4)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(khonoa,itemp)*aird rx2=rates_lut(khonob,itemp) rr(i,j,khono)=rx1/(1.+rx1/rx2)*0.81**(1./(1.+log10(rx1/rx2)**2)) nn=0.75-(1.27*LOG10(0.36)) rx1=rates_lut(KMO2NO2A,itemp)*aird rx2=rates_lut(KMO2NO2B,itemp) rr(i,j,kmo2no2)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.36)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(KCH3O2NO2A,itemp)*aird rx2=rates_lut(KCH3O2NO2B,itemp) rr(i,j,kch3o2no2m)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.36)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(kcooha,itemp)*aird rx2=rates_lut(kcoohb,itemp) rr(i,j,kcooh)=rx1/(1.+rx1/rx2)*rx3**(1./(1.+log10(rx1/rx2)**2)) !JEW: now add the second term for CO + OH rx1=rates_lut(kcoohc,itemp) rx2=rates_lut(kcoohd,itemp)/aird rr(i,j,kcooh)=rr(i,j,kcooh)+(rx1/(1.+rx1/rx2)*rx3**(1./(1.+log10(rx1/rx2)**2))) rx1=rates_lut(KC61A,itemp)*aird rx2=rates_lut(KC61B,itemp) rr(i,j,kc61)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**2)) rx1=rates_lut(kohc3h6a,itemp)*aird rx2=rates_lut(kohc3h6b,itemp) rr(i,j,kohc3h6)=rx1/(1.+rx1/rx2)*0.5**(1./(1.+log10(rx1/rx2)**2)) ! ! photolysis rates and 2 body reactions ! rr(i,j,knoo3)=rates_lut(KNOO3,itemp) rr(i,j,kho2no)=rates_lut(KHO2NO,itemp) rr(i,j,kmo2no)=rates_lut(KMO2NO,itemp) rr(i,j,kno2o3)=rates_lut(KNO2O3,itemp) rr(i,j,knono3)=rates_lut(KNONO3,itemp) rr(i,j,khno4oh)=rates_lut(KHNO4OH,itemp) rr(i,j,kohhono)=rates_lut(KOHHONO,itemp) rr(i,j,kodm)=rates_lut(KODM,itemp) rr(i,j,kh2ood)=rates_lut(KH2OOD,itemp) rr(i,j,ko3ho2)=rates_lut(KO3HO2,itemp) !old formulation rr(i,j,kcooh)=1.5E-13+9E-14*airp/101325. rr(i,j,ko3oh)=rates_lut(KO3OH,itemp) rr(i,j,khpoh)=rates_lut(KHPOH,itemp) rr(i,j,kfrmoh)=rates_lut(KFRMOH,itemp) rr(i,j,kch4oh)=rates_lut(KCH4OH,itemp) latindex=j*180/jm(region) rr(i,j,kh2oh)=rates_lut(kh2oh,itemp)*h2_surface(latindex)*1.0e-9*aird !h2 latitudinal constraint JEW update rr(i,j,kohmper)=rates_lut(KOHMPER,itemp) rr(i,j,kohrooh)=rates_lut(KOHROOH,itemp) rr(i,j,kmo2ho2a)=rates_lut(KMO2HO2a,itemp) rr(i,j,kmo2ho2b)=rates_lut(KMO2HO2b,itemp) rr(i,j,kmo2mo2)=rates_lut(KMO2MO2,itemp) rr(i,j,kho2oh)=rates_lut(KHO2OH,itemp) rr(i,j,kho2ho2)=(rates_lut(KHO2HO2A,itemp) + & rates_lut(KHO2HO2B,itemp)*aird)*(1.+rates_lut(KHO2HO2C,itemp)*h2ox) rr(i,j,kc41)=rates_lut(kc41,itemp) ! JEW for ALD take the average of the C2 and C3 rate co-efficients; measurements suggest ! that CH3CHO only comprises 50% of [higher aldehydes] - grossman et al, JGR, 2003. rr(i,j,kc43)=(rates_lut(kohmcho,itemp)+rates_lut(kohmch2cho,itemp))/2. rr(i,j,kc44)=(rates_lut(kno3mcho,itemp)+ rates_lut(kno3mch2cho,itemp))/2. rr(i,j,kc46)=rates_lut(KC46,itemp) ! ! PAN now reformulated (2014) ! nn=0.75-(1.27*LOG10(0.3)) rx1=rates_lut(KC47A,itemp)*aird rx2=rates_lut(KC47B,itemp) ! rr(i,j,kc47)=RX1/(1.+RX1/RX2)*0.3**(1./(1.+LOG10(RX1/RX2)**2)) ! rr(i,j,kc47)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.3)/(1+(log10(RX1/RX2)/nn)**2)) rx1=rates_lut(KC48A,itemp)*aird rx2=rates_lut(KC48B,itemp) ! rr(i,j,kc48)=RX1/(1.+RX1/RX2)*0.3**(1./(1.+LOG10(RX1/RX2)**2) ! rr(i,j,kc48)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.3)/(1+(log10(RX1/RX2)/nn)**2)) rr(i,j,kc49)=rates_lut(KC49,itemp) rr(i,j,kc50)=rates_lut(KC50,itemp) rr(i,j,kc52)=rates_lut(KC52,itemp) rr(i,j,kc53)=rates_lut(KC53,itemp) rr(i,j,kc54)=rates_lut(KC54,itemp) ! JEW updated average for OLE oxidation reactions rr(i,j,kc57)=(rates_lut(kohmvk,itemp)+rates_lut(kohole,itemp)+rates_lut(kohmacr,itemp))/3. rr(i,j,kc58)=(rates_lut(ko3mvk,itemp)+rates_lut(ko3ole,itemp)+rates_lut(ko3macr,itemp))/3. rr(i,j,kc59)=(rates_lut(kno3ole,itemp)+6.0e-16+3.5e-15)/3. rr(i,j,kc62)=rates_lut(KC62,itemp) rr(i,j,kc73)=rates_lut(KC73,itemp) rr(i,j,kc76)=rates_lut(KC76,itemp) rr(i,j,kc77)=rates_lut(KC77,itemp) rr(i,j,kc78)=rates_lut(KC78,itemp) ! JEW use rates in cb05 for XO2 reactions rr(i,j,kc79)=rates_lut(KC79,itemp) rr(i,j,kc80)=rates_lut(KC80,itemp) rr(i,j,kc81)=rates_lut(KC81,itemp) rr(i,j,kc82)=rates_lut(KC82,itemp) rr(i,j,kc83)=rates_lut(KC83,itemp) rr(i,j,kc84)=rates_lut(kc84,itemp) rr(i,j,kc85)=rates_lut(kc85,itemp) #ifndef without_photolysis rjx(i,j,jo3d)=rjx(i,j,jo3d)*rr(i,j,kh2ood)*h2ox / & ( rr(i,j,kodm)*aird + rr(i,j,kh2ood)*h2ox ) ! this is fraction rjo3d=>OH #endif RX1=rates_lut(kso2oha,itemp)*aird RX2=rates_lut(kso2ohb,itemp) rr(i,j,kso2oh)=RX1/(1.+RX1/RX2)*0.6**(1./(1.+LOG10(RX1/RX2)**2)) ! ! dmsoha => so2 ! dmsohb => 0.75 SO2 + 0.25 MSA ! rr(i,j,kdmsoha)=rates_lut(kdmsoha,itemp) rr(i,j,kdmsohb)=rates_lut(kdmsohb,itemp)*o2/ & (1.+rates_lut(kdmsohc,itemp)*o2) rr(i,j,kdmsno3)=rates_lut(kdmsno3,itemp) ! ammonia chemistry rr(i,j,knh3oh)=rates_lut(knh3oh,itemp) rr(i,j,knh2oh)=rates_lut(knh2oh,itemp) rr(i,j,knh2no)=rates_lut(knh2no,itemp) rr(i,j,knh2no2)=rates_lut(knh2no2,itemp) rr(i,j,knh2ho2)=rates_lut(knh2ho2,itemp) rr(i,j,knh2o2)=rates_lut(knh2o2,itemp)*o2 rr(i,j,knh2o3)=rates_lut(knh2o3,itemp) rr(i,j,knh2o2no)=rates_lut(knh2no,itemp) rr(i,j,knh2o2o3)=rates_lut(knh2o3,itemp) rr(i,j,knh2o2ho2)=rates_lut(knh2ho2,itemp) rr(i,j,krn222)=2.11e-6 ! s-1 decay time Rn222 to Pb210 rr(i,j,ko3po2) = rates_lut(ko3po2,itemp)*aird rr(i,j,ko3po3) = rates_lut(ko3po3,itemp) rr(i,j,kxo2xo2n) = 6.8e-14 rr(i,j,kxo2n) = 6.8e-14 ! ! O3P in molecules cm3 ! ! JEW: Reformulated June 2012 ! ! Runaway O3 above 50 hPa due to missing stratopsheric chemistry ! therefore use pressure as an index ! o3_p=0. if(airp/100 > 50. .and. airp/100. < 350.) o3_p=(2*rjx(i,j,jo2)*o2) ! ! [O2] not used in the ebi ! rjx(i,j,jo2)=(o3_p/(rr(i,j,ko3po2)*o2+rr(i,j,ko3po3)*o3)) * o2 * rr(i,j,ko3po2) rr(i,j,ko3po3)=rr(i,j,ko3po3)*o3_p ! ! additional biogenic reactions rr(i,j,kohhcooh)=4.0e-13 rr(i,j,kohch3oh)=rates_lut(kohch3oh,itemp) ! additional no3 peroxy reactions rr(i,j,kno3ho2)=4.0e-12 rr(i,j,kno3mo2)=1.2e-12 rr(i,j,kno3c2o3)=4.0e-12 rr(i,j,kno3xo2)=2.5e-12 ! Zaveri and Peters ! additional C2 compounds rr(i,j,kohmcooh)=rates_lut(kohmcooh,itemp) rr(i,j,kohc2h6)=rates_lut(kohc2h6,itemp) rr(i,j,kohethoh)=rates_lut(kohethoh,itemp) rr(i,j,kohc3h8)=rates_lut(kohc3h8,itemp) rr(i,j,ko3c3h6)=rates_lut(ko3c3h6,itemp) rr(i,j,kno3c3h6)=rates_lut(kno3c3h6,itemp) ! ! TERPENE reactions ! rr(i,j,kohterp) = rates_lut(kohterp,itemp) rr(i,j,ko3terp) = rates_lut(ko3terp,itemp) rr(i,j,kno3terp) = rates_lut(kno3terp,itemp) ! ! ISPD reactions ! rr(i,j,kohispd)=(rates_lut(kohmvk,itemp)+rates_lut(kohmacr,itemp))/2 rr(i,j,ko3ispd)=(rates_lut(ko3mvk,itemp)+rates_lut(ko3macr,itemp))/2 rr(i,j,kno3ispd)=(6.0e-16+3.4e-15)/2 ! ! Radon decay rr(i,j,krn222)=2.11e-6 ! s-1 decay time Rn222 to Pb210 ! ! acetone rr(i,j,kohacet)=rates_lut(kohaceta,itemp)+rates_lut(kohacetb,itemp) rr(i,j,kaco2ho2)=1.0e-11 rr(i,j,kaco2mo2)=3.8e-12 rr(i,j,kaco2no)=8.0e-12 rr(i,j,kaco2xo2)=6.8e-14 ! CB05 peroxy-loss rr(i,j,kxo2xo2n)=6.8e-14 rr(i,j,kxo2n)=6.8e-14 ! ! rr(i,j,kc3h7o2no)=rates_lut(kc3h7o2no,itemp) rr(i,j,kc3h7o2ho2)=rates_lut(kc3h7o2ho2,itemp) rr(i,j,khypo2no)=rates_lut(khypo2no,itemp) rr(i,j,khypo2ho2)=rates_lut(khypo2ho2,itemp) ! ! end do end do !_lat/lon loop !$OMP END PARALLEL ! ! heterogeneous reaction of N2O5 and H2O -> 2 HNO3 on cloud and aerosol ! included in gas phase chemistry. ! VH, August-December 2013: Methodology updated, and extended with reactions for NO3 and HO2. ! rw_mode and dens_mode contain the average radius [m] and density [kg/m3] for each modes, ! as computed in subroutine m7_averageproperties (pm6rp and prhop) ! !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, ye, rr, reff_cloud,level, & !$OMP g_ho2,g_n2o5,g_no3, sad_cld_out, sad_ice_out,sad_aer_out, & !$OMP mode_nm, mode_tracers, rm,dt_chem, m_n2o5, m_ho2, m_no3) & !$OMP private (i, j, jss, jee, airp, dg, aird, xliq, xice, & !$OMP temp, sad_ice,sad_cloud,iwc,lwc,r_ice, zcc,& !$OMP r_cloud, sad_aer,sad_aert, aer_conc, air_vol, & !$OMP g_n2o5_l_temp,kt_aer_n2o5, kt_aer_ho2,kt_aer_no3, & !$OMP kt_liq_n2o5, kt_ice_n2o5, kt_liq_ho2, kt_ice_ho2 ) & #ifdef with_m7 !$OMP shared (rw_mode, dens_mode) & !$OMP private (aer_rad, aer_dens) #endif do j = j1, j2 do i=i1, i2 sad_cloud = 0. sad_ice = 0. aird = ye(i,j,iairn) ! #/cm3 ! typically the optical reff is somewhat larger than the physical r by 1-2um ! therefore downsize reff by 2.0uM for droplets 9-13uM and 1.0 for those between ! 6-9um ! http://www-das.uwyo.edu/~geerts/cwx/notes/chap08/moist_cloud.html ! set cloud droplet radius [cm] if(reff_cloud(i,j)>=9.0) then r_cloud = (reff_cloud(i,j)-2.0)/1e4 elseif(reff_cloud(i,j)>=6.0) then r_cloud = (reff_cloud(i,j)-1.0)/1e4 elseif(reff_cloud(i,j)>= 4.0) then r_cloud = (reff_cloud(i,j)-0.25)/1e4 else r_cloud = 4.0e-4 ! ensure a minimum - sometime reff_cloud is zero, i.e. not defined... endif ! SAD for ice particles ! JEW: The surface area density for ice particles in now linked to ! the IWC by the parameterization of Heymsfield and McFarquar (1996) ! and the effective radii from Fu (1996) ! VH the computation of sad_ice and r_ice is optimized, taking account of units. ! VH compute SAD representative for cloudy part only (scale with cloud cover) r_ice=5e-3 ! adopt as the default [cm] if(ye(i,j,iciwc)>1.0e-15 .and. ye(i,j,icc) > 1e-5) then ! convert iwc from units [kg/kg] to [gr/m3] iwc=(ye(i,j,iciwc)*aird*xmair*1e6/avog) / ye(i,j,icc) sad_ice=1.0e-4*iwc**0.9 ! [cm2/cm3] ! calculate the r_eff using the relationship of Fu (1996) ! r_ice=(1.73205/(3*0.917))*((iwc/1e6)/sad_ice) filling constants: ! ! Bug fix by VH: corrected 3/(4*rho_ice) ! r_ice=0.8179*((iwc/1e6)/sad_ice) ! [cm] ! the value adopted in von Kuhlmann and Lawrence is too low sad_ice=10*sad_ice ! [cm2/cm3] endif if(ye(i,j,iclwc)>1.0e-15 .and. ye(i,j,icc)>1e-5 ) then ! lwc convert from units [kg/kg] to [kg/m3] lwc=( ye(i,j,iclwc)*aird*xmair*1e3/avog ) / ye(i,j,icc) lwc=lwc*1e-3 ! convert [kg/m3] to [gr/cm3] !VH no_cloud=lwc/sphere_vol !VH sad_cloud=4*pi*r_cloud**2*no_cloud sad_cloud=3.*lwc/(r_cloud) endif ! Assume that loss on cloud particles is not dominating the full loss budget which can ! remove the full trace gas concentration within one time step, within grid cel that is ! partly covered with cloud zcc = min(0.99,ye(i,j,icc)) ! ! We use the original formulation in dentener and crutzen ! of : kt=(r/Dg + 4/vgamma)-1 * A(cm2/cm3) airp=ye(i,j,i_pres) dg=0.1*1e5/airp !simple approximation for diffusion coeff. [cm2/s] !n2o5... (see IUPAC) temp=ye(i,j,i_temp) C_Thermal_N2O5 = SQRT (8*1.38E-23 *temp / (PI*M_N2O5 ) ) * 1e2 ! m/s -> cm/sec C_Thermal_HO2 = SQRT (8*1.38E-23 *temp / (PI*M_HO2 ) ) * 1e2 ! m/s -> cm/sec C_Thermal_NO3 = SQRT (8*1.38E-23 *temp / (PI*M_NO3 ) ) * 1e2 ! m/s -> cm/sec !===================================================================================== ! ! Full parameterization of gamma N2O5 variability of SO4= aerosols from Evans and Jacob, 2005 ! Implemented by JEW July 2014 ! ! RH dependence from Kane et al., Heterogeneous uptake of gaseous N2O5 by (NH4)2SO4 ! NH4HSO4 and H2SO4 aerosols, J. Phys. Chem A, 2001, 105, 6465-6470 ! !===================================================================================== gamma=2.79e-4 + ye(i,j,irh)*(1.3e-4 + ye(i,j,irh)*(-3.43e-6 + ye(i,j,irh)*7.52e-8)) ttemp=MAX(temp, 282.) factor=10.**(-2.0-4.0e-2*(ttemp-294.0))/0.01 g_n2o5_aer_so4=gamma*factor if(zcc > 0.0) then g_n2o5_l_temp = 2.7e-5*exp(1800./temp) kt_liq_n2o5=1./((r_cloud/dg)+(4./(C_Thermal_N2O5*g_n2o5_l_temp))) kt_ice_n2o5=1./((r_ice/dg) +(4./(C_Thermal_N2O5*g_n2o5_i))) ! ! Scaling factor accounting for subgrid-scale mixing within time step ! Assume a lagrangian time scale ZDT_LAG of 3h (i.e.time scale for mixing) ! ZSGS_MIX=0 if CC=0, and 1 if CC=1; ZSGS_MIX is lower than CC for 0 volume (units) 1e-3[cm3] sad_aert=0. #ifdef with_m7 do imode=1,nmod aer_conc = 0. ! Total concentration of aerosol within this mode. Changed to units: [kg]/1e-3[cm3] -> [gr/cm3] do iaer=1,mode_nm(imode) aer_conc=aer_conc+rm(i,j,level,mode_tracers(iaer,imode)) / air_vol enddo if (aer_conc .gt.1e-15) then !>>> TvN ! Avoid using dens_mode, since it does not account for the presence of ! ammonium nitrate and associated water, and MSA- in the soluble accumulation mode. ! Since these components are taken into account in rw_mode, ! the use of dens_mode would introduce inconsistencies. !aer_dens=max(1.0,dens_mode(region,imode)%d3(i,j,level)) *1.e-3 ! particle density [kg/m3] -> [gr/cm3] aer_rad =max(1e-10,rw_mode (region,imode)%d3(i,j,level)) *1e2 ! particle median radius [m] -> [cm] ! To calculate the surface area of each mode, ! first convert the number median radius to the surface mean radius: smr_aer = aer_rad * cmr2ras(imode) ! The suface area can then be calculated as the mean surface area per particle ! times the number of particles in the mode: !sad_aer=max(3.0 * aer_conc / ( aer_dens * aer_rad ),1.e-15) ! units [cm2/cm3] sad_aer = rm(i,j,level,mode_start(imode)) * sad_factor * smr_aer**2 / air_vol sad_aert=sad_aert+sad_aer ! total aerosol SAD for evaluation purposes ! In principle one should average the full expression for ! the rate constant over the lognormal size distributions. ! In practice, an effective value for the radius is used ! in the expression for kt (see Jacob, Atm. Env., 2000), ! in our case the number median radius. ! <<< TvN ! Ensure that g_xxx > 0. !!! kt_aer_n2o5=1./((aer_rad/dg)+(4./(C_Thermal_N2O5 * g_n2o5(imode)))) kt_aer_no3 =1./((aer_rad/dg)+(4./(C_Thermal_NO3 * g_no3 (imode)))) kt_aer_ho2 =1./((aer_rad/dg)+(4./(C_Thermal_HO2 * g_ho2 (imode)))) ! Fill reaction rates rr(i,j,kn2o5_aer)=rr(i,j,kn2o5_aer)+kt_aer_n2o5 * sad_aer rr(i,j,kno3_aer) =rr(i,j,kno3_aer) +kt_aer_no3 * sad_aer rr(i,j,kho2_aer) =rr(i,j,kho2_aer) +kt_aer_ho2 * sad_aer ! aer_rad_max(imode)=max(aer_rad_max(imode),aer_rad) endif enddo !>>> TvN ! Commented the part below, ! because rw_mode already accounts for the presence of NO3_a, NH4, MSA, and aerosol water ! in the soluble accumulation mode (see chemistry.F90). ! ! simple treatment for NO3_a, NH4 and MSA !imode = nmod+1 !aer_conc=(rm(i,j,level,inh4)+rm(i,j,level,ino3_a)+rm(i,j,level,imsa)) / air_vol !if (aer_conc.gt.1e-15) then ! FIXME - need to CHECK AER_RAD_REF - first iteration ! if(ye(i,j,irh) > 69.0) then ! aer_rad_ref = ( rm(i,j,level,ino3_a)*aer_rad_ref_no3_HRH + & ! rm(i,j,level,imsa )*aer_rad_ref_so4_HRH ) /& ! (rm(i,j,level,ino3_a) + rm(i,j,level,imsa)) ! else ! aer_rad_ref = ( rm(i,j,level,ino3_a)*aer_rad_ref_no3_LRH + & ! rm(i,j,level,imsa )*aer_rad_ref_so4_LRH ) /& ! (rm(i,j,level,ino3_a) + rm(i,j,level,imsa)) ! endif ! sad_aer=max(3.0 * aer_conc / ( aer_dens_ref * aer_rad_ref ),1.e-15) ! units [cm2/cm3] ! sad_aert=sad_aert+sad_aer ! total aerosol SAD for evaluation purposes ! Ensure that g_xxx > 0. !!! ! kt_aer_n2o5=1./((aer_rad_ref/dg)+(4./(C_Thermal_N2O5 * g_n2o5(imode)))) ! kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(C_Thermal_NO3 * g_no3 (imode)))) ! kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(C_Thermal_HO2 * g_ho2 (imode)))) ! Fill reaction rates ! rr(i,j,kn2o5_aer)=rr(i,j,kn2o5_aer)+kt_aer_n2o5 * sad_aer ! rr(i,j,kno3_aer) =rr(i,j,kno3_aer) +kt_aer_no3 * sad_aer ! rr(i,j,kho2_aer) =rr(i,j,kho2_aer) +kt_aer_ho2 * sad_aer !endif !<<< TvN #else ! Uptake on on available aerosol: SO4, NO3, NH4, and MSA. ! Total concentration of sulfate aerosol (include others) [gr/cm3] aer_conc=(rm(i,j,level,iso4)+rm(i,j,level,inh4)+rm(i,j,level,ino3_a)+rm(i,j,level,imsa)) / air_vol if (aer_conc.gt.1e-20) then no3frac=rm(i,j,level,ino3_a)/air_vol !IGNORE NH4? so4frac=(rm(i,j,level,iso4)+rm(i,j,level,imsa))/air_vol ! ! JEW (2014) : Introduce a crude fix for swelling of particles at high RH ! NH4+ exists in either SO4-- or NO3- state, so don't double count ! for aerosol volumes. ! if(ye(i,j,irh) < 70.0) then r_no3=(no3frac/aer_conc)*aer_rad_ref_no3_LRH r_so4=(so4frac/aer_conc)*aer_rad_ref_so4_LRH aer_rad_ref=r_no3+r_so4 endif if(ye(i,j,irh) > 69.0) then r_no3=(no3frac/aer_conc)*aer_rad_ref_no3_HRH r_so4=(so4frac/aer_conc)*aer_rad_ref_so4_HRH aer_rad_ref=r_no3+r_so4 endif sad_aer = 3.0 * aer_conc / ( aer_dens_ref * aer_rad_ref ) ! units [cm2/cm3] sad_aert=sad_aer ! Ensure that g_xxx > 0. !!! kt_aer_n2o5=1./((r_no3/dg)+(4./(C_Thermal_N2O5 * g_n2o5_aer))) kt_aer_n2o5=kt_aer_n2o5+(1./((r_so4/dg)+(4./(C_Thermal_N2O5 * g_n2o5_aer_so4)))) kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(C_Thermal_NO3 * g_no3_aer))) kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(C_Thermal_HO2 * g_ho2_aer))) ! Fill reaction rates rr(i,j,kn2o5_aer)=kt_aer_n2o5 * sad_aer rr(i,j,kno3_aer) =kt_aer_no3 * sad_aer rr(i,j,kho2_aer) =kt_aer_ho2 * sad_aer endif #endif /* M7 */ ! output sad_aer_out(i,j)=sad_aert end do end do !$OMP END PARALLEL nullify(rm) status = 0 end subroutine calrates !EOC #endif /*CHEM*/ END MODULE CHEM_RATES