#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 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 JPL2006 rates_lut(knoo3,k)=zfarr(3.e-12,-1500.,ztrec) !JEW: changed to JPL2006 rates_lut(kho2no,k)=zfarr(3.5e-12,250.,ztrec) !JEW: changed to JPL2006 rates_lut(kno2oha,k)=1.8e-30*zt3rec**3.0 rates_lut(kno2ohb,k)=2.8e-11 rates_lut(kohhno3a,k)=zfarr(2.41e-14,460.,ztrec) !!!wp!!! new ravi rates_lut(kohhno3b,k)=zfarr(6.51e-34,1335.,ztrec) !!!wp!!! new ravi rates_lut(kohhno3c,k)=zfarr(2.69e-17,2199.,ztrec) !!!wp!!! new ravi rates_lut(kno2o3,k)=zfarr(1.2e-13,-2450.,ztrec) rates_lut(knono3,k)=zfarr(1.5e-11,170.,ztrec) !JEW: changed to JPL2006 rates_lut(kno2no3a,k)=2.0e-30*zt3rec**4.4 rates_lut(kno2no3b,k)=1.4e-12*zt3rec**0.7 !JEW: changed to JPL2006 rates_lut(kn2o5,k)=zfarr(2.7e-27,11000.,ztrec) rates_lut(khno4oh,k)=zfarr(1.3e-12,380.,ztrec) !JEW: changed to JPL2006 rates_lut(kno2ho2a,k)=2.0e-31*zt3rec**3.4 rates_lut(kno2ho2b,k)=2.9e-12*zt3rec**1.1 rates_lut(khno4m,k)=zfarr(2.1e-27,10900.,ztrec) !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 !JEW: changed to JPL2006 rates_lut(kfrmoh,k)=zfarr(5.5e-12,125.,ztrec) !JEW: changed to JPL2006 rates_lut(kch4oh,k)=zfarr(2.45e-12,-1775.,ztrec) 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 JPL2006 rates_lut(kmo2ho2,k)=zfarr(4.1e-13,750.,ztrec) rates_lut(kmo2no,k)=zfarr(2.8e-12,300.,ztrec) rates_lut(kmo2mo2,k)=zfarr(9.5e-14,390.,ztrec) !JEW: changed to JPL2006 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) ! from IUPAC (Atkinson et al, 2006) rates_lut(kc47a,k)=2.7e-28*zt3rec**7.1 rates_lut(kc47b,k)=1.2e-11*zt3rec**0.9 rates_lut(kc48a,k)=zfarr(4.9e-3,-12100.,ztrec) rates_lut(kc48b,k)=zfarr(5.4e16,-13830.,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 rates_lut(kc61a,k)=1.e-28*zt3rec**4.5 rates_lut(kc61b,k)=8.8e-12*zt3rec**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.04e-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*zt3rec**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(-253./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,500.,ztrec)!at 298 1.01e-12 !JEW: changed to JPL2006 rates_lut(kso2oha,k)=3.3e-31*zt3rec**4.3 rates_lut(kso2ohb,k)= 1.6e-12 ! ! fate of ammonia ! rates_lut(knh3oh,k)= zfarr(1.7e-12,-710.,ztrec)!1.56e-13 at 298k 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 rates_lut(knh2o2,k)= 6.0e-21 rates_lut(knh2o3,k)= zfarr(4.3e-12,-930.,ztrec)!1.89e-13 at 298k ! ! 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.5e-15 rates_lut(kohmvk,k) = zfarr(1.86e-11,175.,ztrec) ! IUPAC rates_lut(kohole,k) = zfarr(8.2e-12,610.,ztrec) ! IUPAC rates_lut(kohmacr,k) = zfarr(2.6e-12,610.,ztrec) ! IUPAC rates_lut(ko3mvk,k) = zfarr(8.5e-16,-1520.,ztrec) ! Poesch et al 2000 rates_lut(ko3ole,k) = 1.0e-17 rates_lut(ko3macr,k) = zfarr(1.4e-15,-2100.,ztrec)! IUPAC rates_lut(kno3ole,k) = zfarr(4.0e-14,-400.,ztrec) ! ! **** 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(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 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 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 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 = 0.18e-4 ! assumed particle radius [cm] ! Parameters to describe subgrid scale mixing real :: zcc ! --- 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 !$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 !------ Fix O3 conc (PLS, June 12, 2012) !PREV ! o3=y(i,j,io3) !NOW o3 = rm(i,j,level,io3) / ye(i,j,iairm) * aird * fscale(io3) !kg ----> #/cm3 !------ ! !**** calculate three body reaction rates ! rx1=rates_lut(KNO2OHA,itemp)*aird rx2=rates_lut(KNO2OHB,itemp) rr(i,j,kno2oh)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**2)) rx1=rates_lut(KOHHNO3C,itemp) rx2=rates_lut(KOHHNO3B,itemp)*aird rr(i,j,kohhno3)=rates_lut(KOHHNO3A,itemp)+rx1*rx2/(rx1+rx2) rx1=rates_lut(KNO2NO3A,itemp)*aird rx2=rates_lut(KNO2NO3B,itemp) rr(i,j,kno2no3)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**2)) rx1=rates_lut(KNO2HO2A,itemp)*aird rx2=rates_lut(KNO2HO2B,itemp) rr(i,j,kno2ho2)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**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)) ! ! 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,kn2o5)=rr(i,j,kno2no3)/rates_lut(KN2O5,itemp) rr(i,j,khno4oh)=rates_lut(KHNO4OH,itemp) rr(i,j,khno4m)=rr(i,j,kno2ho2)/rates_lut(KHNO4M,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) rr(i,j,kh2oh)=rates_lut(kh2oh,itemp)*550.e-9*aird !h2=550ppbv JEW update rr(i,j,kohmper)=rates_lut(KOHMPER,itemp) rr(i,j,kohrooh)=rates_lut(KOHROOH,itemp) rr(i,j,kmo2ho2)=rates_lut(KMO2HO2,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) 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)) 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,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,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,krn222)=2.11e-6 ! s-1 decay time Rn222 to Pb210 ! - BEFORE - 25 Jun 2012 - P. Le Sager ! rr(i,j,ko3po2)=rates_lut(ko3po2,itemp)*o2 ! rr(i,j,ko3po3)=rates_lut(ko3po3,itemp) ! rr(i,j,kxo2xo2n)=6.8e-14 ! rr(i,j,kxo2n)=6.8e-14 ! rjx(i,j,jo2)=(2*rjx(i,j,jo2)*o2)*rr(i,j,ko3po2)/ & ! (rr(i,j,ko3po2) + rr(i,j,ko3po3)*o3) ! ! - NOW - 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 ! - - 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) & !$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) g_n2o5_l_temp = 2.7e-5*exp(1800./temp) kt_liq_n2o5=1./((r_cloud/dg)+(4./(2.4e4*g_n2o5_l_temp))) kt_ice_n2o5=1./((r_ice/dg) +(4./(2.4e4*g_n2o5_i))) rr(i,j,kn2o5l)=(kt_liq_n2o5*sad_cloud+kt_ice_n2o5*sad_ice)*zcc !liquid cloud & ice n2o5 ! ho2... kt_liq_ho2=1./((r_cloud/dg)+(4./(2.4e4*g_ho2_l))) kt_ice_ho2=1./((r_ice/dg) +(4./(2.4e4*g_ho2_i))) ! HO2 uptake on clouds (set to 0 until a better description is available) !rr(i,j,kho2_l) =(kt_liq_ho2*sad_cloud +kt_ice_ho2*sad_ice)*zcc !liquid cloud and ice ho2 rr(i,j,kho2_l) = 0. ! output sad_cld_out(i,j)=sad_cloud sad_ice_out(i,j)=sad_ice ! Aerosol uptake rr(i,j,kn2o5_aer)=0. rr(i,j,kno3_aer) =0. rr(i,j,kho2_aer) =0. air_vol=ye(i,j,iairm)*Avog / (aird * xmair) !air mass [kg] / air density in units 1e-3[kg/m3] -> 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./(2.4e4 * g_n2o5(imode)))) kt_aer_no3 =1./((aer_rad/dg)+(4./(2.4e4 * g_no3 (imode)))) kt_aer_ho2 =1./((aer_rad/dg)+(4./(2.4e4 * 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 ! 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./(2.4e4 * g_n2o5(imode)))) ! kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(2.4e4 * g_no3 (imode)))) ! kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(2.4e4 * 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 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./((aer_rad_ref/dg)+(4./(2.4e4 * g_n2o5_aer))) kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(2.4e4 * g_no3_aer))) kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(2.4e4 * 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) ! ok status = 0 end subroutine calrates !EOC #endif /*CHEM*/ END MODULE CHEM_RATES