SUBROUTINE ch4cloud(ngrid,nlay,naersize, ptimestep, & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, & pq,pdq,pdqcloud,pdqscloud,pdtcloud, & nq,rice_ch4) use comgeomfi_h IMPLICIT NONE c======================================================================= c Treatment of saturation of METHANE 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. CH4 condensee c (si saturation) est remise dans la couche en dessous. c Le methane condensee dans la couche du bas est deposee a la surface c c Melanie Vangvichith (adapted from Mars water clouds) c Tanguy Bertrand c Completely rewritten by Francois Forget (to properly estimate the c latent heat release. 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 CH4(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_ch4(ngrid,nlay) ! CH4 Ice mass mean radius (m) ! (r_c in montmessin_2004) c local: c ------ c REAL Nmix ! Cloud condensation nuclei c parameter (Nmix=1.E5) ! /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 zqsat(ngridmx,nlayermx) ! saturation REAL zt(ngridmx,nlayermx) ! local value of temperature REAL vecnull(ngridmx*nlayermx) LOGICAL,SAVE :: firstcall=.true. ! indexes of methane gas, methane ice and dust tracers: INTEGER,SAVE :: i_ch4=0 ! methane gas INTEGER,SAVE :: i_ice=0 ! methane ice REAL Tfin,qch4_fin REAL temp_fin ! function used to compute Tfin at the end of the timestep c TEMPORAIRE : real pdtcloudmax,pdtcloudmin integer igmin, igmax,lmin,lmax c ** un petit test de coherence c -------------------------- IF (firstcall) THEN IF(ngrid.NE.ngridmx) THEN PRINT*,'STOP dans ch4cloud' PRINT*,'probleme de dimensions :' PRINT*,'ngrid =',ngrid PRINT*,'ngridmx =',ngridmx STOP ENDIF if (nq.gt.nqmx) then write(*,*) 'stop in ch4cloud (nq.gt.nqmx)!' write(*,*) 'nq=',nq,' nqmx=',nqmx stop endif i_ch4=igcm_ch4_gas i_ice=igcm_ch4_ice write(*,*) "methanecloud: i_ch4=",i_ch4 write(*,*) " i_ice=",i_ice firstcall=.false. ENDIF ! of IF (firstcall) c----------------------------------------------------------------------- c 1. initialisation c ----------------- c On "update" la valeur de q(nqmx) (methane vapor) et temperature. c On effectue qqes calculs preliminaires sur les couches : do l=1,nlay do ig=1,ngrid zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep zq(ig,l,i_ch4)=pq(ig,l,i_ch4)+pdq(ig,l,i_ch4)*ptimestep zq(ig,l,i_ch4)=max(zq(ig,l,i_ch4),1.E-30) 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.) 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 q à saturation dans la couche l à temperature prédite zt c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call methanesat(ngridmx*nlayermx,zt,pplay,zqsat,vecnull) c Loop to work where CH4 condensation/sublimation is occuring c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do l=1,nlay do ig=1,ngrid if ((zq(ig,l,i_ch4).gt.zqsat(ig,l)).or. & (zq(ig,l,i_ice).gt.1.E-10) ) then ! phase change if(zq(ig,l,i_ice).lt.(zqsat(ig,l)-zq(ig,l,i_ch4))) then !case when everything is sublimed: qch4_fin = zq(ig,l,i_ch4) + zq(ig,l,i_ice) Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp else !case when CH4 ice is present at the end qch4_fin = zqsat(ig,l) Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp if(abs(Tfin-zt(ig,l)).gt.1.) then ! Improved calculation if the temperature change is significant (1K?) c Estimating the temperature Tfin and condensation at the end of the c timestep due to condensation-sublimation process, taking into acount c latent heat, given by: Tfin - zt = (zq - qsat(Tfin))*L/cp c Tfin, solution, of the equation, is estimated by iteration c We assume that improved Tfin is between zt and the 1st estimation of Tfin: Tfin = temp_fin(zt(ig,l),Tfin,0.01,pplay(ig,l), & zt(ig,l),zq(ig,l,i_ch4),qch4_fin) c Security check for small changes if(abs(zq(ig,l,i_ch4)-qch4_fin).gt. & abs(zq(ig,l,i_ch4)- zqsat(ig,l)) ) then write(*,*) "warning: inaccuracy in CH4 clouds" write(*,*) 'zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)', & zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l) qch4_fin = zqsat(ig,l) end if end if end if c Final tendancies pdqcloud(ig,l,i_ch4)=(qch4_fin-zq(ig,l,i_ch4))/ptimestep !TB16 IF (zq(ig,l,i_ch4)+pdqcloud(ig,l,i_ch4)*ptimestep.lt.0.) & then pdqcloud(ig,l,i_ch4)=-zq(ig,l,i_ch4)/ptimestep END IF pdqcloud(ig,l,i_ice) = - pdqcloud(ig,l,i_ch4) pdtcloud(ig,l)=(Tfin - zt(ig,l))/ptimestep c Ice particle radius zq(ig,l,i_ice)= & zq(ig,l,i_ice)+pdqcloud(ig,l,i_ice)*ptimestep rice_ch4(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ch4_ice & +Nmix_ch4*(4./3.)*pi*rnuclei**3.)/(Nmix_ch4*4./3.*pi)), & rnuclei) end if enddo ! of do ig=1,ngrid enddo ! of do l=1,nlay c A correction if a lot of subliming N2 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_ch4(ig,2)=rice_ch4(ig,3) rice_ch4(ig,1)=rice_ch4(ig,2) end if end do c************ ANALYSE pour Icarus **************** c c pdtcloudmax =0. c pdtcloudmin =0. c igmin=999 c igmax=999 c lmin =999 c lmax =999 c do l=1, nlay c do ig=1,ngridmx c if(pdtcloud(ig,l).gt.pdtcloudmax) then c igmax=ig c lmax = l c pdtcloudmax=pdtcloud(ig,l) c else if(pdtcloud(ig,l).lt.pdtcloudmin) then c igmin=ig c lmin = l c pdtcloudmin=pdtcloud(ig,l) c end if c end do c end do c c write(*,*) c write(*,*) "MAX , ig, l, Tbefore, Tafter, Tfin_ini" c write(*,*) igmax,lmax, zt(igmax,lmax), c & zt(igmax,lmax)+ pdtcloud(igmax,lmax)*ptimestep, c &zt(igmax,lmax)-(zqsat(igmax,lmax)-zq(igmax,lmax,i_ch4))*lw_ch4/cpp c write(*,*) "MAX , Qch4gas_before, Qch4gas_after, Qsat ini" c write(*,*) pq(igmax,lmax,i_ch4)+pdq(igmax,lmax,i_ch4)*ptimestep, c & pq(igmax,lmax,i_ch4)+ c & (pdq(igmax,lmax,i_ch4)+ pdqcloud(igmax,lmax,i_ch4))*ptimestep c & ,zqsat(igmax,lmax) c write(*,*) "MAX , Qch4ice_before, Qch4ice_after" c write(*,*) pq(igmax,lmax,i_ice)+pdq(igmax,lmax,i_ice)*ptimestep, c & pq(igmax,lmax,i_ice)+ c & (pdq(igmax,lmax,i_ice)+ pdqcloud(igmax,lmax,i_ice))*ptimestep c write(*,*) 'Rice=', rice_ch4(igmax,lmax) c c write(*,*) c write(*,*) "MIN , ig, l, Tbefore, Tafter, Tfin_ini" c write(*,*) igmin,lmin, zt(igmin,lmin), c & zt(igmin,lmin)+ pdtcloud(igmin,lmin)*ptimestep, c &zt(igmin,lmin)-(zqsat(igmin,lmin)-zq(igmin,lmin,i_ch4))*lw_ch4/cpp c write(*,*) "MIN , Qch4gas_before, Qch4gas_after, Qsat ini" c write(*,*) pq(igmin,lmin,i_ch4)+pdq(igmin,lmin,i_ch4)*ptimestep, c & pq(igmin,lmin,i_ch4)+ c & (pdq(igmin,lmin,i_ch4)+ pdqcloud(igmin,lmin,i_ch4))*ptimestep c & ,zqsat(igmin,lmin) c write(*,*) "MIN , Qch4ice_before, Qch4ice_after" c write(*,*) pq(igmin,lmin,i_ice)+pdq(igmin,lmin,i_ice)*ptimestep, c & pq(igmin,lmin,i_ice)+ c & (pdq(igmin,lmin,i_ice)+ pdqcloud(igmin,lmin,i_ice))*ptimestep c write(*,*) 'Rice=', rice_ch4(igmin,lmin) c write(*,*) "pdtcloud(igmin,lmin) " , pdtcloud(igmin,lmin) c write(*,*) "pdqcloud(i_ch4)", pdqcloud(igmin,lmin,i_ch4) c write(*,*) "pdqcloud(i_ice)", pdqcloud(igmin,lmin,i_ice) c c************ FIN ANALYSE pour Icarus **************** 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 !================================================================================ FUNCTION temp_fin(X1,X2,XACC,pres,temp,zq,qsat) implicit none #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "tracer.h" real pres, temp, zq,qsat real X1,X2,XACC, F real FMID,D,temp_fin,DX, XMID integer JMAX,J real func PARAMETER (JMAX=40) call methanesat(1,X2,pres,qsat,0.) FMID = X2 -temp -(zq-qsat)*lw_ch4/cpp call methanesat(1,X1,pres,qsat,0.) F = X1 -temp -(zq-qsat)*lw_ch4/cpp IF(F*FMID.GE.0.) write(*,*) 'Fix Tfin firt guesses in ch4cloud' IF(F.LT.0.)THEN temp_fin=X1 DX=X2-X1 ELSE temp_fin=X2 DX=X1-X2 ENDIF DO 11 J=1,JMAX DX=DX*.5 XMID=temp_fin+DX call methanesat(1,XMID,pres,qsat,0.) FMID = XMID -temp -(zq-qsat)*lw_ch4/cpp IF(FMID.LE.0.)temp_fin=XMID IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN 11 CONTINUE ! PAUSE 'too many bisections' END