subroutine n_ethane(ngrid,nq,nbin, * dt,pl,tl,aerad, * q,qprime) implicit none #include "dimensions.h" #include "microtab.h" #include "varmuphy.h" c Arguments c --------- integer ngrid,nq,nbin REAL dt ! physical time step (s) REAL pl(ngrid,nz) ! pressure at each level (mbar) REAL tl(ngrid,nz) ! temperature at each level (K) REAL aerad(nbin) ! Radius array c Tracers : REAL q(ngrid,nz,nq) ! tracer (kg/kg) REAL qprime(ngrid,nz,nbin) ! tracer (kg/kg) c Local variables c --------------- integer ntyp parameter (ntyp=3) real n_aer(nz,nbin,ntyp) real c2h6vap(nz) integer itrac integer ig,i,j,k,l,n ! Loop integers integer ilay,iq c Treatment c --------- DO ig = 1 , NGRID c Set up the aerosol array do j = 1, ntyp do k = 1, nbin itrac = (j-1) * nbin + k do l = 1, nz n_aer(l,k,j) = max(q(ig,l,itrac),0.) enddo enddo enddo c Set up the methane vapor array do l = 1, nz c2h6vap(l) = q(ig,l,nq) enddo call nucleacond2(ngrid,nbin,dt,ig,pl,tl,aerad, & n_aer,qprime,c2h6vap) c Update q arrays do j = 1, ntyp do k = 1, nbin itrac = (j-1) * nbin + k do l = 1, nz q(ig,l,itrac) = n_aer(l,k,j) enddo enddo enddo c Update methane vapor array do l = 1, nz q(ig,l,nq) = c2h6vap(l) c if(q(ig,l,nq).le.0.) print*,'in nuc2h6',ig,l, q(ig,l,nq) enddo ENDDO return END **************************************************************** subroutine nucleacond2(ngrid,nbin,dt,ig, * pl,tl,aerad,n_aer,qprime,c2h6vap) * * * This routine updates species concentrations due * * to both nucleation and condensation-induced variations. * * Gain and loss rates associated to each one of these * * processes are computed separately in other routines. * * * **************************************************************** implicit none #include "dimensions.h" #include "microtab.h" #include "varmuphy.h" integer ng,nalt parameter(ng=1,nalt=llm) real lv * dejà évalué dans n_methane ! REAL x1(nalt) REAL x2(nalt) REAL x3(nalt) REAL icefrac(nalt) REAL pmixch4(nalt) REAL pmixc2h6(nalt) REAL pmixn2(nalt) common/lheat/lv common/mixing/x1,x2,x3,icefrac, & pmixch4,pmixc2h6,pmixn2 c Arguments c --------- integer ngrid,nbin integer ig integer ntyp parameter (ntyp=3) real dt ! Global time step real pl(ngrid,nz),tl(ngrid,nz) real aerad(nbin) real c2h6vap(nz) ! Methane vapor mass mixing ratio (kg/m3) real c2h6vap_old real n_aer(nz,nbin,ntyp) ! number concentrations of particle/each size bin real qprime(ngrid,nz,nbin) ! tracer (kg/kg) REAL total1(nz),total11(nz),total2(nz),total22(nz) REAL dmsm,mtot c Local c ----- integer i,j,k,l,n,iindice,iselec real dQc ! Amount of condensed methane (kg/m3) during timestep real*8 sat_ratio ! Methane saturation ratio over liquid real*8 sat_ratmix ! Methane saturation ratio over liquid real*8 pc2h6 ! Methane partial pressure (Pa) real qsat ! Methane mass mixing ratio at saturation (kg/kg of air) real qsatmix ! Methane mass mixing ratio at saturation (kg/kg of air) real*8 rate(nbin) ! Heterogeneous Nucleation rate (s-1) real*8 elim real nsav(nbin,ntyp) real dn(nbin,ntyp) real rad(nbin) ! Radius of droplets in each size bin real*8 gr(nbin) ! Growth rate in each bin real radius ! Radius of droplets after growth real Qs ! Mass of condensate required to reach saturation real newsat real vol(nbin) real press real sig2,temp,seq(nbin) real Ctot,up,dwn,newvap,gltot real temp0,temp1,temp2,last_temp real qsat1,sat_ratio1,tempf(0:10),sat_ratiof(0:10) real rho_a,cap real tempref real xtime,xtime_prime c Variables for latent heat release real lw data lw / 581.e+3/ save lw c Treatment c --------- do i = 1, nbin vol(i) = 4./3. * pi * aerad(i)**3. enddo do l = 1, nz total1(l)=0. !solide do k = 1, nbin total1(l)=total1(l)+n_aer(l,k,2)*rhoi_c2h6 enddo total2(l)=c2h6vap(l) enddo c Start loop over heights DO 100 l = 1, nz iindice=0 ! mettre l'indice à 0 temp = tl(ig,l) press = pl(ig,l) tempref=temp c Save the values of the particle arrays before condensation do j = 1, ntyp do i = 1, nbin nsav(i,j) = n_aer(l,i,j) enddo enddo 99 continue call c2h6sat(temp,press,qsat) qsatmix=qsat c quantité pmixc2h6(l) déjà calculé dans cnuages.F et passé dans un common c qsat=pmixc2h6(l)*30.0 / (28.0*p) c qsatmix=pmixc2h6(l)*30.0 / (28.0*p) c Get the partial presure of methane vapor and its saturation ratio pc2h6 = c2h6vap(l) * (Mn2/Mc2h6) * press sat_ratio = c2h6vap(l) / qsat sat_ratmix = c2h6vap(l) / qsatmix c Get the rates of nucleation call nuclea2(nbin,aerad,pc2h6,temp,sat_ratio,rate) c Get the growth rates of condensation/sublimation up = c2h6vap(l) dwn = 1. Ctot = c2h6vap(l) DO i = 1, nbin if (n_aer(l,i,3).eq.0) then rad(i) = aerad(i) else rad(i) = ((n_aer(l,i,2)/n_aer(l,i,3) + & qprime(ig,l,i)/n_aer(l,i,3) & +vol(i))*0.75/pi)**(1./3.) endif * Equilibrium saturation ratio (due to curvature effect) seq(i) = exp( 2.*sig2(temp)*Mc2h6 /(rhoi_c2h6*rgp*temp*rad(i))) call growthrate2(dt,temp,press,pc2h6, & sat_ratmix,seq(i),rad(i),gr(i)) up = up + dt * gr(i) * 4. * pi * rhoi_c2h6 * rad(i) * seq(i) * * nsav(i,3) dwn= dwn+ dt * gr(i) * 4. * pi * rhoi_c2h6 * rad(i) / qsat * * nsav(i,3) Ctot= Ctot + rhoi_c2h6 * nsav(i,2) ENDDO newvap = min(up/dwn,Ctot) newvap = max(newvap,0.) gltot = 0. DO i = 1, nbin gr(i) = gr(i) * ( newvap/qsat - seq(i) ) if(nsav(i,2).le.0. .and. gr(i).le.0.) then n_aer(l,i,2) = 0. else n_aer(l,i,2) = nsav(i,2) + dt * gr(i) * 4. * pi * rad(i) * * n_aer(l,i,3) if (n_aer(l,i,2).le.0.) then n_aer(l,i,1) = n_aer(l,i,1) + n_aer(l,i,3) n_aer(l,i,2) = 0. n_aer(l,i,3) = 0. endif gltot=n_aer(l,i,2)*rhoi_c2h6+gltot endif ENDDO c Determine the mass of exchanged methane dQc = 0. DO i = 1, nbin dQc = dQc - rhoi_c2h6 * ( n_aer(l,i,2) - nsav(i,2) ) ENDDO c Update the methane vapor mixing ratio implied by c the cond/eva processes. c Arrays resetted to their initial value before condensation do j = 1, ntyp do i = 1, nbin dn(i,j) = n_aer(l,i,j) - nsav(i,j) n_aer(l,i,j) = nsav(i,j) enddo enddo c Update the c arrays. c nucleation & cond/eva tendencies added together. do i=1,nbin elim = dt * rate(i) n_aer(l,i,1) = n_aer(l,i,1) / (1.+elim) n_aer(l,i,3) = n_aer(l,i,3) + elim * n_aer(l,i,1) + dn(i,3) n_aer(l,i,1) = n_aer(l,i,1) + dn(i,1) n_aer(l,i,2) = n_aer(l,i,2) + dn(i,2) if(n_aer(l,i,2).lt.0.) n_aer(l,i,2)=0. enddo dQc = 0. DO i = 1, nbin dQc = dQc - rhoi_c2h6 * ( n_aer(l,i,2) - nsav(i,2) ) ENDDO c2h6vap(l) = c2h6vap(l) + dQc 100 CONTINUE do l = 1, nz total11(l)=0. do k = 1, nbin total11(l)=total11(l)+n_aer(l,k,2)*rhoi_c2h6 enddo total22(l)=c2h6vap(l) enddo return end ******************************************************* * * subroutine nuclea2(nbin,aerad,pc2h6,temp,sat,nucrate) * This subroutine computes the nucleation rate * * as given in Pruppacher & Klett (1978) in the * * case of water ice forming on a solid substrate. * * Definition refined by Keese (jgr,1989) * * * ******************************************************* implicit none #include "dimensions.h" #include "microtab.h" #include "varmuphy.h" integer nbin real aerad(nbin) real*8 nucrate(nbin) real*8 pc2h6 real temp real*8 sat integer l,i real*8 nc2h6 real sig2 ! Water-ice/air surface tension (N.m) real*8 rstar ! Radius of the critical germ (m) real*8 gstar ! # of molecules forming a critical embryo real*8 x ! Ratio rstar/radius of the nucleating dust particle real fistar ! Activation energy required to form a critical embryo (J) real*8 zeldov ! Zeldovitch factor (no dim) real*8 fshape2 ! function defined at the end of the file real*8 deltaf real nus data nus/1.e+13/ ! Jump frequency of a molecule (s-1) real m0 data m0/4.983e-26/ ! Weight of a methane molecule (kg) real vo1 c data vo1/9.094e-29/ data vo1/5.4746e-5/ ! Volume of 1 mole real desorp data desorp/0.288e-19/ ! Activation energy for desorption of water on a dust-like substrate (J/molecule) real surfdif data surfdif/0.288e-20/! Estimated activation energy for surface diffusion of water molecules (J/molecule) IF (sat .GT. 1.) then ! minimum condition to activate nucleation nc2h6 = pc2h6 / kbz / temp rstar = 2. * sig2(temp) * vo1 / (rgp*temp*log(sat)) gstar = 4. * nav * pi * (rstar**3) / (3.*vo1) c Loop over size bins do i=1,nbin x = aerad(i) / rstar x = aerad(imono) / rstar ! r(5)=monomere fistar = (4./3.*pi) * sig2(temp) * (rstar**2.) & *fshape2(mtetac2h6,x) deltaf = min( max((2.*desorp-surfdif-fistar)/(kbz*temp) & , -100.), 100.) if (deltaf.eq.-100.) then nucrate(i) = 0. else zeldov = sqrt ( fistar / (3.*pi*kbz*temp*(gstar**2.)) ) nucrate(i) = zeldov * kbz* temp * rstar**2. & * 4. * pi * ( nc2h6*aerad(i) )**2. & / ( fshape2(mtetac2h6,x) * nus * m0 ) & * dexp(deltaf) if(i.gt.imono) nucrate(i)= zeldov * kbz* temp * rstar**2. & * 4. * pi * vrat_e**(i-imono)*(nc2h6*aerad(imono) )**2. & / (fshape2(mtetac2h6,x) * nus * m0 ) & * dexp(deltaf) endif enddo ELSE do i=1,nbin nucrate(i) = 0. enddo ENDIF return end ****************************************************************** subroutine c2h6sat(t,p,qsat) * * * cette fonction calcule la pression de vapeur saturante de l' * * ethane a une altitude donnee z par Reid et al., p657 * * * * Compatible avec Barth et al., dans l'intervalle 30-90K * * * * * ****************************************************************** real rgp data rgp/8.3143/ * declaration des variables internes * ---------------------------------- real qsat,t,p c qsat=10.01-1085.0/(t-0.561) c qsat=10.**qsat/760.*1.e5 a=-6.34307 b=1.011630 c=-1.19116 d=-2.03539 pc=48.8*1.e5 tc=305.4 x=(1.-t/tc) if(x.gt. 0.) qsat=(1-x)**(-1)*(a*x+b*x**1.5+c*x**3.+d*x**6.) if(x.le. 0.) qsat=a*x/abs(1.-x) ! approx pour t > tc qsat=pc*exp(qsat) qsat=qsat* 30.0 / (28.0*p) ! kg/kg return end c======================================================================= subroutine growthrate2(timestep,temp,press,pc2h6,sat,seq,r,Cste) c c Determination of the droplet growth rate c c======================================================================= IMPLICIT NONE #include "dimensions.h" #include "microtab.h" #include "varmuphy.h" c----------------------------------------------------------------------- C DECLARATIONS: c ------------- common/lheat/Lv c c arguments: c ---------- REAL timestep REAL temp ! temperature in the middle of the layer (K) REAL press ! pressure in the middle of the layer (K) REAL*8 pc2h6 ! Methane vapor partial pressure (Pa) REAL*8 sat ! saturation ratio REAL r ! crystal radius before condensation (m) REAL seq ! Equilibrium saturation ratio c local: c ------ REAL psat REAL moln2,molc2h6 REAL To,tc2h6,wc2h6 ! Reid et al., (eq 7-9.4 + Appendix compound [168]) REAL fte c Effective gas molecular radius (m) data moln2/1.75e-10/ ! N2 c Effective gas molecular radius (m) data molc2h6/2.22e-10/ ! C2H6 c Temperature critique + omega data tc2h6/305.4/ data wc2h6/9.9e-2/ REAL k,Lv REAL knudsen ! Knudsen number (gas mean free path/particle radius) REAL a,Dv,lambda,Rk,Rd ! Intermediate computations for growth rate REAL*8 Cste c----------------------------------------------------------------------- c Ice particle growth rate by diffusion/impegement of molecules c r.dr/dt = (S-Seq) / (Seq*Rk+Rd) c with r the crystal radius, Rk and Rd the resistances due to c latent heat release and to vapor diffusion respectively c----------------------------------------------------------------------- psat = pc2h6 / sat c - Thermal conductibility of N2 k = ( 2.857e-2 * temp - 0.5428 ) * 4.184e-3 c - Latent heat of c2h6 (J.kg-1) Lv =581.e3 ! eq (7-9.4) Reid et al. fte=(1.-temp/tc2h6) if (fte.le.1.e-3) fte=1.e-3 Lv=8.314*tc2h6*(7.08*fte**0.354+10.95*wc2h6*fte**0.456)/30.e-3 c - Constant to compute gas mean free path c l= (T/P)*a, with a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) a = 0.707*rgp/(4 * pi* (moln2*1.e10)**2 * (nav*1.e-20)) c - Compute Dv, methane vapor diffusion coefficient c accounting for both kinetic and continuum regime of diffusion, c the nature of which depending on the Knudsen number. Dv = 1./3. * sqrt( 8*rgp*temp/(pi*Mc2h6) )* (kbz*1.e20) * temp / & (pi*press*(moln2*1.e10+molc2h6*1.e10)**2 * sqrt(1.+Mc2h6/Mn2) ) knudsen = temp / press * a / r lambda = (1.333+0.71/knudsen) / (1.+1./knudsen) Dv = Dv / (1. + lambda * knudsen) c - Compute Rk Rk = Lv**2 * rhoi_c2h6 * Mc2h6 / (k*rgp*temp**2.) c - Compute Rd Rd = rgp * temp *rhoi_c2h6 / (Dv*psat*Mc2h6) c - Compute: rdr/dt = Cste * (S-Seq) Cste = 1. / (seq*Rk+Rd) RETURN END ********************************************************* real function sig2(t) * this function computes the surface tension (N.m) * * between methane and air as a function of temp. * * Eq (12-3-6 et 12-3-7 Reid et al.) * ********************************************************* real t pc=48.8*1.01325e5 tc=305.4 tb=184.6 tr=t/tc tbr=tb/tc if(t.gt.305.0) then tr=305./tc endif sig2=0.1196*(1.+(tbr*alog(pc/1.01325))/(1.-tbr))-0.279 sig2=pc**(2./3.)*tc**(1./3.)*sig2*(1.-tr)**(11./9.) sig2=sig2*1.e-8 c sig2 = 21.157*((305.4-t)/(305.4-153.2))**(11./9.) c sig2=sig2*1.e-4 c c if(t.gt.305.0) then c sig2 = 21.157*((305.4-305.)/(305.4-153.2))**(11./9.) c sig2=sig2*1.e-4 c endif c c print*,t,sig2,sig3 return end ********************************************************* real*8 function fshape2(cost,rap) * function computing the f(m,x) factor * * related to energy required to form a critical embryo * ********************************************************* implicit none real cost real*8 rap real*8 phi real*8 a,b,c phi = sqrt( 1. - 2.*cost*rap + rap**2 ) a = 1. + ( (1.-cost*rap)/phi )**3 b = (rap**3) * (2.-3.*(rap-cost)/phi+((rap-cost)/phi)**3) c = 3. * cost * (rap**2) * ((rap-cost)/phi-1.) fshape2 = 0.5*(a+b+c) if (rap.gt.3000.) fshape2 = ((2.+cost)*(1.-cost)**2)/4. return end