!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') 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" ! !################################################################# module chem_rates ! ! use GO, only : gol, goPr, goErr ! ! private ! ! public :: rates, calrates, calchetnew1, calchetnew2 ! ! !contains ! ! ! subroutine rates ! !---------------------------------------------------------------------- ! ! ! !**** rates calculation of look up tables for rate constants and ! ! Henry's law constants ! ! ! ! purpose ! ! ------- ! ! calculation of look up tables rate constants/henry coefficients ! ! ! ! interface ! ! --------- ! ! call rates ! ! ! ! method ! ! ------ ! ! use known array of temperatures (integers) to calculate rate constants ! ! 3 body reactions are explicitly calculated in chemistry ! ! ! ! external ! ! -------- ! ! none ! ! ! ! reference ! ! --------- ! ! none ! ! ! !------------------------------------------------------------------ ! use toolbox, only: zfarr ! use chem_param ! ! implicit none ! ! ! local ! integer :: itemp,k,i ! real :: ztrec,zt3rec,temp ! ! ! start ! do k=1,ntemp ! itemp=k+ntlow ! temp=float(itemp) ! ztrec=1./float(itemp) ! zt3rec=300./float(itemp) ! rates_lut(knoo3,k)=zfarr(2.e-12,-1400.,ztrec) ! rates_lut(kho2no,k)=zfarr(3.7e-12,250.,ztrec) ! rates_lut(kno2oha,k)=2.47e-30*zt3rec**2.97 !!!wp!!! new ravi ! rates_lut(kno2ohb,k)=1.45e-11*zt3rec**2.77 !!!wp!!! new ravi ! 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) ! rates_lut(kno2no3a,k)=2.2e-30*zt3rec**3.9 ! rates_lut(kno2no3b,k)=1.5e-12*zt3rec**0.7 ! rates_lut(kn2o5,k)=zfarr(2.7e-27,11000.,ztrec) ! rates_lut(khno4oh,k)=zfarr(1.3e-12,380.,ztrec) ! rates_lut(kno2ho2a,k)=1.8e-31*zt3rec**3.2 ! rates_lut(kno2ho2b,k)=4.7e-12*zt3rec**1.4 ! rates_lut(khno4m,k)=zfarr(2.1e-27,10900.,ztrec) ! rates_lut(kodm,k)=.2095*zfarr(3.2e-11,70.,ztrec)+ & ! .7808*zfarr(1.8e-11,110.,ztrec) ! rates_lut(kh2ood,k)=2.2e-10 ! rates_lut(ko3ho2,k)=zfarr(1.1e-14,-500.,ztrec) ! rates_lut(ko3oh,k)=zfarr(1.6e-12,-940.,ztrec) ! rates_lut(khpoh,k)=zfarr(2.9e-12,-160.,ztrec) ! rates_lut(kfrmoh,k)=1.e-11 ! rates_lut(kch4oh,k)=zfarr(2.65e-12,-1800.,ztrec) ! rates_lut(kohmper,k)=zfarr(3.8e-12,200.,ztrec) ! rates_lut(kohrooh,k)=3.e-12 ! rates_lut(kho2oh,k)=zfarr(4.8e-11,250.,ztrec) ! rates_lut(kmo2ho2,k)=zfarr(3.8e-13,800.,ztrec) ! rates_lut(kmo2no,k)=zfarr(4.2e-12,180.,ztrec) ! rates_lut(kmo2mo2,k)=zfarr(2.5e-13,190.,ztrec) ! rates_lut(kho2ho2a,k)=zfarr(2.3e-13,600.,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(kc43,k)=zfarr(7.0e-12,250.,ztrec) ! rates_lut(kc44,k)=2.5e-15 ! rates_lut(kc46,k)=zfarr(3.5e-11,-180.,ztrec) ! rates_lut(kc47a,k)=2.7e-28*zt3rec**(-7.1) ! rates_lut(kc47b,k)=1.2e-11*zt3rec**0.9 ! rates_lut(kc48,k)=zfarr(2.0e16,-13500.,ztrec) ! rates_lut(kc49,k)=2.e-12 ! rates_lut(kc50,k)=6.5e-12 ! rates_lut(kc52,k)=8.1e-13 ! rates_lut(kc53,k)=zfarr(1.e15,-8000.,ztrec) ! rates_lut(kc54,k)=1.6e3 ! rates_lut(kc57,k)=zfarr(5.2e-12,504.,ztrec) ! rates_lut(kc58,k)=zfarr(4.33e-15,-1800.,ztrec) ! rates_lut(kc59,k)=7.7e-15 ! rates_lut(kc61a,k)=1.e-28*zt3rec**0.8 ! rates_lut(kc61b,k)=8.8e-12 ! rates_lut(kc62,k)=zfarr(9.14e-15,-2580.,ztrec) ! rates_lut(kc73,k)=1.7e-11 ! rates_lut(kc76,k)=zfarr(2.54e-11,410.,ztrec) ! rates_lut(kc77,k)=zfarr(12.3e-15,-2013.,ztrec) ! rates_lut(kc78,k)=7.8e-13 ! rates_lut(kc79,k)=zfarr(4.2e-12,180.,ztrec) ! rates_lut(kc80,k)=zfarr(1.7e-14,1300.,ztrec) ! rates_lut(kc81,k)=6.8e-13 ! rates_lut(kc82,k)=zfarr(3.5e-13,1000.,ztrec) ! rates_lut(kc83,k)=8.e-11 ! rates_lut(kc84,k)=1.78e-12 ! ! ! sulfur and ammonia gas phase reactions ! ! ! the hynes et al. parameterisation is replaced by chin et al. 1996 ! ! rates_lut(kdmsoha,k)= 9.6e-12*exp(-234./temp) ! rates_lut(kdmsohb,k)=1.7e-42*exp(7810./temp) ! rates_lut(kdmsohc,k)=5.5e-31*exp(7460./temp) ! rates_lut(kdmsno3,k)=zfarr(1.9e-13,500.,ztrec)!at 298 1.01e-12 ! rates_lut(kso2oha,k)=3.0e-31*(temp/300.)**(-3.3) ! rates_lut(kso2ohb,k)= 1.5e-12*(temp/300.) ! ! ! ! ! fate of ammonia ! ! ! rates_lut(knh3oh,k)= zfarr(1.7e-12,-710.,ztrec)!1.56e-13 at 298k ! rates_lut(knh2no,k)= zfarr(3.8e-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 ! ! ! ! ! **** solubility Henry equilibrium ! ! HNO3/so4/nh4 just a very high number to take H and ! ! dissociation into account ! ! ! henry(:,k)=0. ! henry(ih2o2,k)=zfarr(1.67e-5,6621.,ZTREC) ! henry(ihno3,k)=1e7 ! henry(ich3o2h,k)=zfarr(1.5e-6,5607.,ZTREC) ! henry(ich2o,k)=zfarr(2.7e-6,6425.,ZTREC) ! henry(irooh,k)=zfarr(1.5e-6,5607.,ZTREC)!(as CH3O2H) ! 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)*3.41e-5 !correction for the 1/298. part ! henry(inh3,k) =75.0*exp(3400.*ZTREC)*1.10e-5 ! henry(io3,k)=1.1e-2*exp(2300.*ZTREC)*4.45e-4 ! end do !k temperature loop ! ! ! ! marked tracers: ! ! ! henry(io3s,:) = henry(io3,:) ! ! ! end subroutine rates ! ! ! ! subroutine calrates(region,rjx,rr,ye) ! !************************************************************** ! ! ! ! CALRATES calculate rate constants using lookup table rates_lut ! ! calculate third bodies ! ! calculate heterogeneous removal on aerosols ! ! ! ! External: CALCHET ! ! ! !************************************************************ ! !debug use hdf ! ! use global_data, only: region_dat ! use chem_param ! use dims, only: isr,ier, jsr,jer, im, jm ! ! implicit none ! ! ! input/output ! integer, intent(in) :: region ! real,dimension(im(region),jm(region),nj) :: rjx ! ! output ! ! rr: reaction rates... ! real,dimension(im(region),jm(region),nreac),intent(out) :: rr ! ! ye: extra 2D fields ! real,dimension(im(region),jm(region),n_extra),intent(inout) :: ye ! ! ! local: ! ! heterogeneous removal rates ! real,dimension(:,:),allocatable :: het_nh3, het_n2o5 ! ! ! help variables ! integer, dimension(:,:), pointer :: zoomed ! integer :: itemp,i,j ! real :: tr, temp, wv, airp, rx1, rx2, rx3 ! real :: dum, h2ox, aird, o2 ! real :: x1, x2, xice, xliq ! ! ! ! cloud chemistry of n2o5 ! real :: dg, kt_liq, kt_ice ! real, parameter :: r_liq=1e-3, r_ice=5e-3 ! cm ! real, parameter :: g_n2o5_i=0.01, g_n2o5_l=0.04 ! ! !debug integer :: io,sfstart,sfend ! ! ! start ! zoomed => region_dat(region)%zoomed ! allocate(het_nh3(im(region),jm(region))) ! allocate(het_n2o5(im(region),jm(region))) ! ! rx3=0.6 ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! 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.24e16) ! ye(i,j,irh) = max(min(ye(i,j,irh),100.0),0.0) !limit rh between 0-100% ! ! o2=0.209476*aird ! ! ! !**** 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(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) ! 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)=1.06*rates_lut(KCH4OH,itemp)*550.e-9*aird !H2=550ppbv ! 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) ! rr(i,j,kc43)=rates_lut(KC43,itemp) ! rr(i,j,kc44)=rates_lut(KC44,itemp) ! rr(i,j,kc46)=rates_lut(KC46,itemp) ! rx1=rates_lut(KC47A,itemp)*aird*0.7808 ! rx2=rates_lut(KC47B,itemp) ! rr(i,j,kc47)=0.96*RX1/(1.+RX1/RX2)*0.3**(1./(1.+LOG10(RX1/RX2)**2)) ! rr(i,j,kc48)=rates_lut(KC48,itemp) ! 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) ! rr(i,j,kc57)=rates_lut(KC57,itemp) ! rr(i,j,kc58)=rates_lut(KC58,itemp) ! rr(i,j,kc59)=rates_lut(KC59,itemp) ! 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) ! 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)=rr(i,j,kc81)*rr(i,j,kc82)/rr(i,j,kc79) ! 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 ! ! 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 ! end do ! end do !_lat/lon loop ! ! ! calculate heterogeneous removal constants of n2o5 ! ! call calchetnew2(region,ye,het_n2o5,1) !n2o5 ! call calchetnew2(region,ye,het_nh3,2) !nh3 ! ! ! ! ! heterogeneous reaction of N2O5 and H2O -> 2 HNO3 on cloud and aerosol ! ! included in gas phase chemistry ! ! ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! ! ! ! kt= (r2/3Dg + 4*r/3vgamma)^-1 ! ! ice r=50 micrometer gamma = 0.01 ! ! water r=10 micrometer gamma = 0.05 ! ! v is 4e5 cm/s and Dg is 0.1 cm2/s at standard press ! airp=ye(i,j,i_pres) ! dg=0.1*1e5/airp !simple approximation for diffusion coeff. [cm2/s] ! kt_liq=1./(r_liq*r_liq/3./dg+4.*r_liq/3./4e5/g_n2o5_l) ! kt_ice=1./(r_ice*r_ice/3./dg+4.*r_ice/3./4e5/g_n2o5_i) ! aird = ye(i,j,iairn) ! xliq = ye(i,j,ilwc) ! xice = ye(i,j,iiwc) ! rr(i,j,kn2o5l)=(kt_liq*xliq+kt_ice*xice) ! cloud ! ! ! ! kn2o5aq and nh3so4 can be done implicitly, ! ! it has occurred that these rates have ! ! become negative over antarctica (aug 1993), ! ! therefore put minimum value of 0. (AJ jul1999) ! ! ! !cmk rr(i,j,kn2o5aq)=max(0.,het_n2o5(i,j)/ ! ! 1e-9*(y(jl,iso4)+y(jl,imsa))/aird ! !cmk multiplication moved to EBI ! rr(i,j,kn2o5aq)=max(0.,het_n2o5(i,j))/1e-9/aird ! ! ! ! knh3so4 is uptake coefficient on H2SO4. ! ! 1 uptake of NH3 consumes 1 acid molecule. ! ! ! rr(i,j,knh3so4)=max(0.,het_nh3(i,j))/1e-9/aird ! ! rr(i,j,knh3so4)=0.0 !CMK sep2003: why was this here? ! ! ! end do ! end do ! ! deallocate(het_n2o5) ! deallocate(het_nh3) ! nullify(zoomed) ! !debug if(level==1) then ! !debug io = sfstart('rr.hdf',dfacc_create) ! !debug call io_write(io,im(region),'im',jm(region),'jm',lm(region),'lm',t,'t') ! !debug call io_write(io,im(region),'im',jm(region),'jm',lm(region),'lm',q,'q') ! !debug call io_write(io,im(region),'im',jm(region),'jm',lm(region),'lm',clwc,'lwc') ! !debug call io_write(io,im(region),'im',jm(region),'jm',lm(region),'lm',ciwc,'iwc') ! !debug call io_write(io,im(region),'im',jm(region),'jm',nrat,'nrat',rr,'rr') ! !debug call io_write(io,nrat,'nrat',ntemp,'ntemp',rates_lut,'rates_lut') ! !debug call io_write(io,ntrace,'ntrace',ntemp,'ntemp',henry,'henry') ! !debug io = sfend(io) ! !debug endif ! end subroutine calrates ! ! ! ! subroutine calchet1(gamma,xmw,icomp) ! !---------------------------------------------------------------------- ! ! ! !**** calchet1 - pre-calculate heterogeneous removal rates on sulfate aerosol ! ! ! ! programmed by frank dentener 01.04.96 ! ! modified by Maarten krol ! ! ! ! purpose ! ! ------- ! ! calculate heterogeneous removal constants for specified species ! ! on sulfate aerosol as function of concentration, ! ! relative humidity and pressure ! ! ! ! interface ! ! --------- ! ! calchet1(gamma,xmw,icomp) ! ! ! ! gamma dimensionless accomodation coefficient ! ! xmw molar weight [g/mol] ! ! icomp (compound number) ! ! ! ! method ! ! ------ ! ! use Whitby sulfate distribution, and Fuchs' rate expression ! ! to integrate rate coefficient on aerosol distribution ! ! ! ! external ! ! -------- ! ! none ! ! ! ! reference ! ! --------- ! ! Dentener (1993) Ph.D. thesis ! ! ! !------------------------------------------------------------------ ! use binas, xgamma=> gamma ! use chem_param ! ! implicit none ! ! ! input ! real,intent(in) :: xmw,gamma ! integer,intent(in):: icomp ! ! ! local ! integer :: ip,isat,i ! real :: press,temp,dxm,dn2o5,vsp,xl,aird,aervol ! real :: hsat1,hsat2,raer,rx1,zlogs,rx2 ! real :: FN1,FN2,FR1,FR2,FA1,FA2,FV1,FV2,rmean,qi30,xkn,xlab,xfac ! real,parameter :: RG=8.314E3,VENT=1.0, FLN10=2.302585,W2PI=2.506638 ! real,parameter :: xmnso4=xmso4+xmh+xmnh4,p1=1.,t1=288.,g=1.40,conc=1e-9 ! real,dimension(3),parameter :: apar=(/1.,3.4e-8,0.301/) ! ! quantities of integration (e.g.number surface, volume and rate coefficient ! integer,parameter :: nt=4 ! integer,parameter :: nint1=2000 ! number of integration intervals ! ! rint: integration stepsizes [m],0-1 um,1-100 um ! real,dimension(nint1) :: rint = & ! (/ (.001E-6, i=1,1000), (0.1E-6, i=1,1000) /) ! real,dimension(nt) :: qi ! ! ! ! ! 1 particle (unity)/cm3, radius and log(sigma) measurements from Whitby ! ! the radius is assumed to be 'dry radius ! ! We take aerosol size from Whitby accumulation mode (1978) ! ! Numbermean radius: 0.034um, sigma=2, 1 (unity) particles cm-3 ! ! Molecular weight NH4HSO4 111 g/mol ! ! aerosol density of dry NH4HSO4 1.8 E3 kg/m3= 1.8 gcm-3 ! ! temperature is not a determining factor is implicitly accounted for ! ! as a function of pressure. ! ! temperature is assumed to follow an adiabatic lapse rate: ! ! (T2/T1)=(P2/P1)^{(g-1)/g} function of pressure with g=Cp/Cv=ca. 1.40 ! ! ! ! ! start ! ! print *,'calchet1: initialize heterogeneous rem. rates' ! ! ! pressure from 1000 to 0 mbar ! do ip=1,11 ! press=max(0.001,1.1-ip*0.1) !atmosphere (minimum is 1 hPa) ! temp=max(210.,t1*(press/p1)**((g-1)/g)) ! !this estimate of temp is not very accurate ! DN2O5=4.6e-6*TEMP**1.75/PRESS*1E-4 ! diffusion coefficient for n2o5 [m2/s] ! dxm=dn2o5*xmw/xmn2o5 ! diffusion coefficient for other component ! VSP=SQRT(8.*RG*TEMP/PI/XMW) ! Molecular speed [m/s] ! XL=3.*DXM/VSP ! free molecular path length [m] ! ! aird=press/(rg*temp)*1e2 ! (mole/cm3) ! aervol=conc*aird*xmnso4/aerdens ! (mole/cm3) * (g/mole)/(g/cm3)=>[cm3/cm3] ! ! aervol is the volume of 1 pbbv dry nh4hso4 at temp and press ! do isat=1,11 ! hsat1=-10.+10.*isat ! HSAT1=AMIN1(HSAT1,95.) !max RH of 95 % ! HSAT2=HSAT1*HSAT1 ! RAER=AMAX1(1.,1.0129231-0.0041328044*hsat1+0.00070336143*hsat2& ! -1.4388956e-05*hsat1*hsat2+9.1359802e-08*hsat2*hsat2) ! ! growth of aerosol radius due to deliquescence for NH4HSO4 ! ! ! actual integration ! RX1=0.0 ! QI(:)=0. ! ZLOGS=APAR(3) ! rmean=apar(2)*raer ! do I=1,NINT1-1 ! RX1=RX1+RINT(I) ! RX2=RX1+RINT(I+1) !RX2-RX1 is integration size interval ! XKN= 3.*DXM/RX1/VSP !Knudsen number ! XLAB=(XKN*4./3.+0.71)/(XKN+1.) ! FN1=(LOG10(RX1/rmean))**2 ! FN2=(LOG10(RX2/rmean))**2 ! FN1= EXP(-FN1/2./ZLOGS/ZLOGS)/RX1 !number integration ! FN2= EXP(-FN2/2./ZLOGS/ZLOGS)/RX2 !number integration ! FR1=1./(1.+XKN*(XLAB+(4.*(1.-GAMMA)/3./GAMMA)))*RX1*FN1 ! !reactivity integration ! FR2=1./(1.+XKN*(XLAB+(4.*(1.-GAMMA)/3./GAMMA)))*RX2*FN2 ! FA1= RX1*RX1*FN1 !surface integration ! FA2= RX2*RX2*FN2 !surface integration ! FV1= RX1*RX1*RX1*FN1 !volume integration ! FV2= RX2*RX2*RX2*FN2 !volume integration ! QI(1)=QI(1)+RINT(I)/2.*(FN1+FN2) !EULER INTEGRATION ! QI(2)=QI(2)+RINT(I)/2.*(FA1+FA2) ! QI(3)=QI(3)+RINT(I)/2.*(FV1+FV2) ! QI(4)=QI(4)+RINT(I)/2.*(FR1+FR2) ! end do !I=1,NINT1-1 ! xfac=APAR(1)*1.e6/FLN10/W2PI/zlogs ! constant integration factor ! QI(1)=QI(1)*xfac*1.E-6 ! conversion cm3=>m3 number ! QI(2)=QI(2)*4.*PI*xfac*1.e-2 ! m=>cm surface ! QI(3)=QI(3)*4./3.*PI*xfac !volume ! QI(4)=QI(4)*4.*PI*DXM*xfac*vent !reactivity ! if (isat == 1) qi30=qi(3) ! dry volume ! hetrem(isat,ip,icomp)=aervol/qi30*qi(4)! removal coefficient ! end do !isat ! end do !ip ! ! end subroutine calchet1 ! ! ! ! subroutine calchet2(region,ye,het,icomp) ! !---------------------------------------------------------------------- ! ! ! !**** calchet2 - calculate heterogeneous removal rates on sulfate aerosol ! ! ! ! programmed by frank dentener 01.04.96 ! ! modified for TM5 by maarten krol jan 2002 ! ! ! ! purpose ! ! ------- ! ! calculate heterogeneous removal constants for specified species ! ! on sulfate aerosol as function of concentration, ! ! relative humidity and pressure ! ! ! ! interface ! ! --------- ! ! calchet2(region,ye,het,icomp) ! ! ! ! region to indicate for which zoom region ! ! ye : extra fields (for rh,pressure). ! ! het heterogeneous removal constant [s-1/ppbv] ! ! icomp (1,2) compound (N2O5, NH3) ! ! ! ! method ! ! ------ ! ! use Whitby sulfate distribution, and Fuchs' rate expression ! ! to integrate rate coefficient on aerosol distribution ! ! ! ! external ! ! -------- ! ! none ! ! ! ! reference ! ! --------- ! ! Dentener (1993) Ph.D. thesis ! ! ! use global_data, only: region_dat ! use binas, xgamma=>gamma ! use dims, only: isr,ier, jsr,jer, im, jm ! use chem_param ! ! implicit none ! ! ! input ! integer,intent(in) :: region,icomp ! real,dimension(im(region),jm(region)) :: het ! result ! real,dimension(im(region),jm(region),n_extra) :: ye ! ! ye: extra fields (rh, pressure) ! ! ! local ! integer,dimension(:,:),pointer :: zoomed ! real :: pres,px,hx,hp1,hp2 ! integer :: np1,np2,nh1,nh2,i,j ! ! ! relative humidity should be between 0-100 % and ! ! pressure between 105000 and 0 Pa ! ! actual interpolation.... ! ! hetrem field in in module..... ! ! zoomed => region_dat(region)%zoomed ! ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! pres = ye(i,j,i_pres) ! np1=min(11,1+nint(10.-pres/10000.)) !pressure ! np1=max(1,np1) ! np2=min(11,np1+1) ! nh1=max(1,nint(ye(i,j,irh)/10.+0.5)) !relative humidity ! nh1=min(11,nh1) ! nh2=min(11,nh1+1) ! px=((11-np1)*10000.-pres)/10000. ! hx=(ye(i,j,irh)-(nh1-1)*10.)/10. ! hp1=px*hetrem(nh1,np2,icomp)+(1.-px)*hetrem(nh1,np1,icomp) ! hp2=px*hetrem(nh2,np2,icomp)+(1.-px)*hetrem(nh2,np1,icomp) ! het(i,j)=hx*hp2+(1.-hx)*hp1 ! end do ! end do ! ! nullify(zoomed) ! ! end subroutine calchet2 ! ! ! ! subroutine calchetnew1(gamma,xmw,icomp) ! !---------------------------------------------------------------------- ! ! ! !**** calchetnew1 - pre- calculate heterogeneous removal rates ! ! on sulfate aerosol ! ! ! ! programmed by frank dentener 01.04.96 ! ! modified MK oct 2003: splitted in two. ! ! ! ! purpose ! ! ------- ! ! calculate heterogeneous removal constants for specified species ! ! on sulfate aerosol as function of concentration, ! ! relative humidity and pressure ! ! ! ! interface ! ! --------- ! ! calchetnew1(gamma,xmw,icomp) ! ! ! ! gamma dimensionless accomodation coefficient ! ! xmw molar weight [g/mol] ! ! icomp component number: 1: n2o5 2:nh3 ! ! ! ! method ! ! ------ ! ! use Whitby sulfate distribution, and Fuchs' rate expression ! ! to integrate rate coefficient on aerosol distribution ! ! ! ! external ! ! -------- ! ! none ! ! ! ! reference ! ! --------- ! ! Dentener (1993) Ph.D. thesis ! ! ! !------------------------------------------------------------------ ! use toolbox, only : escape_tm ! use reaction_data, only : ncomponent, hetrem, aerdens, rra ! use reaction_data, only : nr_interval, np_interval ! ! implicit none ! ! ! input/output ! real, intent(in) :: gamma ! real, intent(in) :: xmw ! integer, intent(in) :: icomp ! ! !local ! integer :: ip,i,iaero ! ! quantities of integration e.g. number surface, volume and rate coefficient ! integer,parameter :: nt=4 ! integer,parameter :: nint1=2000 !number of integration intervals ! real :: rx1,rx2,fn1,fn2,fr1,fr2,fa1,fa2,fv1,fv2 ! real :: zlogs,rmean,qi(nt),qi30 ! real,dimension(nint1) :: rint = & ! (/ (.001E-6, i=1,1000), (0.1E-6, i=1,1000) /) ! ! ! ! needed for calculation of reaction constants ! real,parameter :: xmn2o5 = 108. ! real,parameter :: rg = 8.314e3 ! real,parameter :: vent = 1.0 ! real,parameter :: pi = 3.14159 ! real,parameter :: fln10 = 2.302585 ! real,parameter :: w2pi = 2.506638 ! real,parameter :: avo = 6.0e23 ! real,parameter :: xmnso4 = 111. ! real,parameter :: p1 = 1. ! real,parameter :: t1 = 288. ! real,parameter :: g = 1.40 ! real,parameter :: conc = 1e-9 ! ! ! real :: temp, press, dn2o5, vsp, xl, xfac, xkn, xlab, dxm ! real :: raer,aervol,aird ! real,dimension(3) :: apar=(/1.,3.4e-8,0.301/) ! ! ! 1 particle (unity)/cm3, radius and log(sigma) measurements from Whitby ! ! the radius is assumed to be 'dry radius ! ! We take aerosol size from Whitby accumulation mode (1978) ! ! Numbermean radius: 0.034um, sigma=2, 1 (unity) particles cm-3 ! ! Molecular weight NH4HSO4 115 g/mol ! ! aerosol density of dry NH4HSO4 1.8 E3 kg/m3= 1.8 gcm-3 ! ! temperature is not a determining factor is implicitly accounted for ! ! as a function of pressure. ! ! temperature is assumed to follow an adiabatic lapse rate: ! ! (T2/T1)=(P2/P1)^{(g-1)/g} function of pressure with g=Cp/Cv=ca. 1.40 ! ! ! ! if ( icomp > ncomponent ) then ! call escape_tm('calchetnew1: Check component in calchetnew1') ! end if ! print *,'calchetnew1: initialize heterogeneous rem. rates ', icomp ! ! do ip=1,np_interval !pressure from 1000 to 0 mbar ! ! press=max(0.001,1.1-ip*0.1) !atmosphere (minimum is 1 hPa) ! temp=max(210.,t1*(press/p1)**((g-1)/g)) ! !this estimate of temp is not very accurate ! DN2O5=4.6e-6*TEMP**1.75/PRESS*1E-4 ! diffusion coefficient for n2o5 [m2/s] ! dxm=dn2o5*xmw/xmn2o5 ! diffusion coefficient for other component ! VSP=SQRT(8.*RG*TEMP/PI/XMW) ! Molecular speed [m/s] ! XL=3.*DXM/VSP ! free molecular path length [m] ! ! aird=press/(rg*temp)*1e2 ! (mole/cm3) ! aervol=conc*aird*xmnso4/aerdens ! (mole/cm3) * (g/mole)/(g/cm3)=>[cm3/cm3] ! ! aervol is the volume of 1 pbbv dry nh4hso4 at temp and press ! ! do iaero=1,nr_interval ! aerosol increased radius loop ! ! ! RAER: growth of aerosol radius due to deliquescence for NH4HSO4 ! RAER=rra(iaero) ! ! ! actual integration ! RX1=0.0 ! QI(:)=0. ! ZLOGS=APAR(3) ! rmean=apar(2)*raer ! ! do I=1,NINT1-1 ! ! RX1=RX1+RINT(I) ! RX2=RX1+RINT(I+1) !RX2-RX1 is integration size interval ! XKN= 3.*DXM/RX1/VSP !Knudsen number ! XLAB=(XKN*4./3.+0.71)/(XKN+1.) ! FN1=(LOG10(RX1/rmean))**2 ! FN2=(LOG10(RX2/rmean))**2 ! FN1= EXP(-FN1/2./ZLOGS/ZLOGS)/RX1 !number integration ! FN2= EXP(-FN2/2./ZLOGS/ZLOGS)/RX2 !number integration ! FR1=1./(1.+XKN*(XLAB+(4.*(1.-GAMMA)/3./GAMMA)))*RX1*FN1 ! !reactivity integration ! FR2=1./(1.+XKN*(XLAB+(4.*(1.-GAMMA)/3./GAMMA)))*RX2*FN2 ! FA1= RX1*RX1*FN1 !surface integration ! FA2= RX2*RX2*FN2 !surface integration ! FV1= RX1*RX1*RX1*FN1 !volume integration ! FV2= RX2*RX2*RX2*FN2 !volume integration ! QI(1)=QI(1)+RINT(I)/2.*(FN1+FN2) !EULER INTEGRATION ! QI(2)=QI(2)+RINT(I)/2.*(FA1+FA2) ! QI(3)=QI(3)+RINT(I)/2.*(FV1+FV2) ! QI(4)=QI(4)+RINT(I)/2.*(FR1+FR2) ! ! end do !I=1,NINT1-1 ! ! xfac=APAR(1)*1.e6/FLN10/W2PI/zlogs ! constant integration factor ! QI(1)=QI(1)*xfac*1.E-6 ! conversion cm3=>m3 number ! QI(2)=QI(2)*4.*PI*xfac*1.e-2 ! m=>cm surface ! QI(3)=QI(3)*4./3.*PI*xfac !volume ! QI(4)=QI(4)*4.*PI*DXM*xfac*vent !reactivity ! if (iaero == 1) qi30=qi(3) ! dry volume ! hetrem(iaero,ip,icomp)=aervol/qi30*qi(4)! removal coefficient ! ! end do !iaero ! write(*,'(a,i3,19(1x,1pe9.1))') ' calchetnew1: ',ip,hetrem(:,ip,icomp) ! ! end do !ip ! ! end subroutine calchetnew1 ! ! ! ! subroutine calchetnew2(region,ye,het,icomp) ! !---------------------------------------------------------------------- ! ! ! !**** calchetnew2 - calculate heterogeneous removal rates on sulfate aerosol ! ! ! ! programmed by frank dentener 01.04.96 ! ! modified for TM5 by maarten krol jan 2002 ! ! ! ! purpose ! ! ------- ! ! calculate heterogeneous removal constants for specified species ! ! on sulfate aerosol as function of concentration, ! ! relative humidity and pressure ! ! ! ! interface ! ! --------- ! ! calchetnew2(region,ye,het,icomp) ! ! ! ! region to indicate for which zoom region ! ! ye : extra fields (for rh,pressure). ! ! het heterogeneous removal constant [s-1/ppbv] ! ! icomp (1,2) compound (N2O5, NH3) ! ! ! ! ! ! method ! ! ------ ! ! use Whitby sulfate distribution, and Fuchs' rate expression ! ! to integrate rate coefficient on aerosol distribution ! ! 1 particle (unity)/cm3, radius and log(sigma) measurements from Whitby ! ! the radius is assumed to be 'dry radius ! ! We take aerosol size from Whitby accumulation mode (1978) ! ! Numbermean radius: 0.034um, sigma=2, 1 (unity) particles cm-3 ! ! Molecular weight NH4HSO4 115 g/mol ! ! aerosol density of dry NH4HSO4 1.8 E3 kg/m3= 1.8 gcm-3 ! ! temperature is not a determining factor is implicitly accounted for ! ! as a function of pressure. ! ! temperature is assumed to follow an adiabatic lapse rate: ! ! (T2/T1)=(P2/P1)^{(g-1)/g} function of pressure with g=Cp/Cv=ca. 1.40 ! ! ! ! external ! ! -------- ! ! none ! ! ! ! reference ! ! --------- ! ! Dentener (1993) Ph.D. thesis ! !---------------------------------------------------------------------- ! ! use global_data, only : region_dat ! use binas, only : xgamma=>gamma ! use dims, only : isr,ier, jsr,jer, im, jm ! use reaction_data, only : ncomponent, hetrem, nr_interval ! use reaction_data, only : np_interval, aerdens, rra ! use chem_param, only : n_extra, irinc, i_pres ! ! implicit none ! ! ! input ! integer, intent(in) :: region ! integer, intent(in) :: icomp ! real,dimension(im(region),jm(region)) :: het !result (time scale) ! real,dimension(im(region),jm(region),n_extra) :: ye ! ! ye: extra fields ( pressure, rinc) ! ! ! local ! integer,dimension(:,:),pointer :: zoomed ! ! real :: pres,px,hx,hp1,hp2 ! integer :: np1,np2,nh1,nh2,i,j,jr, nr1, nr2 ! ! ! start ! ! zoomed => region_dat(region)%zoomed ! ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! pres = ye(i,j,i_pres) ! np1=min(np_interval,1+nint(10.-pres/10000.)) ! np1=max(1,np1) ! np2=min(np_interval,np1+1) ! nr1=1 ! do jr=1,nr_interval ! if(ye(i,j,irinc).ge.rra(jr)) nr1=jr ! lower bound of rinc array ! end do ! nr2=min(nr1+1,nr_interval) ! upper bound of rinc ! px=((np_interval-np1)*10000.-pres)/10000. ! hx=1. ! if (nr1.ne.nr2) hx=(ye(i,j,irinc)-rra(nr1))/(rra(nr2)-rra(nr1)) ! hp1=px*hetrem(nr1,np2,icomp)+(1.-px)*hetrem(nr1,np1,icomp) ! hp2=px*hetrem(nr2,np2,icomp)+(1.-px)*hetrem(nr2,np1,icomp) ! het(i,j)=hx*hp2+(1.-hx)*hp1 ! end do ! end do ! nullify(zoomed) ! ! end subroutine calchetnew2 ! ! end module chem_rates