#include #ifdef STRAT MODULE HETCHEM ! This module contains the heterogeneous chemistry routines for ! the stratosphere from the REPROBUS CTM provided by F. Lefevre. ! They were adapted to LMDz-INCA by L. Jourdain and reworked by ! by D. Hauglustaine and F. Jegou. DH, LSCE, 06/2005 ! ! Updated version 2018, DH. USE INCA_DIM IMPLICIT NONE INTEGER, PARAMETER :: nlat = 1 INTEGER, PARAMETER :: nivtop = 32 INTEGER, PARAMETER :: nivbot = 15 INTEGER, PARAMETER :: nbcon = PCNST REAL, SAVE, ALLOCATABLE :: parthno3(:,:) REAL, SAVE, ALLOCATABLE :: parthno3s(:,:) REAL, SAVE, ALLOCATABLE :: parthno3l(:,:) REAL, SAVE, ALLOCATABLE :: parthno3sl(:,:) REAL, SAVE, ALLOCATABLE :: parth2o(:,:) REAL, SAVE, ALLOCATABLE :: parth2os(:,:) REAL, SAVE, ALLOCATABLE :: parth2ol(:,:) REAL, SAVE, ALLOCATABLE :: parthcl(:,:) REAL, SAVE, ALLOCATABLE :: parthcll(:,:) REAL, SAVE, ALLOCATABLE :: parthcls(:,:) REAL, SAVE, ALLOCATABLE :: parthbr(:,:) REAL, SAVE, ALLOCATABLE :: parthbrl(:,:) REAL, SAVE, ALLOCATABLE :: parthbrs(:,:) REAL, SAVE, ALLOCATABLE :: vsed(:,:) REAL, SAVE, ALLOCATABLE :: surfarealiq(:,:) REAL, SAVE, ALLOCATABLE :: surfareanat(:,:) REAL, SAVE, ALLOCATABLE :: surfareaice(:,:) REAL, SAVE, ALLOCATABLE :: qh2so4(:,:) !$OMP THREADPRIVATE(parthno3,parthno3s,parthno3l,parthno3sl,parth2o, parth2os,parth2ol) !$OMP THREADPRIVATE(parthcl, parthcll,parthcls, parthbr,parthbrl,parthbrs,vsed,qh2so4) !$OMP THREADPRIVATE(surfarealiq,surfareanat,surfareaice) CONTAINS SUBROUTINE INIT_HETCHEM USE INCA_DIM IMPLICIT NONE ALLOCATE ( vsed (PLON,PLEV)) ALLOCATE ( parthno3 (PLON,PLEV)) ALLOCATE ( parthno3s (PLON,PLEV)) ALLOCATE ( parthno3l (PLON,PLEV)) ALLOCATE ( parthno3sl (PLON,PLEV)) ALLOCATE ( parth2o (PLON,PLEV)) ALLOCATE ( parth2os (PLON,PLEV)) ALLOCATE ( parth2ol (PLON,PLEV)) ALLOCATE ( parthcl (PLON,PLEV)) ALLOCATE ( parthcll (PLON,PLEV)) ALLOCATE ( parthcls (PLON,PLEV)) ALLOCATE ( parthbr (PLON,PLEV)) ALLOCATE ( parthbrl (PLON,PLEV)) ALLOCATE ( parthbrs (PLON,PLEV)) ALLOCATE ( qh2so4 (PLON,PLEV)) ALLOCATE ( surfarealiq (PLON,PLEV)) ALLOCATE ( surfareanat (PLON,PLEV)) ALLOCATE ( surfareaice (PLON,PLEV)) vsed(:,:) = 0. parthno3(:,:) = 1. parthno3s(:,:) = 1. parthno3l(:,:) = 1. parthno3sl(:,:) = 1. parth2o(:,:) = 1. parth2os(:,:) = 1. parth2ol(:,:) = 1. parthcl(:,:) = 1. parthcll(:,:) = 1. parthcls(:,:) = 1. parthbr (:,:) = 1. parthbrl(:,:) = 1. parthbrs(:,:) = 1. qh2so4(:,:) = 0. surfarealiq(:,:) = 0. surfareanat(:,:) = 0. surfareaice(:,:) = 0. END SUBROUTINE INIT_HETCHEM SUBROUTINE ANALYTIC( & temp, & pmid, & hnm, & qj, & k1, & k2, & k3, & k4, & k5, & k6, & k7, & k8, & k9, & h2o, & mair ) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c F. LEFEVRE (LATMOS/CNRS) from REPROBUS !c !c Adapted to INCA: !c !c D. HAUGLUSTAINE (LSCE) !c L. JOURDAIN (LATMOS) !c F. JEGOU (LSCE) !c R. VALORSO (LIVE) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc USE SPECIES_NAMES USE SAD_COM, ONLY : MASS,VOL USE INCA_DIM USE AIRPLANE_SRC, ONLY : ptrop REAL, PARAMETER :: r=8.205e-5 REAL, PARAMETER :: ctoa=7.336e21 REAL, PARAMETER :: sigma=1.6 INTEGER :: i, ilon, iniv, j REAL :: t, ph2o, zh2o, zh2oeq REAL :: pi, amu REAL :: qh2o, qhno3, qhcl REAL :: zmt, zbt, zhno3eq REAL :: pn0t, pn REAL :: hno3s, h2os REAL :: xsb, hsb, xnb, hnb REAL :: tf, tt, tr, pr REAL :: phi REAL :: wsf, wnf, partf REAL :: vliq, rmode REAL :: ph2o0, phcl1 REAL :: sclono2, shocl REAL :: rpart, rhopart REAL :: ynre REAL :: vsedliq, shapek, vsedsol REAL :: wts REAL temp(PLON,PLEV), ptot(PLON,PLEV), hnm(PLON,PLEV) REAL ck(PLON,PLEV,nbcon),qj(PLON,PLEV,nbcon) REAL k1(PLON,PLEV), k2(PLON,PLEV), k3(PLON,PLEV) REAL k4(PLON,PLEV), k5(PLON,PLEV), k6(PLON,PLEV) REAL k7(PLON,PLEV), k8(PLON,PLEV), k9(PLON,PLEV) REAL h2o(PLON,PLEV),pmid(PLON,PLEV) REAL a,b,c,det,msb,mnb,mcl,mnf,msf REAL qn(10),qs(10),kn(7),ks(7) SAVE qn,qs,kn,ks !$OMP THREADPRIVATE(qn,qs,kn,ks) REAL, INTENT(in) :: mair(PLON,PLEV) ! total atm density REAL :: VOL_S(PLON,PLEV) REAL dens(PLON,PLEV) REAL ns(PLON,PLEV) REAL pn0(PLON,PLEV), phcl0(PLON,PLEV) REAL pw(PLON,PLEV), ms(PLON,PLEV), mn(PLON,PLEV) REAL ws(PLON,PLEV), wn(PLON,PLEV) REAL hhcl(PLON,PLEV), hhocl(PLON,PLEV) REAL hhobr(PLON,PLEV), hclono2(PLON,PLEV) REAL hhbr(PLON,PLEV) REAL tice(PLON,PLEV), rice(PLON,PLEV) REAL rnat(PLON,PLEV), aliq(PLON,PLEV) REAL rmean(PLON,PLEV) REAL mh2so4c, bh2so4, ph2so4, muh2so4 INTEGER condliq(PLON,PLEV) INTEGER lnat(PLON,PLEV), lice(PLON,PLEV) !c !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c ajouts jpl 2000 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c REAL aw(PLON,PLEV) REAL wt(PLON,PLEV), mh2so4(PLON,PLEV) REAL m, y1, y2, bigx REAL a1(3), b1(3), c1(3), d1(3) REAL a2(3), b2(3), c2(3), d2(3) DATA a1/12.37208932, 11.820654354, -180.06541028/ DATA b1/-0.16125516114, -0.20786404244, -0.38601102592/ DATA c1/-30.490657554, -4.807306373, -93.317846778/ DATA d1/-2.1133114241, -5.1727540348, 273.88132245/ DATA a2/13.455394705, 12.891938068, -176.95814097/ DATA b2/-0.1921312255,-0.232338847708,-0.36257048154/ DATA c2/-34.285174607, -6.4261237757, -90.469744201/ DATA d2/-1.7620073078, -4.9005471319, 267.45509988/ REAL ndrop(PLON,PLEV) REAL nice(PLON,PLEV), nnat(PLON,PLEV) REAL nucleimin REAL nnat_small_max, vnat_small_max REAL nnats(PLON,PLEV) REAL nnatl(PLON,PLEV) REAL rnats, rnatl REAL anat(PLON,PLEV), vnat(PLON,PLEV), vnatl(PLON,PLEV) REAL aice(PLON,PLEV), vice(PLON,PLEV) REAL vsed(PLON,PLEV) REAL rverif(PLON*PLEV) REAL rhoair, rhoice, rhonat, rholiq REAL netha, lambda, lambda0, nre, us, cdnre2 REAL bnre(0:6), tcel !c !c pruppacher and klett, 1988 !c DATA bnre /-3.18657,0.992696,-0.153193e-2, & -0.987059e-3,-0.578878e-3,0.855176e-4,& -0.327815e-5/ !c !c carslaw et al., 1995 !c DATA qn/14.5734,0.0615994,-1.14895,0.691693, -0.098863,& 0.0051579,0.123472,-0.115574,0.0110113,0.0097914/ DATA qs/14.4700,0.0638795,-3.29597,1.778224,-0.223244, & 0.0086486,0.536695,-0.335164,0.0265153,0.0157550/ DATA kn/-39.136,6358.4,83.29,-17650.0,198.53,-11948,-28.469/ DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/ LOGICAL pscliq, pscnat, pscice !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! choix du schema de psc: ! ! pscnat = true : autoriser les nat ! pscnat + pscice = true : autoriser les nat et la glace ! pscnat + pscice + pscliq = true : autoriser nat, glace et liquide ! !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc pscnat = .true. pscice = .true. pscliq = .true. !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c initialisations !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c nnat_small_max = 1. ! densite maximale de nat (cm-3) rnats = 0.5e-4 ! rayon du small mode de nat (cm) rnatl = 6.5e-4 ! rayon du large mode de nat (cm) rice = 10.e-4 ! rayon de la glace nucleimin = 1.e-3 pi = 2.*ASIN(1.0) amu = 1.66e-24 ptot(:,:)=pmid(:,:)/100. vsed = 0. ws(:,:) = 0. ! Remove values below tropopause for sulfates DO i=1,PLON DO j=1,PLEV IF (pmid(i,j).GT.ptrop(i)) THEN VOL_S(i,j)=0 ELSE VOL_S(i,j)=VOL(i,j) ENDIF ENDDO ENDDO DO ilon = 1, PLON DO iniv = nivbot, nivtop qh2o = h2o(ilon,iniv) pw(ilon,iniv) = qh2o*ptot(ilon,iniv)/1013.0 pw(ilon,iniv) = MIN(pw(ilon,iniv),2.e-3/1013.) qhno3 = qj(ilon,iniv,id_HNO3) pn0(ilon,iniv)= ptot(ilon,iniv)/1013.0*qhno3 qhcl = qj(ilon,iniv,id_HCl) phcl0(ilon,iniv)=ptot(ilon,iniv)/1013.0*qhcl !c !c all h2so4 is assumed to be in the droplets. you might !c want to input directly the moles h2so4/m3, which is a more !c normal unit for h2so4 amounts, but remember that this will !c vary with pressure and temperature. pinatubo can be simulated !c by multiplying qh2so4 by different factors. !c !c DH. Note: so4 is volume_aerosols/volume_air and is based on aerosol surface !c density from Considine provided from SAGE and in um2/cm3. Then we transform !c the volume into nb moles below. !qh2so4 = MAX((MASS(ilon,iniv)/mair(ilon,iniv)),1.e-20) !qh2so4 = MAX(VOL_S(ilon,iniv)*1e-12,1.e-20) !qh2so4 = MAX(so4(ilon,iniv),1.e-20) !qh2so4 = 1.e-20 !c The line below replaced by the calculation below: !c ns(ilon,iniv)=ptot(ilon,iniv)*100.0*qh2so4/8.314/t ndrop(ilon,iniv) = 10. END DO END DO DO ilon = 1, PLON DO iniv = nivbot, nivtop lnat(ilon,iniv) = 0 lice(ilon,iniv) = 0 rmean(ilon,iniv) = 0. aliq(ilon,iniv) = 0. nnats(ilon,iniv) = 0. nnatl(ilon,iniv) = 0. nnat(ilon,iniv) = 0. nice(ilon,iniv) = 0. rnat(ilon,iniv) = 0. vnat(ilon,iniv) = 0. vnatl(ilon,iniv) = 0. anat(ilon,iniv) = 0. vice(ilon,iniv) = 0. aice(ilon,iniv) = 0. vsed(ilon,iniv) = 0. parthno3s(ilon,iniv) = 1. parthno3sl(ilon,iniv) = 1. parthno3l(ilon,iniv) = 1. parthno3(ilon,iniv) = 1. parth2o(ilon,iniv) = 1. parthcl(ilon,iniv) = 1. END DO END DO !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c traitement des aerosols solides !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF (pscnat) THEN !c !cc traitement du nat !c DO ilon = 1, PLON DO iniv = nivbot, nivtop t = MIN(temp(ilon,iniv),273.) pn0t = pn0(ilon,iniv)*1013.*0.75 zmt = -2.7836 - 0.00088*t zbt = 38.9855 - 11397.0/t + 0.009179*t zhno3eq = 10.0**(zmt*LOG10(pw(ilon,iniv)*1013.*0.75) + zbt) IF (pn0t/zhno3eq .GE. 1.) THEN lnat(ilon,iniv) = 1 pn = zhno3eq/(0.75*1013.) hno3s = (pn0(ilon,iniv) - pn)*1013. & *hnm(ilon,iniv)/ptot(ilon,iniv) hno3s = MAX(hno3s,0.) vnat(ilon,iniv) = hno3s*amu*117./1.35 ! volume maximal du small mode vnat_small_max = nnat_small_max*(4./3.*pi*rnats**3.) ! densite de nat (cm-3) pour les deux modes if (vnat(ilon,iniv) .gt. vnat_small_max) then nnats(ilon,iniv) = nnat_small_max nnatl(ilon,iniv) = (vnat(ilon,iniv) - vnat_small_max)/(4./3.*pi*rnatl**3.) vnatl(ilon,iniv) = vnat(ilon,iniv) - vnat_small_max else nnats(ilon,iniv) = vnat(ilon,iniv)/(4./3.*pi*rnats**3.) nnatl(ilon,iniv) = 0. vnatl(ilon,iniv) = 0. end if ! surface du nat anat(ilon,iniv) = nnats(ilon,iniv)*4*pi*rnats**2& + nnatl(ilon,iniv)*4*pi*rnatl**2 ! parthno3s: part en phase gazeuse parthno3s(ilon,iniv) = 1.0 - (pn0(ilon,iniv)-pn)& /pn0(ilon,iniv) parthno3s(ilon,iniv) = max(parthno3s(ilon,iniv),0.0) parthno3s(ilon,iniv) = min(parthno3s(ilon,iniv),1.0) ! parthno3sl: part non sedimentee (gaz + petit mode de nat) parthno3sl(ilon,iniv) = 1.0& - (pn0(ilon,iniv)-pn)*vnatl(ilon,iniv)/vnat(ilon,iniv)/pn0(ilon,iniv) parthno3sl(ilon,iniv) = max(parthno3sl(ilon,iniv),0.0) parthno3sl(ilon,iniv) = min(parthno3sl(ilon,iniv),1.0) ENDIF ENDDO ENDDO ENDIF !pscnat IF ( pscice ) THEN DO ilon = 1, PLON DO iniv = nivbot, nivtop !cc traitement de la glace t = MIN(temp(ilon,iniv),273.) tice(ilon,iniv) = 2668.70/(10.4310-(log(pw(ilon,iniv))+log(760.0))/log(10.0)) IF (t .LE. tice(ilon,iniv) ) THEN lice(ilon,iniv) = 1 lnat(ilon,iniv) = 0 anat(ilon,iniv) = 0. zh2oeq = 10.0**(-2668.7/t + 10.4310) ph2o = zh2oeq/760.0 h2os = MAX(pw(ilon,iniv) - ph2o,0.) & *1013.*hnm(ilon,iniv)/ptot(ilon,iniv) vice(ilon,iniv) = h2os*amu*18./0.928 nice(ilon,iniv) = vice(ilon,iniv)/(4./3.*pi*rice(ilon,iniv)**3.) aice(ilon,iniv) = nice(ilon,iniv)*4*pi*rice(ilon,iniv)**2 ! parth2o : part en phase gazeuse parth2o(ilon,iniv) = (1.0-MAX(pw(ilon,iniv)-ph2o,0.)/pw(ilon,iniv)) parth2o(ilon,iniv) = MAX(parth2o(ilon,iniv),0.) parth2o(ilon,iniv) = MIN(parth2o(ilon,iniv),1.) ! On retire hno3 de la phase gazeuse sous forme de nat zmt = -2.7836 - 0.00088*t zbt = 38.9855 - 11397.0/t + 0.009179*t zhno3eq = 10.0**(zmt*log10(pw(ilon,iniv)& *parth2o(ilon,iniv)*1013.*0.75)+ zbt) pn = zhno3eq/(0.75*1013.) ! parthno3s : part en phase gazeuse parthno3s(ilon,iniv) = 1.0 - (pn0(ilon,iniv)-pn)& /pn0(ilon,iniv) parthno3s(ilon,iniv) = MAX(parthno3s(ilon,iniv),0.) parthno3s(ilon,iniv) = MIN(parthno3s(ilon,iniv),1.) parthno3sl(ilon,iniv) = parthno3s(ilon,iniv) ENDIF END DO END DO ENDIF !pscice !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c traitement des aerosols liquides !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF (pscliq) THEN !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c Calcul du nb de moles par unite de volume !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc DO ilon = 1, PLON DO iniv = nivbot, nivtop t = MAX(temp(ilon,iniv),MAX(tice(ilon,iniv)-3.0,185.)) qh2o = h2o(ilon,iniv)*parth2o(ilon,iniv) pw(ilon,iniv) = qh2o*ptot(ilon,iniv)/1013.0 pw(ilon,iniv) = MIN(pw(ilon,iniv),2.e-3/1013.) pw(ilon,iniv) = MAX(pw(ilon,iniv),2.e-5/1013.) xsb = 1.0/(2.0*(ks(3)+ks(4)/t))*( -ks(1)-ks(2)/t- & ((ks(1)+ks(2)/t)**2-4.0*(ks(3)+ks(4)/t) & *(ks(5)+ks(6)/t+ks(7)*LOG(t)-LOG(pw(ilon,iniv))))& **0.5) msb = 55.51*xsb/(1.0-xsb) ms(ilon,iniv) = msb mn(ilon,iniv) = 0.0 ws(ilon,iniv) = msb*0.098076/(1.0+msb*0.098076) wn(ilon,iniv) = 0.0 ws(ilon,iniv) = MAX(ws(ilon,iniv), 0.005) !dilution minimale 0.005% END DO END DO CALL density(ws,wn,temp,dens) DO ilon = 1, PLON DO iniv = nivbot, nivtop qh2so4(ilon,iniv) = MAX(VOL_S(ilon,iniv)*1e-12,1.e-20) ns(ilon,iniv)=qh2so4(ilon,iniv)*1.e6*ws(ilon,iniv)*dens(ilon,iniv)/98.076 END DO END DO !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc aerosols liquides binaires h2so4/h2o !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc DO ilon = 1, PLON DO iniv = nivbot, nivtop condliq(ilon,iniv) = 1 t = MAX(temp(ilon,iniv),MAX(tice(ilon,iniv)-3.0,185.)) qh2o = h2o(ilon,iniv)*parth2o(ilon,iniv) pw(ilon,iniv) = qh2o*ptot(ilon,iniv)/1013.0 pw(ilon,iniv) = MIN(pw(ilon,iniv),2.e-3/1013.) pw(ilon,iniv) = MAX(pw(ilon,iniv),2.e-5/1013.) qhno3 = qj(ilon,iniv,id_HNO3)*parthno3s(ilon,iniv) pn0(ilon,iniv)= ptot(ilon,iniv)/1013.0*qhno3 a = ks(3) + ks(4)/t b = ks(1) + ks(2)/t c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw(ilon,iniv)) det = b**2 - 4.*a*c if (det .gt. 0.) then xsb = (-b - sqrt(det))/(2.*a) msb = 55.51*xsb/(1.0 - xsb) else msb = 0. condliq(ilon,iniv) = 0 end if ms(ilon,iniv) = msb ws(ilon,iniv) = msb*0.098076/(1.0 + msb*0.098076) mn(ilon,iniv) = 0.0 wn(ilon,iniv) = 0.0 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! pression saturante de h2so4 ! d apres ayers et al., grl, 1980 ! et giauque et al., j. am. chem. soc., 1960 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc wts = max(41., ws(ilon,iniv)*100.) muh2so4 = 4.184*(1.514e4 - 286.*(wts - 40.)& + 1.08*(wts - 40.)**2& - 3941./(wts - 40.)**0.1) ! pression saturante en atmosphere ph2so4 = exp(-10156./t + 16.2590 - muh2so4/(8.314*t)) if (qh2so4(ilon,iniv)*ptot(ilon,iniv)/1013. .le. ph2so4) then ms(ilon,iniv) = 0. ws(ilon,iniv) = 0. condliq(ilon,iniv) = 0 end if !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc aerosols liquides ternaires hno3/h2so4/h2o !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IF ( condliq(ilon,iniv) .eq. 1 .and. t .lt. 215.) THEN tf = t tt = r*tf*ns(ilon,iniv) tr = 1.0e4/tf-43.4782608 pr = LOG(pw(ilon,iniv))+18.4 hsb = qs(1)+qs(2)*tr**2+(qs(3)+qs(4)*tr+qs(5)*tr**2 & +qs(6)*tr**3)*pr + (qs(7)+qs(8)*tr+qs(9)*tr**2) & *pr**2+qs(10)*tr*pr**3 hsb = EXP(hsb) xnb = 1.0/(2.0*(kn(3)+kn(4)/tf))*(-kn(1)-kn(2)/tf- & ((kn(1)+kn(2)/tf)**2-4.0*(kn(3)+kn(4)/tf) & *(kn(5)+kn(6)/tf +kn(7)*LOG(tf) & -LOG(pw(ilon,iniv))))**0.5) mnb = 55.51*xnb/(1.0-xnb) hnb = qn(1)+qn(2)*tr**2+(qn(3)+qn(4)*tr+qn(5)*tr**2 & +qn(6)*tr**3)*pr + (qn(7)+qn(8)*tr+qn(9)*tr**2) & *pr**2+qn(10)*tr*pr**3 hnb = EXP(hnb) a = (tt*hnb*mnb**2 - tt*hsb*mnb*msb - 2.0*mnb**2*msb & + mnb*msb**2 + hnb*mnb*msb*pn0(ilon,iniv) & - hsb*msb**2*pn0(ilon,iniv))/(mnb**2 - mnb*msb) b = msb*(-2.0*tt*hnb*mnb+tt*hsb*msb+mnb*msb & -hnb*msb*pn0(ilon,iniv))/(mnb-msb) c = (tt*hnb*mnb*msb**2)/(mnb - msb) phi = ATAN(SQRT(MAX(4.0*(a**2-3.0*b)**3-(-2.0*a**3+9.0*a*b & -27.0*c)**2,0.))/(-2.0*a**3+9.0*a*b-27.0*c) ) IF (phi .LT. 0.) THEN phi = phi + pi ENDIF msf = -1.0/3.0*(a+2.0*SQRT(a**2-3.0*b)*COS((pi+phi)/3.0)) msf = MAX(msf,0.) mnf = mnb*(1.0-msf/msb) mnf = MAX(mnf,0.) pn = mnf/(hnb*mnf/(mnf+msf)+hsb*msf/(mnf+msf)) wsf = msf*0.098076/(1.0+msf*0.098076+mnf*0.063012) wnf = mnf*0.063012/(1.0+msf*0.098076+mnf*0.063012) ms(ilon,iniv) = msf mn(ilon,iniv) = mnf ws(ilon,iniv) = wsf wn(ilon,iniv) = wnf parthno3l(ilon,iniv) = (1.0-(pn0(ilon,iniv)-pn)/pn0(ilon,iniv)) parthno3l(ilon,iniv) = MIN(parthno3l(ilon,iniv),1.) ENDIF END DO END DO !c dilution minimale de h2so4: 0.5% WHERE ( ws < .005 ) ws = 0.005 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc densite des aerosols ternaires (g/cm3) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc CALL density(ws,wn,temp,dens) DO ilon = 1, PLON DO iniv = nivbot, nivtop !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc total specific liquid aerosol volume (dimensionless) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c IF ( condliq(ilon,iniv) .EQ. 1) THEN vliq = ns(ilon,iniv)*98.076e-6/ws(ilon,iniv)/dens(ilon,iniv) !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc liquid aerosol radius and surface area !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c !c rmean=mean radius of liquid droplets (units = cm) !c this varies with the total liquid volume per unit volume of air !c (called vliq in main program). rmean can be calculated !c from vliq using !c !c rmean=rmode*exp(0.5*(log(sigma))**2), where !c !c rmode=(3.0*vliq/(4.*pi*ndrop)* !c exp(-9.0/2.0*(log(sigma)**2)))**(1./3.) !c !c ndrop=number of droplets per cm3 air !c sigma=width of lognormal (use sigma=1.8). !c rmode=(3.0*vliq/(4.*pi*ndrop(ilon,iniv))*EXP(-9.0/2.0*(LOG(sigma)**2)))**(1./3.) rmean(ilon,iniv)=rmode*exp(0.5*(log(sigma))**2) aliq(ilon,iniv)=4.0*pi*rmode**2*ndrop(ilon,iniv) & *EXP(2.0*(LOG(sigma))**2) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc the solubility of hcl !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c !c hhcl est calcule directement a partir de jpl 2000, !c lui-meme utilisant carslaw et al., j. phys. chem., 1995. !c l'ensemble des calculs suppose un aerosol binaire: !c la molalite de h2so4 (wt) est donc recalculee dans !c ces conditions. !c t = MAX(temp(ilon,iniv),MAX(tice(ilon,iniv)-3.0,185.)) !c !c ph2o0 : saturation water vapor pressure, hpa !c ph2o0 = EXP(18.452406985 - 3505.1578807/t & - 330918.55082/t**2 + 12725068.262/t**3) !c !c aw : water activity !c aw(ilon,iniv) = pw(ilon,iniv)*1013./ph2o0 aw(ilon,iniv) = MIN(aw(ilon,iniv), 1.) IF (aw(ilon,iniv) .LE. 0.05) THEN y1 = a1(1)*aw(ilon,iniv)**b1(1) + c1(1)*aw(ilon,iniv) + d1(1) y2 = a2(1)*aw(ilon,iniv)**b2(1) + c2(1)*aw(ilon,iniv) + d2(1) else if (aw(ilon,iniv) .lt. 0.85) then y1 = a1(2)*aw(ilon,iniv)**b1(2) + c1(2)*aw(ilon,iniv) + d1(2) y2 = a2(2)*aw(ilon,iniv)**b2(2) + c2(2)*aw(ilon,iniv) + d2(2) else y1 = a1(3)*aw(ilon,iniv)**b1(3) + c1(3)*aw(ilon,iniv) + d1(3) y2 = a2(3)*aw(ilon,iniv)**b2(3) + c2(3)*aw(ilon,iniv) + d2(3) endif !c !c m: h2so4 molality, mol/kg !c m = y1 + (t - 190.)*(y2 - y1)/70. !c !c wt: h2so4 weight percentage, % !c wt(ilon,iniv) = 9800.*m/(98.*m + 1000.) wt(ilon,iniv) = max(wt(ilon,iniv), 0.5) !c !c mh2so4: h2so4 molarity, mol l-1 !c mh2so4(ilon,iniv) = dens(ilon,iniv)*wt(ilon,iniv)/9.8 bigx = wt(ilon,iniv)/(wt(ilon,iniv) & + (100. - wt(ilon,iniv))*98./18.) hhcl(ilon,iniv) = (0.094 - 0.61*bigx + 1.2*bigx**2) & *exp(-8.68 + (8515. - 10718.*(bigx**0.7))/t) !c !c phcl1 : pression de hcl restant en phase gazeuse !c ns/ms : kilogrammes (litres) d'eau par metre cube !c phcl1 = (phcl0(ilon,iniv)/(r*t)) & /(1./(r*t) + hhcl(ilon,iniv)*ns(ilon,iniv)/ms(ilon,iniv)) parthcl(ilon,iniv) = 1.0-(phcl0(ilon,iniv)-phcl1)/phcl0(ilon,iniv) !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc the solubility of clono2 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c !c on peut negliger la diminution de la concentration !c en phase gazeuse. !c !c sclono2 : setchenow coefficient m-1 !c sclono2 = 0.306 + 24./t !c !c hclono2 : shi et al., m atm-1 !c hclono2(ilon,iniv) = 1.6e-6*exp(4710./t) & *exp(-sclono2*mh2so4(ilon,iniv)) !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc the solubility of hocl !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c !c hhocl d'apres jpl 2000 (shi et al.) !c on peut negliger la diminution de la concentration !c en phase gazeuse. !c !c shocl : setchenow coefficient m-1 !c shocl = 0.0776 + 59.18/t !c !c hhocl : shi et al. !c hhocl(ilon,iniv) = 1.91e-6*exp(5862.4/t) & *exp(-shocl*mh2so4(ilon,iniv)) !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc the solubility of hobr !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c !c hhobr d'apres waschewsky and abbatt, j. phys. chem. a, !c 5312-5320, 1999. !c on peut negliger la diminution de la concentration !c en phase gazeuse. !c hhobr(ilon,iniv) = 4.6e-4*exp(4.5e3/t) !c !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccc the solubility of hbr !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c !c hhbr d'apres kleffmann et al., j. phys. chem. a, !c 8489-8495, 2000. attention a l'erreur typographique !c contenue dans la version imprimee de l'article !c (m3 = - 4.445 imprime au lieu de m3 = + 4.445) !c mh2so4c = -1.977e-4*wt(ilon,iniv)**2. & -2.096e-2*wt(ilon,iniv) + 4.445 bh2so4 = -8.979e-5*wt(ilon,iniv)**2. & +2.141e-2*wt(ilon,iniv) - 6.067 hhbr(ilon,iniv) = 10.**(mh2so4c*1000./t + bh2so4) ENDIF !condliq ENDDO ENDDO ENDIF !pscliq !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c surfaces des PSC !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do ilon = 1,PLON do iniv = nivbot,nivtop surfarealiq(ilon,iniv) = aliq(ilon,iniv) surfareanat(ilon,iniv) = anat(ilon,iniv)*lnat(ilon,iniv) surfareaice(ilon,iniv) = aice(ilon,iniv)*lice(ilon,iniv) end do end do !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c partitions !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do ilon = 1, PLON do iniv = nivbot,nivtop parthno3(ilon,iniv) = parthno3l(ilon,iniv) * parthno3s(ilon,iniv) parthno3(ilon,iniv) = MIN ( parthno3(ilon,iniv) , 1. ) end do end do !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c sedimentation !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc rhoice = 0.928*1000.0 rhonat = 1.350*1000.0 do ilon = 1, PLON do iniv = nivbot,nivtop IF ( lice(ilon,iniv) .EQ. 1 .OR. lnat(ilon,iniv) .EQ. 1 ) THEN t = temp(ilon,iniv) rpart=(lice(ilon,iniv)*rice(ilon,iniv) & +lnat(ilon,iniv)*rnat(ilon,iniv)) *1.e-2 rpart=min(rpart,535.e-6) if (rpart .ge. 5.e-7) then rhopart=lice(ilon,iniv)*rhoice +lnat(ilon,iniv)*rhonat rhoair=(hnm(ilon,iniv)*28.97/6.022e23) *1000.0 tcel = t - 273.15 netha = (1.718 + 0.0049*tcel -1.2e-5*tcel*tcel)*1.e-5 cdnre2 = 32.0*(rpart)**3*(rhopart-rhoair) & *rhoair*9.80665 /(3.0*(netha)**2) !c pruppacher and klett 1998 ynre = bnre(0)+bnre(1)*log(cdnre2) & +bnre(2)*(log(cdnre2))**2 & +bnre(3)*(log(cdnre2))**3 & +bnre(4)*(log(cdnre2))**4 & +bnre(5)*(log(cdnre2))**5 & +bnre(6)*(log(cdnre2))**6 nre = exp(ynre) if (nre .ge. 1.e-2) then !c !c big particles !c (1.e-2 < reynolds number < 300) !c (10.0 < radius < 535 microns) !c vsed(ilon,iniv) = netha*nre & /(2*rhoair*rpart) else !c !c small particules !c (1.e-6 < reynolds number < 1.e-2) !c (0.5 < radius < 10 microns) !c !c lambda0 : libre parcours moyen (m) a 1013.25 hpa !c et 293.15k. voir pruppacher and klett !c page 416. !c lambda0 = 6.6e-8 lambda = lambda0 & *(1013.25/ptot(ilon,iniv)) & *(t/293.15) us = 2.0*rpart**2*9.80665 & *(rhopart-rhoair)/(9.0*netha) !c !ccc solid particles: !c !c niels larsen - 1991 !c columnar or elipsoid shape !c with a shape factor of 1.2 (shapek) !c shapek = 1.2 vsedsol = us/shapek & *(1.0+1.246*lambda/rpart & +0.42*lambda/rpart & *exp(-0.87*rpart/lambda)) vsed(ilon,iniv) = vsedsol endif endif ENDIF end do end do !c !c only the parameters that are needed in 'hetero' are !c included in the variable list. output variables should !c be included in the list as needed! !c call hetero(aw,temp,ptot,hnm,qj,ws, & hhcl,hhocl,hhobr,hclono2,hhbr, & rice,nice,rnat,nnat, & aliq,rmean,lnat,lice,condliq, & k1,k2,k3,k4,k5,k6,k7,k8,k9,wt,h2o) END SUBROUTINE analytic SUBROUTINE HETERO(aw,temp,ptot,hnm,qj,ws, & hhcl,hhocl,hhobr,hclono2,hhbr, & rice,nice,rnat,nnat, & aliq,rmean,lnat,lice,condliq, & k1,k2,k3,k4,k5,k6,k7,k8,k9,wt,h2o) !C ROUTINE TO CALCULATE UPTAKE COEFFICIENTS (GAMMA VALUES). !C GAMMA VALUES ARE INDICATED BY VARIABLES WITH PREFIX 'G', FOR !C EXAMPLE GHOCLHCL IS THE GAMMA VALUE OF HOCL DUE TO REACTION WITH !C HCL IN THE DROPLETS. !C FROM THE GAMMA VALUES, SECOND ORDER RATE CONSTANTS ARE CALCULATED. !C THESE HAVE THE PREFIX 'R' AND HAVE UNITS CM3 MOLECULE-1 S-1. FOR !C EXAMPLE, THE LOSS OF CLNO3 AND HCL DUE TO THE HETEROGENEOUS REACTION !C CLNO3+HCL -> CL2+HNO3 IS D(CLNO3)/DT (UNITS MOLECULE CM-3 S-1) = !C -RCLNO3HCL.[CLNO3].[HCL], WHERE [CLNO3] AND [HCL] ARE THE !C ****TOTAL**** AMOUNTS OF THESE SPECIES IN UNITS MOLECULE CM-3. !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !c F. LEFEVRE (SA/CNRS) from REPROBUS !c !c Adapted to LMDz-INCA: !c !c L. JOURDAIN (SA) 2004 !c D. HAUGLUSTAINE (LSCE) 06/2005 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc use SPECIES_NAMES INTEGER :: i, ilon, iniv REAL, PARAMETER :: PI=3.1415 REAL, PARAMETER :: ZERO = 1.e-20 REAL :: t0, t, dt REAL :: zh2o, zhcl, zhocl, zhobr REAL :: biga, eta, ah REAL :: dclono2, cclono2 REAL :: gbh2o, grxn, gbhcl, fhcl, gs, gps, gbhclp, gb REAL :: gclono2, gclono2hcl, gclono2h2o REAL :: rclono2h2o, rclono2hcl REAL :: chocl, dhocl, fhocl, ghoclhcl, rhoclhcl REAL :: ratio REAL :: dhobr, fhobr, chobr REAL :: ghobrhcl, rhobrhcl REAL :: dhcl, chcl REAL :: cbar REAL :: gxrn REAL :: gn2o5h2o, rn2o5h2o REAL :: gbrno3h2o, rbrno3h2o REAL :: gclno3hclnat, gclno3h2onat REAL :: gn2o5h2onat, gn2o5hclnat, ghoclhclnat, gbrno3h2onat, ghobrhclnat, ghobrhbrnat, ghoclhbrnat REAL :: rclno3h2onat, rclno3hclnat, rn2o5h2onat, rbrno3h2onat, rn2o5hclnat, rhoclhclnat, rhobrhclnat REAL :: gclno3hclice, gclno3h2oice, gn2o5h2oice, gn2o5hclice, ghoclhclice, gbrno3h2oice, ghobrhclice, ghobrhbrice, ghoclhbrice REAL :: rclno3hclice, rclno3h2oice, rn2o5h2oice, rn2o5hclice, rhoclhclice, rbrno3h2oice, rhobrhclice, rhobrhbrice REAL TEMP(PLON,PLEV), PTOT(PLON,PLEV), HNM(PLON,PLEV) REAL CK(PLON,PLEV,NBCON), qj(PLON,PLEV,NBCON) REAL, DIMENSION(PLON,PLEV) :: k1 REAL, DIMENSION(PLON,PLEV) :: k2 REAL, DIMENSION(PLON,PLEV) :: k3 REAL, DIMENSION(PLON,PLEV) :: k4 REAL, DIMENSION(PLON,PLEV) :: k5 REAL, DIMENSION(PLON,PLEV) :: k6 REAL, DIMENSION(PLON,PLEV) :: k7 REAL, DIMENSION(PLON,PLEV) :: k8 REAL, DIMENSION(PLON,PLEV) :: k9 REAL, DIMENSION(PLON,PLEV) :: k1l REAL, DIMENSION(PLON,PLEV) :: k2l REAL, DIMENSION(PLON,PLEV) :: k3l REAL, DIMENSION(PLON,PLEV) :: k4l REAL, DIMENSION(PLON,PLEV) :: k5l REAL, DIMENSION(PLON,PLEV) :: k6l REAL, DIMENSION(PLON,PLEV) :: k7l REAL, DIMENSION(PLON,PLEV) :: k8l REAL, DIMENSION(PLON,PLEV) :: k9l REAL, DIMENSION(PLON,PLEV) :: k1s REAL, DIMENSION(PLON,PLEV) :: k2s REAL, DIMENSION(PLON,PLEV) :: k3s REAL, DIMENSION(PLON,PLEV) :: k4s REAL, DIMENSION(PLON,PLEV) :: k5s REAL, DIMENSION(PLON,PLEV) :: k6s REAL, DIMENSION(PLON,PLEV) :: k7s REAL, DIMENSION(PLON,PLEV) :: k8s REAL, DIMENSION(PLON,PLEV) :: k9s REAL H2O(PLON,PLEV), H2OCONC(PLON,PLEV) REAL WS(PLON,PLEV) REAL HHCL(PLON,PLEV), HHOCL(PLON,PLEV) REAL HHOBR(PLON,PLEV), HCLONO2(PLON,PLEV) REAL HHBR(PLON,PLEV) REAL RICE(PLON,PLEV) REAL RNAT(PLON,PLEV), ALIQ(PLON,PLEV) REAL RMEAN(PLON,PLEV) REAL NNAT(PLON,PLEV), NICE(PLON,PLEV) INTEGER CONDLIQ(PLON,PLEV) INTEGER LNAT(PLON,PLEV), LICE(PLON,PLEV) !c !c additions pour jpl 2000 !c real aw(PLON,PLEV) real wt(PLON,PLEV) real kh, kh2o, khydr real phcl, mhcl, khcl real pclono2, lclono2, fclono2 real phobr, mhobr, k1hobrhcl, k2hobrhcl real phbr, mhbr, lhbr, k1hoclhbr, k2hoclhbr real khoclhcl, lhocl, phocl, mhocl real lhobr, lhcl ! Initialisation k1l(:,:) = 0. k2l(:,:) = 0. k3l(:,:) = 0. k4l(:,:) = 0. k5l(:,:) = 0. k6l(:,:) = 0. k7l(:,:) = 0. k8l(:,:) = 0. k9l(:,:) = 0. k1s(:,:) = 0. k2s(:,:) = 0. k3s(:,:) = 0. k4s(:,:) = 0. k5s(:,:) = 0. k6s(:,:) = 0. k7s(:,:) = 0. k8s(:,:) = 0. k9s(:,:) = 0. h2oconc(:,:) = h2o(:,:)*hnm(:,:) DO i = 1, nbcon ck(:,:,i) = qj(:,:,i)*hnm(:,:) ENDDO !C !C ****************************************************** !C THE FOLLOWING PARAMETERS ARE NEEDED IN THIS ROUTINE!! !C ****************************************************** !C !C RNAT=MEAN RADIUS OF NAT PARTICLES (CM) !C NNAT=NUMBER OF NAT PARTICLES PER CM3 AIR !C RICE=MEAN RADIUS OF ICE PARTICLES (CM) !C NICE=NUMBER OF ICE PARTICLES PER CM3 AIR !C RMEAN=MEAN RADIUS OF LIQUID AEROSOL DISTRIBUTION (CM) !C ALIQ=TOTAL AREA OF LIQUID AEROSOLS (CM2 CM-3 AIR) DO ILON = 1, PLON DO INIV = nivbot, nivtop IF ( condliq(ilon,iniv) .EQ. 1 ) THEN T = TEMP(ILON,INIV) !C !C Z VARIABLES ARE SET TO AVOID DIVISION BY ZERO !C ZH2O = MAX(h2oconc(ILON,INIV),1.E-08*HNM(ILON,INIV)) ZHCL = MAX(CK(ILON,INIV,id_HCl),1.E-13*HNM(ILON,INIV)) ZHOCL = MAX(CK(ILON,INIV,id_HOCl),1.E-13*HNM(ILON,INIV)) ZHOBR = MAX(CK(ILON,INIV,id_HOBr),1.E-13*HNM(ILON,INIV)) !c !c pressions de hcl, clono2, hobr, et hbr (atm) !c phcl = (ck(ilon,iniv,id_HCl)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. pclono2 = (ck(ilon,iniv,id_ClONO2)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. phocl = (ck(ilon,iniv,id_HOCl)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. phobr = (ck(ilon,iniv,id_HOBr)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. phbr = (ck(ilon,iniv,id_HBr)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. !c !C ===================================================================== !C ====================== CLONO2+HCL -> CL2+HNO3 ======================= !C ====================== CLONO2+H2O -> HOCL+HNO3 ====================== !C ========================= ON LIQUID AEROSOL ========================= !C ===================================================================== !c !c d'apres jpl 2000. On ne tient pas compte de la deposition !c de hno3 (solutions ternaires) sur les probabilites de reaction. !c t0 = 144.11 + 0.166*wt(ilon,iniv) - 0.015*wt(ilon,iniv)**2 & + 2.18e-4*wt(ilon,iniv)**3 biga = 169.5 + 5.18*wt(ilon,iniv) - 0.0825*wt(ilon,iniv)**2 & + 3.27e-3*wt(ilon,iniv)**3 !c !c eta : viscosity of h2so4 solution, cp !c ! Temperature check dt = t - t0 IF ( dt < 1.) THEN dt = 1. PRINT *, 'Warning : temperature problem in hetchem (ILON, T, T0) : ', ilon, t, t0 END IF eta = biga*t**(-1.43)*exp(448./dt) !c !c ah : acid activity in molarity !c ah = exp(60.51 - 0.095*wt(ilon,iniv) + 0.0077*wt(ilon,iniv)**2 & -1.61e-5*wt(ilon,iniv)**3 - (1.76 + & 2.52e-4*wt(ilon,iniv)**2)*sqrt(t) & + (-805.89 + 253.05*wt(ilon,iniv)**0.076)/sqrt(t)) !c !c dclono2 : clono2 diffusivity, klassen et al., cm2 cm-1 !c dclono2 = 5.e-8*t/eta !c !c cbar !c cclono2 = 1474.*sqrt(t) !c !c kh2o : shi et al., s-1 !c kh2o = 1.95e10*exp(-2800./t) !c !c kh : shi et al., s-1 !c kh = 1.22e12*exp(-6200./t) !c !c khydr : shi et al., s-1 !c khydr = kh2o*aw(ilon,iniv) + kh*ah*aw(ilon,iniv) gbh2o = 4*hclono2(ilon,iniv)*0.082*t*sqrt(dclono2*khydr) & / cclono2 mhcl = hhcl(ilon,iniv)*phcl !c !c khcl : shi et al., s-1 !c khcl = 7.9e11*ah*dclono2*mhcl !c !c lclono2 : reacto-diffusive length, cm !c lclono2 = sqrt(dclono2/(khydr + khcl)) !c !c fclono2 !c fclono2 = 1./tanh(rmean(ilon,iniv)/lclono2) & - lclono2/rmean(ilon,iniv) grxn = fclono2*gbh2o*sqrt(1. + khcl/khydr) gbhcl = grxn*khcl/(khcl + khydr) gs = 66.12*hclono2(ilon,iniv)*mhcl*exp(-1374./t) fhcl = 1. & /(1 + 0.612*(gs + gbhcl)*pclono2/phcl) gps = fhcl*gs gbhclp = fhcl*gbhcl gb = gbhclp + grxn*khydr/(khcl + khydr) gclono2 = 1./(1. + 1./(gps + gb)) gclono2hcl = gclono2*(gps + gbhclp)/(gps + gb) gclono2h2o = gclono2 - gclono2hcl rclono2h2o = 0.25*gclono2h2o*cclono2*aliq(ilon,iniv)/zh2o k1l(ilon,iniv) = rclono2h2o rclono2hcl = 0.25*gclono2hcl*cclono2*aliq(ilon,iniv)/zhcl k2l(ilon,iniv) = rclono2hcl !c !C ===================================================================== !C ========================== HOCL + HCL =============================== !C ======================= ON LIQUID AEROSOL =========================== !C ===================================================================== !c !c d'apres jpl 2000 !c !c cbar !c chocl = 2009.*SQRT(t) !c !c dhocl : hocl diffusivity, klassen et al., cm2 cm-1 !c dhocl = 6.4e-8*t/eta !c !c khoclhcl : shi et al., s-1 !c khoclhcl = 1.25e9*ah*dhocl*mhcl !c !c lhocl : reacto-diffusive length, cm !c lhocl = SQRT(dhocl/khoclhcl) !c !c fhocl !c fhocl = 1./TANH(rmean(ilon,iniv)/lhocl) & - lhocl/rmean(ilon,iniv) grxn = 4.*hhocl(ilon,iniv)*0.082*t*SQRT(dhocl*khoclhcl)/chocl !c !c modif line ghoclhcl commenter !c ghoclhcl = 1./(1. + 1./(fhocl*grxn*fhcl)) !cc a la place IF ( (fhocl.LE.zero) .OR. (fhcl.LE.zero) .OR. (grxn.LE.zero) ) THEN ghoclhcl = 0. ELSE ghoclhcl = 1./(1. + 1./(fhocl*grxn*fhcl)) ENDIF rhoclhcl = 0.25*ghoclhcl*chocl*aliq(ilon,iniv)/zhcl k5l(ilon,iniv) = rhoclhcl !C !C ==================================================================== !C =========================== HOBR + HCL ============================= !C ======================= ON LIQUID AEROSOL ========================== !C ==================================================================== !c !c d'apres waschewsky and abbatt, j. phys. chem. a, !c 5312-5320, 1999. !c !c mhobr : concentration de hobr en phase liquide !c mhobr = hhobr(ilon,iniv)*phobr k2hobrhcl = EXP(0.542*wt(ilon,iniv) - 6.44e3/t + 10.3) ratio = mhobr/mhcl IF (ratio .LT. 1.) THEN !c !c hcl est en exces dans l'aerosol. !c hobr est le reactif limitant. !c k1hobrhcl = k2hobrhcl*mhcl !c !c dhobr : hobr diffusivity, klassen et al., cm2 cm-1 !c dhobr = 6.2e-8*t/eta !c !c lhobr : reacto-diffusive length, cm !c lhobr = SQRT(dhobr/k1hobrhcl) fhobr = 1./TANH(rmean(ilon,iniv)/lhobr) & - lhobr/rmean(ilon,iniv) chobr = 1477.*SQRT(t) ghobrhcl = 4.*hhobr(ilon,iniv)*0.082*t & *SQRT(dhobr*k1hobrhcl)*fhobr/chobr rhobrhcl = 0.25*ghobrhcl*chobr*aliq(ilon,iniv)/zhcl ELSE !c !c hobr est en exces dans l'aerosol. !c hcl est le reactif limitant. !c k1hobrhcl = k2hobrhcl*mhobr !c !c dhcl : hcl diffusivity, klassen et al., cm2 cm-1 !c dhcl = 7.8e-8*t/eta !c !c lhcl : reacto-diffusive length, cm !c lhcl = SQRT(dhcl/k1hobrhcl) fhcl = 1./TANH(rmean(ilon,iniv)/lhcl) & - lhcl/rmean(ilon,iniv) chcl = 2408.*SQRT(t) ghobrhcl = 4.*hhcl(ilon,iniv)*0.082*t & *SQRT(dhcl*k1hobrhcl)*fhcl/chcl rhobrhcl = 0.25*ghobrhcl*chcl*aliq(ilon,iniv)/zhobr ENDIF k7l(ilon,iniv) = rhobrhcl !C !C ==================================================================== !C =========================== HOBR + HBR ============================= !C ======================= ON LIQUID AEROSOL ========================== !C ==================================================================== !c k8l(ilon,iniv) = 0. !C !C ==================================================================== !C =========================== HOCL + HBR ============================= !C ======================= ON LIQUID AEROSOL ========================== !C ==================================================================== !c !c mhbr : concentration de hbr en phase liquide !c mhocl : concentration de hocl en phase liquide !c !c mhbr = hhbr(ilon,iniv)*phbr !c mhocl = hhocl(ilon,iniv)*phocl !c !c k2hoclhbr: abbatt and nowak, j. phys. chem. a, 2131, 1997. !c valeur a 228k et 69.3% h2so4 !c !c k2hoclhbr = 2.e6 !c !c ratio = mhbr/mhocl !c !c if (ratio .lt. 1.) then !c !c hocl est en exces dans l'aerosol. !c hbr est le reactif limitant. !c !c k1hoclhbr = k2hoclhbr*mhocl !c !c dhbr : hbr diffusivity, klassen et al., cm2 cm-1 !c !c dhbr = 7.9e-8*t/eta !c chbr = 1616.*sqrt(t) !c !c lhbr : reacto-diffusive length, cm !c !c lhbr = sqrt(dhbr/k1hoclhbr) !c !c fhbr = 1./tanh(rmean(ilon,iniv)/lhbr) !c $ - lhbr/rmean(ilon,iniv) !c !c ghoclhbr = 4.*hhbr(ilon,iniv)*0.082*t !c $ *sqrt(dhbr*k1hoclhbr)*fhbr/chbr !c !c else !c !c hbr est en exces dans l'aerosol. !c hocl est le reactif limitant. !c !c k1hoclhbr = k2hoclhbr*mhbr !c !c lhocl : reacto-diffusive length, cm !c !c lhocl = sqrt(dhocl/k1hoclhbr) !c !c fhocl = 1./tanh(rmean(ilon,iniv)/lhocl) !c $ - lhocl/rmean(ilon,iniv) !c !c ghoclhbr = 4.*hhocl(ilon,iniv)*0.082*t !c $ *sqrt(dhocl*k1hoclhbr)*fhocl/chocl !c endif !c k9l(ilon,iniv) = 0. !C !C ===================================================================== !C ======================== N2O5 + H2O ================================= !C ====================== ON LIQUID AEROSOL ============================ !C ===================================================================== !C GN2O5H2O=0.1 CBAR=1400.1*SQRT(T) RN2O5H2O=0.25*GN2O5H2O*CBAR*ALIQ(ILON,INIV)/ZH2O K3L(ILON,INIV) = RN2O5H2O !c !C ===================================================================== !C ======================== N2O5 + HCl ================================= !C ====================== ON LIQUID AEROSOL ============================ !C ===================================================================== !C K4L(ILON,INIV) = 0. !C !C ===================================================================== !C ======================== BrONO2 + H2O =============================== !C ====================== ON LIQUID AEROSOL ============================ !C ===================================================================== !C !c parametrisation jpl 2000 (hanson, private comm.) !c gxrn = exp(29.24 - 0.396*100.*ws(ilon,iniv)) + 0.114 gbrno3h2o = 1./(1./0.805 + 1./gxrn) CBAR=1221.4*SQRT(T) RBRNO3H2O=0.25*GBRNO3H2O*CBAR*ALIQ(ILON,INIV)/ZH2O K6L(ILON,INIV) = RBRNO3H2O ELSE !condliq K1L(ILON,INIV) = 0. K2L(ILON,INIV) = 0. K3L(ILON,INIV) = 0. K4L(ILON,INIV) = 0. K5L(ILON,INIV) = 0. K6L(ILON,INIV) = 0. K7L(ILON,INIV) = 0. K8L(ILON,INIV) = 0. K9L(ILON,INIV) = 0. ENDIF !condliq END DO END DO DO ILON = 1, PLON DO iniv = nivbot, nivtop IF (LNAT(ILON,INIV) .EQ. 1) THEN T = TEMP(ILON,INIV) !C !C Z VARIABLES ARE SET TO AVOID DIVISION BY ZERO !C ZH2O = MAX(h2oconc(ILON,INIV),1.E-08*HNM(ILON,INIV)) ZHCL = MAX(CK(ILON,INIV,id_HCl),1.E-13*HNM(ILON,INIV)) ZHOCL = MAX(CK(ILON,INIV,id_HOCl),1.E-13*HNM(ILON,INIV)) ZHOBR = MAX(CK(ILON,INIV,id_HOBr),1.E-13*HNM(ILON,INIV)) !C !C ===================================================================== !C ====================== CLNO3 + HCL ================================== !C ====================== CLNO3 + H2O ================================== !C ====================== N2O5 + H2O ================================== !C ====================== N2O5 + HCL ================================== !C ====================== HOCL + HCL ================================== !C ====================== BRNO3 + H2O ================================= !C ====================== HOBR + HCL ================================= !C ====================== HOBR + HBR ================================= !C ====================== HOCL + HBR ================================= !C ======================= ON NAT ====================================== !C ===================================================================== !C !C GAMMA VALUES !C !DH 2018 !GCLNO3HCLNAT=0.1 GCLNO3HCLNAT=0.2 GCLNO3H2ONAT=0.004 GN2O5H2ONAT=4.0E-4 GN2O5HCLNAT=3.0E-3 GHOCLHCLNAT=0.1 GBRNO3H2ONAT=0. GHOBRHCLNAT=0. GHOBRHBRNAT=0. GHOCLHBRNAT=0. !C !C BECAUSE NAT PARTICLES CAN BE VERY LARGE, GAS DIFFUSION LIMITATION !C IS TAKEN INTO ACCOUNT. THE SECOND ORDER RATE CONSTANTS ARE !C CALCULATED FROM THE EQUATION R=G*PI*R**2*CBAR*N/(1+3G.R/(4.L)), WHERE !C G=GAMMA, R=PARTICLE RADIUS, CBAR=MEAN MOLECULAR SPEED, L=MEAN FREE !C PATH (TURCO ET AL, JGR, 94, 16493, 1989). !C RCLNO3H2ONAT = & GCLNO3H2ONAT*4.56E4*SQRT(T/97.45)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GCLNO3H2ONAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZH2O K1S(ILON,INIV) = RCLNO3H2ONAT RCLNO3HCLNAT= & GCLNO3HCLNAT*4.56E4*SQRT(T/97.45)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GCLNO3HCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K2S(ILON,INIV) = RCLNO3HCLNAT RN2O5H2ONAT = & GN2O5H2ONAT*4.56E4*SQRT(T/108.0)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GN2O5H2ONAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZH2O K3S(ILON,INIV) = RN2O5H2ONAT RBRNO3H2ONAT = & GBRNO3H2ONAT*4.56E4*SQRT(T/142.0)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GBRNO3H2ONAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZH2O K6S(ILON,INIV) = RBRNO3H2ONAT RN2O5HCLNAT= & GN2O5HCLNAT*4.56E4*SQRT(T/108.0)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GN2O5HCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K4S(ILON,INIV) = RN2O5HCLNAT RHOCLHCLNAT= & GHOCLHCLNAT*4.56E4*SQRT(T/52.5)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GHOCLHCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K5S(ILON,INIV) = RHOCLHCLNAT RHOBRHCLNAT= & GHOBRHCLNAT*4.56E4*SQRT(T/96.9)*RNAT(ILON,INIV)**2 & *NNAT(ILON,INIV)/ & (1.0+3.3E4*GHOBRHCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K7S(ILON,INIV) = RHOBRHCLNAT K8S(ILON,INIV) = 0. K9S(ILON,INIV) = 0. ENDIF END DO END DO DO ILON = 1, PLON DO iniv = nivbot, nivtop IF (LICE(ILON,INIV) .EQ. 1) THEN T = TEMP(ILON,INIV) !C !C Z VARIABLES ARE SET TO AVOID DIVISION BY ZERO !C ZH2O = MAX(h2oconc(ILON,INIV),1.E-08*HNM(ILON,INIV)) ZHCL = MAX(CK(ILON,INIV,id_HCl),1.E-13*HNM(ILON,INIV)) ZHOCL = MAX(CK(ILON,INIV,id_HOCl),1.E-13*HNM(ILON,INIV)) ZHOBR = MAX(CK(ILON,INIV,id_HOBr),1.E-13*HNM(ILON,INIV)) !C !C ===================================================================== !C ====================== CLNO3 + HCL ================================== !C ====================== CLNO3 + H2O ================================== !C ====================== N2O5 + H2O ================================== !C ====================== N2O5 + HCL ================================== !C ====================== HOCL + HCL ================================== !C ====================== BRNO3 + H2O ================================= !C ====================== HOBR + HCL ================================= !C ====================== HOBR + HBR ================================= !C ====================== HOCL + HBR ================================= !C ======================= ON ICE ====================================== !C ===================================================================== !C GAMMA VALUES !C GCLNO3HCLICE=0.3 GCLNO3H2OICE=0.3 GN2O5H2OICE=0.02 GN2O5HCLICE=0.03 GHOCLHCLICE=0.2 GBRNO3H2OICE=0.3 GHOBRHCLICE=0.3 GHOBRHBRICE=0.1 GHOCLHBRICE=0.1 !FJegou 13/10/06 JPL 2006 !C !C BECAUSE ICE PARTICLES CAN BE VERY LARGE, GAS DIFFUSION LIMITATION !C IS TAKEN INTO ACCOUNT. THE SECOND ORDER RATE CONSTANTS ARE !C CALCULATED FROM THE EQUATION R=G*PI*R**2*CBAR*N/(1+3G.R/(4.L)), WHERE !C G=GAMMA, R=PARTICLE RADIUS, CBAR=MEAN MOLECULAR SPEED, L=MEAN FREE !C PATH (TURCO ET AL, JGR, 94, 16493, 1989). !C RCLNO3HCLICE= & GCLNO3HCLICE*4.56E4*SQRT(T/97.45)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GCLNO3HCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K2S(ILON,INIV) = RCLNO3HCLICE RCLNO3H2OICE= & GCLNO3H2OICE*4.56E4*SQRT(T/97.45)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GCLNO3H2OICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZH2O K1S(ILON,INIV) = RCLNO3H2OICE RN2O5H2OICE= & GN2O5H2OICE*4.56E4*SQRT(T/108.0)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GN2O5H2OICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZH2O K3S(ILON,INIV) = RN2O5H2OICE RN2O5HCLICE= & GN2O5HCLICE*4.56E4*SQRT(T/108.0)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GN2O5HCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K4S(ILON,INIV) = RN2O5HCLICE RHOCLHCLICE= & GHOCLHCLICE*4.56E4*SQRT(T/52.5)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GHOCLHCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K5S(ILON,INIV) = RHOCLHCLICE RBRNO3H2OICE= & GBRNO3H2OICE*4.56E4*SQRT(T/142.0)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GBRNO3H2OICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZH2O K6S(ILON,INIV) = RBRNO3H2OICE RHOBRHCLICE= & GHOBRHCLICE*4.56E4*SQRT(T/96.9)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GHOBRHCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHCL K7S(ILON,INIV) = RHOBRHCLICE RHOBRHBRICE= & GHOBRHBRICE*4.56E4*SQRT(T/96.9)*RICE(ILON,INIV)**2 & *NICE(ILON,INIV)/ & (1.0+3.3E4*GHOBRHBRICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & /ZHOBR K8S(ILON,INIV) = RHOBRHBRICE END IF END DO END DO DO ILON = 1, PLON DO iniv = nivbot, nivtop K1(ILON,INIV) = K1L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K1S(ILON,INIV) K2(ILON,INIV) = K2L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K2S(ILON,INIV) K3(ILON,INIV) = K3L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K3S(ILON,INIV) K4(ILON,INIV) = K4L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K4S(ILON,INIV) K5(ILON,INIV) = K5L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K5S(ILON,INIV) K6(ILON,INIV) = K6L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K6S(ILON,INIV) K7(ILON,INIV) = K7L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K7S(ILON,INIV) K8(ILON,INIV) = K8L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K8S(ILON,INIV) K9(ILON,INIV) = 0. END DO END DO DO ILON = 1, PLON DO iniv = nivbot, nivtop K1(ILON,INIV) = MAX(K1(ILON,INIV),0.) K2(ILON,INIV) = MAX(K2(ILON,INIV),0.) K3(ILON,INIV) = MAX(K3(ILON,INIV),0.) K4(ILON,INIV) = MAX(K4(ILON,INIV),0.) K5(ILON,INIV) = MAX(K5(ILON,INIV),0.) K6(ILON,INIV) = MAX(K6(ILON,INIV),0.) K7(ILON,INIV) = MAX(K7(ILON,INIV),0.) K8(ILON,INIV) = MAX(K8(ILON,INIV),0.) K9(ILON,INIV) = MAX(K9(ILON,INIV),0.) END DO END DO RETURN END SUBROUTINE hetero SUBROUTINE density(WS,WN,TEMP,DENS) !C !C DENSITY OF TERNARY SOLUTION IN G/CM3 !C WS ,WN ARE WT FRACTION, !C FITTED TO 0.05