SUBROUTINE watercloud(ngrid,nlay, ptimestep, & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, & pq,pdq,pdqcloud,pdqscloud,pdtcloud, & nq,naersize,tau, & ccn,rdust,rice,nuice) IMPLICIT NONE c======================================================================= c Treatment of saturation of water vapor c c c Modif de zq si saturation dans l'atmosphere c si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l) c Le test est effectue de bas en haut. L'eau condensee c (si saturation) est remise dans la couche en dessous. c L'eau condensee dans la couche du bas est deposee a la surface c c Modification: Franck Montmessin water ice scheme c Francois Forget : change nuclei density & outputs c Ehouarn Millour: sept.2008, tracers are now handled c by name (and not fixed index) c c======================================================================= c----------------------------------------------------------------------- c declarations: c ------------- #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" #include "tracer.h" #include "comgeomfi.h" c Inputs: c ------ INTEGER ngrid,nlay 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 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 pq(ngrid,nlay,nq) ! traceur (kg/kg) real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1) integer nq ! nombre de traceurs integer naersize ! nombre de traceurs radiativement actifs (=naerkind) REAL tau(ngridmx,naersize) REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei ! (particules kg-1) real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m) c Outputs: c ------- real pdqcloud(ngrid,nlay,nq) ! tendance 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 rice(ngrid,nlay) ! Ice mass mean radius (m) ! (r_c in montmessin_2004) REAL nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution c local: c ------ REAL CBRT EXTERNAL CBRT INTEGER ig,l REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers REAL zqsat(ngridmx,nlayermx) ! saturation REAL zt(ngridmx,nlayermx) ! local value of temperature REAL masse (ngridmx,nlayermx) REAL epaisseur (ngridmx,nlayermx) ! REAL rfinal ! Ice crystal radius after condensation(m) REAL*8 seq ! Equilibrium saturation ration (accounting for curvature effect) REAL*8 dzq ! masse de glace echangee (kg/kg) REAL lw !Latent heat of sublimation (J.kg-1) REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K REAL Ctot REAL*8 ph2o,satu REAL*8 gr,Cste,up,dwn,newvap LOGICAL,SAVE :: firstcall=.true. ! To use more refined microphysics, set improved to .true. LOGICAL,PARAMETER :: improved=.true. c Pour diagnostique : c ~~~~~~~~~~~~~~~~~ c REAL icetot(ngridmx) ! Total mass of water ice (kg/m2) c REAL rave(ngridmx) ! Mean crystal radius in a column (m) ! indexes of water vapour, water ice and dust tracers: INTEGER,SAVE :: i_h2o=0 ! water vapour INTEGER,SAVE :: i_ice=0 ! water ice 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 i_h2o=igcm_h2o_vap i_ice=igcm_h2o_ice write(*,*) "watercloud: i_h2o=",i_h2o write(*,*) " i_ice=",i_ice firstcall=.false. ENDIF ! of IF (firstcall) c----------------------------------------------------------------------- c 1. initialisation c ----------------- c On "update" la valeur de q(nqmx) (water vapor) et temperature. c On effectue qqes calculs preliminaires sur les couches : c masse (kg.m-2), epaisseur(m). do l=1,nlay do ig=1,ngrid zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004 zq0(ig,l,i_h2o)=zq(ig,l,i_h2o) zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004 zq0(ig,l,i_ice)=zq(ig,l,i_ice) c This typical profile is not used anymore; rdust is c set up in updatereffrad.F. c rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9) enddo enddo do l=1,nlay do ig=1,ngrid c Calcul du rayon moyen des particules de glace. c Hypothese : Dans une couche, la glace presente se c repartit uniformement autour du nbre de poussieres c dont le rayon moyen est prescrit par rdust. c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+ & ccn(ig,l)*(4./3.)*pi*rdust(ig,l)**3. ) & / (ccn(ig,l)*4./3.*pi) ) rice(ig,l)=max(rice(ig,l),rdust(ig,l)) c Effective variance of the size distribution nuice(ig,l)=nuice_ref enddo ! of do ig=1,ngrid enddo ! of dol=1,nlay pdqscloud(1:ngrid,1:nq)=0 pdqcloud(1:ngrid,1:nlay,1:nq)=0 pdtcloud(1:ngrid,1:nlay)=0 c icetot(1:ngrid)=0 c rave(1:ngrid)=0 c ---------------------------------------------- c c c Rapport de melange a saturation dans la couche l : ------- c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call watersat(ngridmx*nlayermx,zt,pplay,zqsat) c taux de condensation (kg/kg/s-1) dans les differentes couches c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c Iceparty is not used anymore: water=>iceparty (JBM). c if(iceparty) then do l=1,nlay do ig=1,ngrid IF (improved) then c Improved microphysics scheme c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice) ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l) satu = zq(ig,l,i_h2o) / zqsat(ig,l) call growthrate(ptimestep,zt(ig,l),pplay(ig,l), & ph2o,ph2o/satu,seq,rice(ig,l),gr) Cste = ptimestep * 4. * pi * rice(ig,l) * * rho_ice * ccn(ig,l) up = zq(ig,l,i_h2o) + Cste * gr * seq dwn = 1. + Cste * gr / zqsat(ig,l) newvap = min(up/dwn,Ctot) gr = gr * ( newvap/zqsat(ig,l) - seq ) dzq = min( max( Cste * gr,-zq(ig,l,i_ice) ) * , zq(ig,l,i_h2o) ) c Nucleation (sat ratio must be larger than a critical value) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (satu.gt.1.) then if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8) * dzq = 0. endif ELSE c Old version c ~~~~~~~~~~~ if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then ! Condensation dzq=zq(ig,l,i_h2o)-zqsat(ig,l) elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then ! Sublimation dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice)) endif ENDIF ! of IF (improved) c Water Mass change c ~~~~~~~~~~~~~~~~~ zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice & +ccn(ig,l)*(4./3.)*pi*rdust(ig,l)**3.) & /(ccn(ig,l)*4./3.*pi) ), rdust(ig,l)) enddo ! of do ig=1,ngrid enddo ! of do l=1,nlay c The following part have been commented because iceparty c is not used anymore: water=>iceparty (JBM). c else ! if not iceparty c Saturation couche nlay a 2 : c ~~~~~~~~~~~~~~~~~~~~~~~~~~ c do l=nlay,2, -1 c do ig=1,ngrid c if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then c zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+ c & (zq(ig,l,i_h2o)-zqsat(ig,l)) c & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l)) c zq(ig,l,i_h2o)=zqsat(ig,l) c endif c enddo c enddo c Saturation couche l=1 si pas iceparty c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c do ig=1,ngridmx c if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then c pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1)) c & *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep) c zq(ig,1,i_h2o)=zqsat(ig,1) c endif c enddo c endif ! of if (iceparty) c Tendance finale c ~~~~~~~~~~~~~~~ do l=1, nlay do ig=1,ngridmx pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o) & -zq0(ig,l,i_h2o))/ptimestep pdqcloud(ig,l,i_ice) = & (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3 pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp end do end do 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************************************************** c Output c************************************************** ! NB: for diagnostics use zq(), the updated value of tracers c do ig=1,ngridmx c do l=1 ,nlay c masse de glace d'eau dans la couche l c icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice) c rayon moyen des cristaux dans la colonne ig c rave(ig)=rave(ig)+masse(ig,l)*zq(ig,l,i_ice)*rice(ig,l) c enddo c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) c if (icetot(ig)*1000.lt.0.01) rave(ig)=0. c enddo ! (ngridmx) c************************************************** RETURN END