SUBROUTINE cocloud(ngrid,nlay,naersize, ptimestep, & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, & pq,pdq,pdqcloud,pdqscloud,pdtcloud, & nq,rice_co) use comgeomfi_h IMPLICIT NONE c======================================================================= c Treatment of saturation of CARBON MONOXIDE 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. CO condensee c (si saturation) est remise dans la couche en dessous. c CO condensee dans la couche du bas est deposee a la surface c c c======================================================================= c----------------------------------------------------------------------- c declarations: c ------------- #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" #include "tracer.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 naersize ! nombre de traceurs radiativement actifs (=naerkind) integer nq ! nombre de traceurs c Outputs: c ------- real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation CO(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_co(ngrid,nlay) ! Ice mass mean radius (m) ! (r_c in montmessin_2004) c local: c ------ c REAL Nmix ! Cloud condensation nuclei c parameter (Nmix=1.E2) ! /kg ! parameter (Nmix=1) ! /kg real rnuclei ! Nuclei geometric mean radius (m) parameter (rnuclei=2.E-7) ! m 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 vecnull(ngridmx*nlayermx) REAL masse (ngridmx,nlayermx) REAL epaisseur (ngridmx,nlayermx) ! REAL rfinal ! Ice crystal radius after condensation(m) REAL*8 dzq ! masse de glace echangee (kg/kg) REAL lw !Latent heat of sublimation (J.kg-1) LOGICAL,SAVE :: firstcall=.true. ! indexes of co gas, co ice and dust tracers: INTEGER,SAVE :: i_co=0 ! co gas INTEGER,SAVE :: i_ice=0 ! co ice c ** un petit test de coherence c -------------------------- IF (firstcall) THEN IF(ngrid.NE.ngridmx) THEN PRINT*,'STOP dans cocloud' PRINT*,'probleme de dimensions :' PRINT*,'ngrid =',ngrid PRINT*,'ngridmx =',ngridmx STOP ENDIF if (nq.gt.nqmx) then write(*,*) 'stop in cocloud (nq.gt.nqmx)!' write(*,*) 'nq=',nq,' nqmx=',nqmx stop endif c MELANIE : change these line i_co=igcm_co_gas i_ice=igcm_co_ice write(*,*) "cocloud: i_co=",i_co write(*,*) " i_ice=",i_ice firstcall=.false. ENDIF ! of IF (firstcall) c----------------------------------------------------------------------- c 1. initialisation c ----------------- c On "update" la valeur de q(nqmx) (co 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_co)=pq(ig,l,i_co)+pdq(ig,l,i_co)*ptimestep zq(ig,l,i_co)=max(zq(ig,l,i_co),1.E-30) ! FF 12/2004 zq0(ig,l,i_co)=zq(ig,l,i_co) 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) enddo enddo pdqscloud(1:ngrid,1:nq)=0 pdqcloud(1:ngrid,1:nlay,1:nq)=0 pdtcloud(1:ngrid,1:nlay)=0 c ---------------------------------------------- c c c Rapport de melange a saturation dans la couche l : ------- c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c do l=1,nlay c do ig=1,ngrid c call cosat(zt(ig,l),pplay(ig,l),zqsat(ig,l)) c write(101,*)'qsat',qsat(ig,l) c enddo c enddo call cosat(ngridmx*nlayermx,zt,pplay,zqsat,vecnull,vecnull) c TEMPORAIRE : c test sans condensation atmospherique c do l=1,nlay c do ig=1,ngrid c zqsat(ig,l) = zqsat(ig,l) *1000. c end do c end do c call WRITEDIAGFI(ngrid,"qsat_co","qsat_co","unit",3,zqsat) c do l=1,nlay c do ig=1,ngrid c zqsat(ig,l)=0.117*exp((16*568.7/8.314)*(1/90.7 c & -1/zt(ig,l)))*100000 c zqsat(ig,l)=(zqsat(ig,l)/pplay(ig,l))*(16/28) c write(106,*)'zqsat',zqsat(ig,l) c enddo ! of do ig=1,ngrid c enddo ! of do l=1,nlay c taux de condensation (kg/kg/s-1) dans les differentes couches c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do l=1,nlay do ig=1,ngrid c call cosat(zt(ig,l),pplay(ig,l),zqsat(ig,l)) if (zq(ig,l,i_co).ge.zqsat(ig,l))then ! Condensation dzq=zq(ig,l,i_co)-zqsat(ig,l) elseif(zq(ig,l,i_co).lt.zqsat(ig,l))then ! Sublimation dzq=-min(zqsat(ig,l)-zq(ig,l,i_co),zq(ig,l,i_ice)) endif c CO Mass change c ~~~~~~~~~~~~~~~~~ zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq zq(ig,l,i_co)=zq(ig,l,i_co)-dzq rice_co(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_co_ice & +Nmix_co*(4./3.)*pi*rnuclei**3.) & /(Nmix_co*4./3.*pi) ), rnuclei) ! CBRT=cube root enddo ! of do ig=1,ngrid enddo ! of do l=1,nlay c Saturation couche nlay a 2 : c ~~~~~~~~~~~~~~~~~~~~~~~~~~ c do l=nlay,2, -1 c do ig=1,ngrid c if (zq(ig,l,i_co).gt.zqsat(ig,l))then c zq(ig,l-1,i_co)= zq(ig,l-1,i_co)+ c & (zq(ig,l,i_co)-zqsat(ig,l)) c & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l)) c zq(ig,l,i_co)=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_co).gt.zqsat(ig,1))then c pdqscloud(ig,i_ice)=(zq(ig,1,i_co)-zqsat(ig,1)) c & *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep) c zq(ig,1,i_co)=zqsat(ig,1) c endif c enddo c Tendance finale c ~~~~~~~~~~~~~~~ do l=1, nlay do ig=1,ngridmx pdqcloud(ig,l,i_co)=(zq(ig,l,i_co) & -zq0(ig,l,i_co))/ptimestep pdqcloud(ig,l,i_ice) = & (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep lw=lw_co pdtcloud(ig,l)=-pdqcloud(ig,l,i_co)*lw/cpp end do end do c A correction if a lot of subliming co 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_co(ig,2)=rice_co(ig,3) rice_co(ig,1)=rice_co(ig,2) end if end do c************************************************** c Output --- removed c************************************************** ! NB: for diagnostics use zq(), the updated value of tracers c Computing ext visible optical depth in each layer RETURN END