Index: trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F (revision 3353) +++ trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F (revision 3356) @@ -1,5 +1,5 @@ SUBROUTINE callsedim(ngrid,nlay, ptimestep, & pplev,zlev, pt, pdt, - & pq, pdqfi, pdqsed,pdqs_sed,nq) + & pq, pdqfi, pdqsed,pdqs_sed,nq,pphi) use radinc_h, only : naerkind @@ -39,4 +39,5 @@ real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature real,intent(in):: zlev(ngrid,nlay+1) ! altitude at layer boundaries + real,intent(in):: pphi(ngrid,nlay) ! geopotential integer,intent(in) :: nq ! number of tracers real,intent(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg) @@ -120,5 +121,5 @@ call newsedim(ngrid,nlay,1,ptimestep, & pplev,masse,epaisseur,zt,radius(iq),rho_q(iq), - & zqi(1,1,iq),wq,iq) + & zqi(1,1,iq),wq,iq,pphi) ! endif Index: trunk/LMDZ.PLUTO/libf/phypluto/newsedim.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/newsedim.F (revision 3353) +++ trunk/LMDZ.PLUTO/libf/phypluto/newsedim.F (revision 3356) @@ -1,7 +1,7 @@ SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep, - & pplev,masse,epaisseur,pt,rd,rho,pqi,wq,iq) + & pplev,masse,epaisseur,pt,rd,rho,pqi,wq,iq,pphi) use ioipsl_getin_p_mod, only: getin_p - use comcstfi_mod, only: r, g + use comcstfi_mod, only: r, g, rad use gases_h ! use tracer_h, only : igcm_h2o_ice @@ -41,4 +41,5 @@ real,intent(out) :: wq(ngrid,nlay+1) ! flux of tracer during timestep (?/m-2) integer,intent(in) :: iq ! tracer index + real,intent(in):: pphi(ngrid,nlay) ! geopotential c local: @@ -90,7 +91,7 @@ do igas=1, ngasmx if(gfrac(igas).ge.0.0) then - if(igas.eq.igas_N2) then - molrad = molrad + gfrac(igas)*2.2e-10 ! N2 - visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5 ! N2 + if(igas.eq.igas_CO2) then + molrad = molrad + gfrac(igas)*2.2e-10 ! CO2 + visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5 ! CO2 elseif(igas.eq.igas_N2) then molrad = molrad + gfrac(igas)*1.8e-10 ! N2 (Kunze et al. 2022) @@ -110,6 +111,6 @@ visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5 ! CH4 else - molrad = molrad + gfrac(igas)*2.2e-10 ! N2 by default - visc(:,:) = visc(:,:) + gfrac(igas)*1.e-5 ! N2 by default + molrad = molrad + gfrac(igas)*1.93e-10 ! N2 by default + visc(:,:) = visc(:,:) + 6.67e-6 ! N2 by default write(*,*) trim(gnom(igas))," is not included in" & ," newsedim, N2 is used by default" @@ -181,4 +182,7 @@ rsurf=rfall ! endif + + ! b = 2./9. * g + b = 2./9. * ((g*rad-pphi(ig,l))**2/(g*(rad**2))) ! AF24: from Pluto.old vstokes(ig,l) = b / visc(ig,l) * rho * rfall**3 / rsurf * Index: trunk/LMDZ.PLUTO/libf/phypluto/newsedim_pluto.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/newsedim_pluto.F (revision 3356) +++ trunk/LMDZ.PLUTO/libf/phypluto/newsedim_pluto.F (revision 3356) @@ -0,0 +1,184 @@ + SUBROUTINE newsedim_pluto(ngrid,nlay,naersize,ptimestep, + & pplev,masse,epaisseur,pt,rd,rho,pqi,wq,pphi) + + use comcstfi_mod, only: r, g, rad + IMPLICIT NONE +!================================================================== +! +! Purpose +! ------- +! Calculates sedimentation of 1 tracer +! of radius rd (m) and density rho (kg.m-3) +! +!================================================================== + +c----------------------------------------------------------------------- +c declarations: +c ------------- + +c +c arguments: +c ---------- + + INTEGER ngrid,nlay,naersize + REAL ptimestep ! pas de temps physique (s) + REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) + REAL pt(ngrid,nlay) ! temperature au centre des couches (K) + real masse (ngrid,nlay) ! masse d'une couche (kg) + real epaisseur (ngrid,nlay) ! epaisseur d'une couche (m) + real rd(naersize) ! particle radius (m) + real rho ! particle density (kg.m-3) + real pphi(ngrid,nlay) ! geeopotential + + +c Traceurs : + real pqi(ngrid,nlay) ! traceur (e.g. ?/kg) + real wq(ngrid,nlay+1) ! flux de traceur durant timestep (?/m-2) + +c local: +c ------ + + INTEGER l,ig, k, i + REAL rfall + + LOGICAL firstcall + SAVE firstcall + +c Traceurs : +c ~~~~~~~~ + real traversee(ngrid,nlay) + real vstokes(ngrid,nlay) + real w(ngrid,nlay) + real ptop, dztop, Ep, Stra + + DATA firstcall/.true./ + +c Physical constant +c ~~~~~~~~~~~~~~~~~ + REAL visc, molrad +c Gas molecular viscosity (N.s.m-2) + data visc/6.67e-6/ ! N2 TB14 +c Effective gas molecular radius (m) + data molrad/1.93e-10/ ! N2 TB14 +c local and saved variable + real a,b + save a,b + real b2 + +c ** un petit test de coherence +c -------------------------- + + IF (firstcall) THEN + IF(ngrid.NE.ngrid) THEN + PRINT*,'STOP dans newsedim' + PRINT*,'probleme de dimensions :' + PRINT*,'ngrid =',ngrid + PRINT*,'ngrid =',ngrid + STOP + ENDIF + firstcall=.false. + + +!======================================================================= +! Preliminary calculations for sedimenation velocity + +! - Constant to compute stokes speed simple formulae +! (Vstokes = b * rho* r**2 avec b= (2/9) * rho * g / visc + b = 2./9. * g / visc + +! - Constant to compute gas mean free path +! l= (T/P)*a, avec a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) + a = 0.707*8.31/(4*3.1416* molrad**2 * 6.023e23) + +c - Correction to account for non-spherical shape (Murphy et al. 1990) +c (correction = 0.85 for irregular particles, 0.5 for disk shaped particles) +c a = a * 0.85 + ENDIF + +c write(*,*) 'TB16 : stokes : g,visc,b,a,molrad ',g,visc,b,a,molrad + + + +c----------------------------------------------------------------------- +c 1. initialisation +c ----------------- + +c Sedimentation velocity (m/s) +c ~~~~~~~~~~~~~~~~~~~~~~ +c (stokes law corrected for low pressure by the Cunningham +c slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) + do l=1,nlay + do ig=1, ngrid + if (naersize.eq.1) then + rfall=rd(1) + else + i=ngrid*(l-1)+ig + rfall=rd(i) + endif + !vstokes(ig,l) = b * rho * rfall**2 * + b2=((g*rad-pphi(ig,l))**2/(g*(rad**2)))*2./9./visc + vstokes(ig,l) = b2 * rho * rfall**2 * + & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) + +c Layer crossing time (s) : + traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l) + end do + end do + +c Calcul de la masse d'atmosphere correspondant a q transferee +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c (e.g. on recherche le niveau en dessous de laquelle le traceur +c va traverser le niveau intercouche l : "dztop" est sa hauteur +c au dessus de l (m), "ptop" est sa pression (Pa)) + + do l=1,nlay + do ig=1, ngrid + + dztop = vstokes(ig,l)* ptimestep + Ep=0 + k=0 + +c ************************************************************** +c Simple Method + w(ig,l) = + & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g +cc write(*,*) 'OK simple method l,w =', l, w(ig,l) +cc write(*,*) 'OK simple method dztop =', dztop +c ************************************************************** +cccc Complex method : + if (dztop.gt.epaisseur(ig,l)) then +cccc Cas ou on "epuise" la couche l : On calcule le flux +cccc Venant de dessus en tenant compte de la variation de Vstokes + + Ep= epaisseur(ig,l) + Stra= traversee(ig,l) + do while(dztop.gt.Ep.and.l+k+1.le.nlay) + k=k+1 + dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra) + Ep = Ep + epaisseur(ig,l+k) + Stra = Stra + traversee(ig,l+k) + enddo + Ep = Ep - epaisseur(ig,l+k) + end if + ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) + w(ig,l) = (pplev(ig,l) -Ptop)/g +c +cc write(*,*) 'OK new method l,w =', l, w(ig,l) +cc write(*,*) 'OK new method dztop =', dztop +cc if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop +cc if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep +cc if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l) +cc if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l) +cc if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l) +c ************************************************************** + end do + end do + + call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq) +c write(*,*) ' newsed: wq(6), wq(7), q(6)', +c & wq(1,6),wq(1,7),pqi(1,6) + + + RETURN + END + Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3353) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3356) @@ -1512,5 +1512,5 @@ call callsedim(ngrid,nlayer,ptimestep, & pplev,zzlev,pt,pdt,pq,pdq, & - zdqsed,zdqssed,nq) + zdqsed,zdqssed,nq,pphi) endif