SUBROUTINE ch4cloud(ngrid,nlay,naersize, ptimestep, & pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, & pq,pdq,pdqcloud,pdqscloud,pdtcloud, & nq,rice_ch4) use comgeomfi_h use comcstfi_mod, only: pi, g, cpp use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4 use callkeys_mod, only: Nmix_ch4 IMPLICIT NONE !======================================================================= ! Treatment of saturation of METHANE ! ! ! Modif de zq si saturation dans l'atmosphere ! si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l) ! Le test est effectue de bas en haut. CH4 condensee ! (si saturation) est remise dans la couche en dessous. ! Le methane condensee dans la couche du bas est deposee a la surface ! ! Melanie Vangvichith (adapted from Mars water clouds) ! Tanguy Bertrand ! Completely rewritten by Francois Forget (to properly estimate the ! latent heat release. ! !======================================================================= !----------------------------------------------------------------------- ! declarations: ! ------------- ! Inputs: ! ------ 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 ! Outputs: ! ------- 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) ! local: ! ------ ! REAL Nmix ! Cloud condensation nuclei ! 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(ngrid,nlay,nq) ! local value of tracers REAL zqsat(ngrid,nlay) ! saturation REAL zt(ngrid,nlay) ! local value of temperature REAL vecnull(ngrid*nlay) 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 ! TEMPORAIRE : real pdtcloudmax,pdtcloudmin integer igmin, igmax,lmin,lmax ! ** un petit test de coherence ! -------------------------- IF (firstcall) THEN IF(ngrid.NE.ngrid) THEN PRINT*,'STOP dans ch4cloud' PRINT*,'probleme de dimensions :' PRINT*,'ngrid =',ngrid PRINT*,'ngrid =',ngrid STOP ENDIF if (nq.gt.nq) then write(*,*) 'stop in ch4cloud (nq.gt.nq)!' write(*,*) 'nq=',nq,' nq=',nq stop endif i_ch4=igcm_ch4_gas i_ice=igcm_ch4_ice if (i_ch4.eq.0) then write(*,*) "stop in ch4cloud: i_ch4=0. " write(*,*) "Maybe put ch4 in traceur.def?" stop endif write(*,*) "methanecloud: i_ch4=",i_ch4 write(*,*) " i_ice=",i_ice firstcall=.false. ENDIF ! of IF (firstcall) !----------------------------------------------------------------------- ! 1. initialisation ! ----------------- ! On "update" la valeur de q(nq) (methane vapor) et temperature. ! 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 ! ---------------------------------------------- ! ! ! q à saturation dans la couche l à temperature prédite zt ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call methanesat(ngrid*nlay,zt,pplay,zqsat,vecnull) ! Loop to work where CH4 condensation/sublimation is occuring ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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?) ! Estimating the temperature Tfin and condensation at the end of the ! timestep due to condensation-sublimation process, taking into acount ! latent heat, given by: Tfin - zt = (zq - qsat(Tfin))*L/cp ! Tfin, solution, of the equation, is estimated by iteration ! 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) ! 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 ! 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 ! 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 ! A correction if a lot of subliming N2 fills the 1st layer FF04/2005 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Then that should not affect the ice particle radius do ig=1,ngrid 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 !************ ANALYSE pour Icarus **************** ! ! pdtcloudmax =0. ! pdtcloudmin =0. ! igmin=999 ! igmax=999 ! lmin =999 ! lmax =999 ! do l=1, nlay ! do ig=1,ngrid ! if(pdtcloud(ig,l).gt.pdtcloudmax) then ! igmax=ig ! lmax = l ! pdtcloudmax=pdtcloud(ig,l) ! else if(pdtcloud(ig,l).lt.pdtcloudmin) then ! igmin=ig ! lmin = l ! pdtcloudmin=pdtcloud(ig,l) ! end if ! end do ! end do ! ! write(*,*) ! write(*,*) "MAX , ig, l, Tbefore, Tafter, Tfin_ini" ! write(*,*) igmax,lmax, zt(igmax,lmax), ! & zt(igmax,lmax)+ pdtcloud(igmax,lmax)*ptimestep, ! &zt(igmax,lmax)-(zqsat(igmax,lmax)-zq(igmax,lmax,i_ch4))*lw_ch4/cpp ! write(*,*) "MAX , Qch4gas_before, Qch4gas_after, Qsat ini" ! write(*,*) pq(igmax,lmax,i_ch4)+pdq(igmax,lmax,i_ch4)*ptimestep, ! & pq(igmax,lmax,i_ch4)+ ! & (pdq(igmax,lmax,i_ch4)+ pdqcloud(igmax,lmax,i_ch4))*ptimestep ! & ,zqsat(igmax,lmax) ! write(*,*) "MAX , Qch4ice_before, Qch4ice_after" ! write(*,*) pq(igmax,lmax,i_ice)+pdq(igmax,lmax,i_ice)*ptimestep, ! & pq(igmax,lmax,i_ice)+ ! & (pdq(igmax,lmax,i_ice)+ pdqcloud(igmax,lmax,i_ice))*ptimestep ! write(*,*) 'Rice=', rice_ch4(igmax,lmax) ! ! write(*,*) ! write(*,*) "MIN , ig, l, Tbefore, Tafter, Tfin_ini" ! write(*,*) igmin,lmin, zt(igmin,lmin), ! & zt(igmin,lmin)+ pdtcloud(igmin,lmin)*ptimestep, ! &zt(igmin,lmin)-(zqsat(igmin,lmin)-zq(igmin,lmin,i_ch4))*lw_ch4/cpp ! write(*,*) "MIN , Qch4gas_before, Qch4gas_after, Qsat ini" ! write(*,*) pq(igmin,lmin,i_ch4)+pdq(igmin,lmin,i_ch4)*ptimestep, ! & pq(igmin,lmin,i_ch4)+ ! & (pdq(igmin,lmin,i_ch4)+ pdqcloud(igmin,lmin,i_ch4))*ptimestep ! & ,zqsat(igmin,lmin) ! write(*,*) "MIN , Qch4ice_before, Qch4ice_after" ! write(*,*) pq(igmin,lmin,i_ice)+pdq(igmin,lmin,i_ice)*ptimestep, ! & pq(igmin,lmin,i_ice)+ ! & (pdq(igmin,lmin,i_ice)+ pdqcloud(igmin,lmin,i_ice))*ptimestep ! write(*,*) 'Rice=', rice_ch4(igmin,lmin) ! write(*,*) "pdtcloud(igmin,lmin) " , pdtcloud(igmin,lmin) ! write(*,*) "pdqcloud(i_ch4)", pdqcloud(igmin,lmin,i_ch4) ! write(*,*) "pdqcloud(i_ice)", pdqcloud(igmin,lmin,i_ice) ! !************ FIN ANALYSE pour Icarus **************** !************************************************** ! Output --- removed !************************************************** ! NB: for diagnostics use zq(), the updated value of tracers ! Computing ext visible optical depth in each layer RETURN END !================================================================================ FUNCTION temp_fin(X1,X2,XACC,pres,temp,zq,qsat) use comcstfi_mod, only: pi, g, cpp use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4 implicit none 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