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 c======================================================================= c Water-ice cloud formation c c Includes two different schemes: c - A simplified scheme (see simpleclouds.F) 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, Thomas Navarro c c 2004 - 2012 c======================================================================= c----------------------------------------------------------------------- c declarations: c ------------- #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" #include "tracer.h" #include "comgeomfi.h" #include "dimradmars.h" c Inputs: c ------ INTEGER ngrid,nlay 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) ! 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) ! tendence temperature des autres param. real pq(ngrid,nlay,nq) ! traceur (kg/kg) 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) ! 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) ! tendence temperature due ! a la chaleur latente REAL rice(ngrid,nlay) ! Ice mass mean radius (m) ! (r_c in montmessin_2004) REAL nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3) c local: c ------ ! 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. c ** un petit test de coherence c -------------------------- IF (firstcall) THEN IF(ngrid.NE.ngridmx) THEN PRINT*,'STOP dans watercloud' PRINT*,'probleme de dimensions :' PRINT*,'ngrid =',ngrid PRINT*,'ngridmx =',ngridmx STOP ENDIF if (nq.gt.nqmx) then write(*,*) 'stop in watercloud (nq.gt.nqmx)!' write(*,*) 'nq=',nq,' nqmx=',nqmx stop endif 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-----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,microtimestep, & pplev,pplay,pzlev,pzlay,pt,subpdt, & pq,subpdq,subpdqcloud,subpdtcloud, & nq,tau,rice,nuice,rsedcloud) 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) if ((Mo.lt.1.e-15) .or. (No.le.50)) then rice(ig,l) = 1.e-8 else !! AS: only perform computations if conditions not met rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.) endif 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)*rdust(ig,l)*rdust(ig,l)) & /(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 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c Then that should not affect the ice particle radius do ig=1,ngridmx if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) & rice(ig,2)=rice(ig,3) rice(ig,1)=rice(ig,2) end if end do c======================================================================= END