Index: trunk/LMDZ.MARS/libf/phymars/watercloud.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/watercloud.F (revision 520) +++ trunk/LMDZ.MARS/libf/phymars/watercloud.F (revision 522) @@ -1,7 +1,9 @@ - SUBROUTINE watercloud(ngrid,nlay, ptimestep, + SUBROUTINE watercloud(ngrid,nlay,ptimestep, & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, & pq,pdq,pdqcloud,pdqscloud,pdtcloud, & nq,tau,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) +! to use 'getin' + USE ioipsl_getincom IMPLICIT NONE @@ -13,8 +15,11 @@ c - An improved microphysical scheme (see improvedclouds.F) c +c There is a time loop specific to cloud formation and sedimentation +c due to timescales smaller than the GCM integration timestep. +c c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, -c J.-B. Madeleine +c J.-B. Madeleine, Thomas Navarro c -c 2004 - Oct. 2011 +c 2004 - 2012 c======================================================================= @@ -35,28 +40,28 @@ INTEGER ngrid,nlay - integer nq ! nombre de traceurs + INTEGER nq ! nombre de traceurs REAL ptimestep ! pas de temps physique (s) REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) - REAL pdpsrf(ngrid) ! tendance surf pressure + REAL pdpsrf(ngrid) ! tendence surf pressure REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K) - REAL pdt(ngrid,nlay) ! tendance temperature des autres param. + REAL pdt(ngrid,nlay) ! tendence temperature des autres param. real pq(ngrid,nlay,nq) ! traceur (kg/kg) - real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1) - - REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point - REAL tauscaling(ngridmx) ! Convertion factor for dust amount - real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m) + real pdq(ngrid,nlay,nq) ! tendence avant condensation (kg/kg.s-1) + + REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point + REAL tauscaling(ngridmx) ! Convertion factor for dust amount + real rdust(ngridmx,nlay) ! Dust geometric mean radius (m) c Outputs: c ------- - real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1) + real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1) real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) - REAL pdtcloud(ngrid,nlay) ! tendance temperature due - ! a la chaleur latente + REAL pdtcloud(ngrid,nlay) ! tendence temperature due + ! a la chaleur latente REAL rice(ngrid,nlay) ! Ice mass mean radius (m) @@ -64,11 +69,42 @@ REAL nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution - real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius - real rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3) + real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius + real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3) c local: c ------ - INTEGER ig,l + ! for sedimentation + REAL zq(ngridmx,nlay,nqmx) ! local value of tracers + real masse (ngridmx,nlay) ! Layer mass (kg.m-2) + real epaisseur (ngridmx,nlay) ! Layer thickness (m) + real wq(ngridmx,nlay+1) ! displaced tracer mass (kg.m-2) + + ! for ice radius computation + REAL Mo,No + REAl ccntyp + real beta ! correction for the shape of the ice particles (cf. newsedim) + save beta + + ! for time loop + INTEGER microstep ! time subsampling step variable + INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation + SAVE imicro + REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation + SAVE microtimestep + + ! tendency given by clouds (inside the micro loop) + REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud + REAL subpdtcloud(ngrid,nlay) ! cf. pdtcloud + + ! global tendency (clouds+sedim+physics) + REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud + REAL subpdt(ngrid,nlay) ! cf. pdtcloud + + REAL CBRT + EXTERNAL CBRT + + + INTEGER iq,ig,l LOGICAL,SAVE :: firstcall=.true. @@ -93,21 +129,298 @@ write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap write(*,*) " igcm_h2o_ice=",igcm_h2o_ice + + + write(*,*) "correction for the shape of the ice particles ?" + beta=0.75 ! default value + call getin("ice_shape",beta) + write(*,*) " ice_shape = ",beta + + write(*,*) "time subsampling for microphysic ?" + imicro = 1 + call getin("imicro",imicro) + write(*,*)"imicro = ",imicro + + microtimestep = ptimestep/float(imicro) + write(*,*)"Physical timestep is",ptimestep + write(*,*)"Microphysics timestep is",microtimestep firstcall=.false. ENDIF ! of IF (firstcall) - -c Main call to the different cloud schemes: - IF (microphys) THEN - CALL improvedclouds(ngrid,nlay,ptimestep, - & pplev,pplay,pt,pdt, - & pq,pdq,pdqcloud,pdqscloud,pdtcloud, + +c-----Initialization + subpdq(1:ngrid,1:nlay,1:nq) = 0 + subpdt(1:ngrid,1:nlay) = 0 + pdqscloud(1:ngrid,1:nq) = 0 + zq(1:ngrid,1:nlay,1:nq) = 0 + +c-----Computing the different layer properties for clouds sedimentation + do l=1,nlay + do ig=1, ngrid + masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g + epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) + enddo + enddo + +c------------------------------------------------------------------ +c Time subsampling for coupled microphysics and sedimentation +c------------------------------------------------------------------ + DO microstep=1,imicro + +c------------------------------------------------------------------- +c 1. Tendencies: +c------------------ + + +c------ Temperature tendency subpdt + ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt + ! If imicro=1 subpdt is the same as pdt + DO l=1,nlay + DO ig=1,ngrid + subpdt(ig,l) = subpdt(ig,l) + & + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry + ENDDO + ENDDO +c------ Traceurs tendencies subpdq +c------ At each micro timestep we add pdq in order to have a stepped entry + IF (microphys) THEN + DO l=1,nlay + DO ig=1,ngrid + subpdq(ig,l,igcm_dust_mass) = + & subpdq(ig,l,igcm_dust_mass) + & + pdq(ig,l,igcm_dust_mass) + subpdq(ig,l,igcm_dust_number) = + & subpdq(ig,l,igcm_dust_number) + & + pdq(ig,l,igcm_dust_number) + subpdq(ig,l,igcm_ccn_mass) = + & subpdq(ig,l,igcm_ccn_mass) + & + pdq(ig,l,igcm_ccn_mass) + subpdq(ig,l,igcm_ccn_number) = + & subpdq(ig,l,igcm_ccn_number) + & + pdq(ig,l,igcm_ccn_number) + ENDDO + ENDDO + ENDIF + DO l=1,nlay + DO ig=1,ngrid + subpdq(ig,l,igcm_h2o_ice) = + & subpdq(ig,l,igcm_h2o_ice) + & + pdq(ig,l,igcm_h2o_ice) + subpdq(ig,l,igcm_h2o_vap) = + & subpdq(ig,l,igcm_h2o_vap) + & + pdq(ig,l,igcm_h2o_vap) + ENDDO + ENDDO + + +c------------------------------------------------------------------- +c 2. Main call to the different cloud schemes: +c------------------------------------------------ + IF (microphys) THEN + CALL improvedclouds(ngrid,nlay,microtimestep, + & pplev,pplay,pzlev,pt,subpdt, + & pq,subpdq,subpdqcloud,subpdtcloud, & nq,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) - ELSE - CALL simpleclouds(ngrid,nlay,ptimestep, - & pplev,pplay,pzlev,pzlay,pt,pdt, - & pq,pdq,pdqcloud,pdqscloud,pdtcloud, + ELSE + CALL simpleclouds(ngrid,nlay,microtimestep, + & pplev,pplay,pzlev,pzlay,pt,subpdt, + & pq,subpdq,subpdqcloud,subpdtcloud, & nq,tau,rice,nuice,rsedcloud) - ENDIF + ENDIF + +c-------------------------------------------------------------------- +c 3. CCN & ice sedimentation: +c-------------------------------- +! Sedimentation is done here for water ice and its CCN only +! callsedim in physics is done for dust (and others species if any) + + DO l=1,nlay + DO ig=1,ngrid + subpdt(ig,l) = + & subpdt(ig,l) + subpdtcloud(ig,l) + ENDDO + ENDDO + +c------- water ice update before sedimentation (radius is done in the cloud scheme routine) + DO l=1,nlay + DO ig=1, ngrid + zq(ig,l,igcm_h2o_ice)= max(1e-30, + & pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice) + & + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep) +! zq(ig,l,igcm_h2o_vap)= max(1e-30, +! & pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap) +! & + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep) + ENDDO + ENDDO + + +c------- CCN update before sedimentation + IF (microphys) THEN + DO l=1,nlay + DO ig=1,ngrid + zq(ig,l,igcm_ccn_number)= + & pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number) + & + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep + zq(ig,l,igcm_ccn_number)= max(1e-30, + & zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ? + zq(ig,l,igcm_ccn_mass)= + & pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass) + & + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep + zq(ig,l,igcm_ccn_mass)=max(1e-30, + & zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ? + ENDDO + ENDDO + ENDIF + + + + IF (microphys) THEN + +c------- CCN number sedimentation + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, + & rhocloud,zq(1,1,igcm_ccn_number),wq,beta) + do ig=1,ngrid + ! matters if one would like to know ccn surface deposition + pdqscloud(ig,igcm_ccn_number)= + & pdqscloud(ig,igcm_ccn_number) + wq(ig,1) + enddo + +c------- CCN mass sedimentation + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, + & rhocloud,zq(1,1,igcm_ccn_mass),wq,beta) + do ig=1,ngrid + ! matters if one would like to know ccn surface deposition + pdqscloud(ig,igcm_ccn_mass)= + & pdqscloud(ig,igcm_ccn_mass) + wq(ig,1) + enddo + +c------- Water ice sedimentation -- improved microphys + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, + & rhocloud,zq(1,1,igcm_h2o_ice),wq,beta) + ELSE + +c------- Water ice sedimentation -- simple microphys + call newsedim(ngrid,nlay,ngrid*nlay,1, + & microtimestep,pplev,masse,epaisseur,pt,rsedcloud, + & rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta) + + ENDIF + + +c------- Surface ice tendency update + DO ig=1,ngrid + pdqscloud(ig,igcm_h2o_ice)= + & pdqscloud(ig,igcm_h2o_ice) + wq(ig,1) + ENDDO + + +c------------------------------------------------------------------- +c 5. Updating tendencies after sedimentation: +c----------------------------------------------- + + DO l=1,nlay + DO ig=1,ngrid + + subpdq(ig,l,igcm_h2o_ice) = + & (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice)) + & /microtimestep + + + subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap) + & +subpdqcloud(ig,l,igcm_h2o_vap) + + ENDDO + ENDDO + + + IF (microphys) then + DO l=1,nlay + DO ig=1,ngrid + subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number) + & -pq(ig,l,igcm_ccn_number))/microtimestep + subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass) + & -pq(ig,l,igcm_ccn_mass))/microtimestep + ENDDO + ENDDO + ENDIF + + + + ENDDO ! of DO microstep=1,imicro + +c------------------------------------------------------------------- +c 6. Compute final tendencies after time loop: +c------------------------------------------------ + +c------ Whole temperature tendency pdtcloud + DO l=1,nlay + DO ig=1,ngrid + pdtcloud(ig,l) = + & subpdt(ig,l)/imicro-pdt(ig,l) + ENDDO + ENDDO +c------ Traceurs tendencies pdqcloud + DO iq=1,nq + DO l=1,nlay + DO ig=1,ngrid + pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro + & - pdq(ig,l,iq) + ENDDO + ENDDO + ENDDO +c------ Traceurs surface tendencies pdqscloud + DO iq=1,nq + DO ig=1,ngrid + pdqscloud(ig,iq) = + & pdqscloud(ig,iq)/ptimestep + ENDDO + ENDDO + + + +c------Update the ice particle size "rice" for output or photochemistry. +c------It is not used again in the water cycle. + IF(scavenging) THEN + DO l=1, nlay + DO ig=1,ngrid + Mo = pq(ig,l,igcm_h2o_ice) + & + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep + & + (pq(ig,l,igcm_ccn_mass) + & + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep) + & *tauscaling(ig) + 1.e-30 + No = (pq(ig,l,igcm_ccn_number) + & + pdqcloud(ig,l,igcm_ccn_number)*ptimestep) + & *tauscaling(ig) + 1.e-30 + rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + + & pdqcloud(ig,l,igcm_h2o_ice)*ptimestep) + & / Mo * rho_ice + & + (pq(ig,l,igcm_ccn_mass) + & + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep) + & *tauscaling(ig)/ Mo * rho_dust + rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) + rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.) + if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8 + ENDDO + ENDDO + ELSE + DO l=1,nlay + DO ig=1,ngrid + ccntyp = + & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.) + ccntyp = ccntyp /ccn_factor + rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice) + & + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice + & +ccntyp*(4./3.)*pi*rdust(ig,l)**3.) + & /(ccntyp*4./3.*pi) ), rdust(ig,l)) + ENDDO + ENDDO + ENDIF ! of IF(scavenging) + +!-------------------------------------------------------------- +!-------------------------------------------------------------- + c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 @@ -122,4 +435,33 @@ end do +c======================================================================= + +!!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice +!! if (photochem) then +!!c computation of dust and ice surface area (micron2/cm3) +!!c for heterogeneous chemistry +!! +!! do l = 1,nlay +!! do ig = 1,ngrid +!!c +!!c npart: number density of ccn in #/cm3 +!!c +!! npart(ig,l) = 1.e-6*ccn(ig,l) +!! $ *masse(ig,l)/epaisseur(ig,l) +!!c +!!c dust and ice surface area +!!c +!! surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2 +!!c +!! if (rice(ig,l) .ge. rdust(ig,l)) then +!! surfice(ig,l) = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2 +!! surfdust(ig,l) = 0. +!! else +!! surfice(ig,l) = 0. +!! end if +!! end do ! of do ig=1,ngrid +!! end do ! of do l=1,nlay +!! end if ! of photochem + RETURN END