MODULE cloudth_mod IMPLICIT NONE CONTAINS SUBROUTINE cloudth(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) IMPLICIT NONE !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "nuage.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Gestion de deux versions de cloudth !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (iflag_cloudth_vert.GE.1) THEN CALL cloudth_vert(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:,:)=0. sigma2(:,:)=0. qlth(:,:)=0. qlenv(:,:)=0. qltot(:,:)=0. rneb(:,:)=0. qcloud(:)=0. cth(:,:)=0. cenv(:,:)=0. ctot(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) ! zqenv(ind1)=po(ind1) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*pi) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) endif ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif enddo return ! end END SUBROUTINE cloudth !=========================================================================== SUBROUTINE cloudth_vert(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== USE ioipsl_getin_p_mod, ONLY : getin_p IMPLICIT NONE #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "nuage.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv,sqrt2,sqrtpi REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL,SAVE :: vert_alpha !$OMP THREADPRIVATE(vert_alpha) LOGICAL, SAVE :: firstcall = .TRUE. !$OMP THREADPRIVATE(firstcall) REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf !------------------------------------------------------------------------------ ! Initialisation des variables r?elles !------------------------------------------------------------------------------ sigma1(:,:)=0. sigma2(:,:)=0. qlth(:,:)=0. qlenv(:,:)=0. qltot(:,:)=0. rneb(:,:)=0. qcloud(:)=0. cth(:,:)=0. cenv(:,:)=0. ctot(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) IF (firstcall) THEN vert_alpha=0.5 CALL getin_p('cloudth_vert_alpha',vert_alpha) WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha firstcall=.FALSE. ENDIF !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) ! zqenv(ind1)=po(ind1) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi=sqrt(2.*pi) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) ! deltasenv=aenv*0.01*po(ind1) ! deltasth=ath*0.01*zqta(ind1,ind2) xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlenv(ind1,ind2)=IntJ ! print*, qlenv(ind1,ind2),'VERIF EAU' IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) ! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2)) ! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1)) IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=IntJ ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ELSE IF (iflag_cloudth_vert == 2) THEN !------------------------------------------------------------------------------- ! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha*sigma2s xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) ! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) ! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2) IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2)) ! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) ! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) ! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlenv(ind1,ind2)=IntJ ! print*, qlenv(ind1,ind2),'VERIF EAU' IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2) IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=IntJ ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert==1 or 2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif enddo return ! end END SUBROUTINE cloudth_vert SUBROUTINE cloudth_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) IMPLICIT NONE !=========================================================================== ! Author : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "nuage.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL cth_vol(ngrid,klev) REAL cenv_vol(ngrid,klev) REAL ctot_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2 REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf IF (iflag_cloudth_vert.GE.1) THEN CALL cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:,:)=0. sigma2(:,:)=0. qlth(:,:)=0. qlenv(:,:)=0. qltot(:,:)=0. rneb(:,:)=0. qcloud(:)=0. cth(:,:)=0. cenv(:,:)=0. ctot(:,:)=0. cth_vol(:,:)=0. cenv_vol(:,:)=0. ctot_vol(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) !------------------------------------------------------------------------------- ! Cloud fraction in the thermals and standard deviation of the PDFs !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 !po = qt de l'environnement ET des thermique !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 ath=1./(1.+(alth*Lv/cppd)) !al, p84 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviations !------------------------------------------------------------------------------ ! tests ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) ! final option sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) !------------------------------------------------------------------------------ ! Condensed water and cloud cover !------------------------------------------------------------------------------ xth=sth/(sqrt2*sigma2s) xenv=senv/(sqrt2*sigma1s) cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) endif else ! Environnement only, follow the if l.110 zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sigma1s=ratqs(ind1,ind2)*zqenv(ind1) xenv=senv/(sqrt2*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183 enddo ! from the loop on ngrid l.108 return ! end END SUBROUTINE cloudth_v3 !=========================================================================== SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,t) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== USE ioipsl_getin_p_mod, ONLY : getin_p USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, & & cloudth_sigmath,cloudth_sigmaenv IMPLICIT NONE #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "nuage.h" INTEGER itap,ind1,ind2 INTEGER ngrid,klev,klon,l,ig REAL ztv(ngrid,klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid,klev) REAL fraca(ngrid,klev+1) REAL zpspsk(ngrid,klev) REAL paprs(ngrid,klev+1) REAL pplay(ngrid,klev) REAL ztla(ngrid,klev) REAL zthl(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zqsatenv(ngrid,klev) REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) REAL cth_vol(ngrid,klev) REAL cenv_vol(ngrid,klev) REAL ctot_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL,SAVE :: vert_alpha, vert_alpha_th !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th) REAL,SAVE :: sigma1s_factor=1.1 REAL,SAVE :: sigma1s_power=0.6 REAL,SAVE :: sigma2s_factor=0.09 REAL,SAVE :: sigma2s_power=0.5 REAL,SAVE :: cloudth_ratqsmin=-1. !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin) INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0 !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs) LOGICAL, SAVE :: firstcall = .TRUE. !$OMP THREADPRIVATE(firstcall) REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf REAL rhodz(ngrid,klev) REAL zrho(ngrid,klev) REAL dz(ngrid,klev) DO ind1 = 1, ngrid !Layer calculation rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2 zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3 dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre END DO !------------------------------------------------------------------------------ ! Initialize !------------------------------------------------------------------------------ sigma1(:,:)=0. sigma2(:,:)=0. qlth(:,:)=0. qlenv(:,:)=0. qltot(:,:)=0. rneb(:,:)=0. qcloud(:)=0. cth(:,:)=0. cenv(:,:)=0. ctot(:,:)=0. cth_vol(:,:)=0. cenv_vol(:,:)=0. ctot_vol(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) IF (firstcall) THEN vert_alpha=0.5 CALL getin_p('cloudth_vert_alpha',vert_alpha) WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha ! The factor used for the thermal is equal to that of the environment ! if nothing is explicitly specified in the def file vert_alpha_th=vert_alpha CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th) WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th ! Factor used in the calculation of sigma1s CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor) WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor ! Power used in the calculation of sigma1s CALL getin_p('cloudth_sigma1s_power',sigma1s_power) WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power ! Factor used in the calculation of sigma2s CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor) WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor ! Power used in the calculation of sigma2s CALL getin_p('cloudth_sigma2s_power',sigma2s_power) WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power ! Minimum value for the environmental air subgrid water distrib CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin) WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin ! Remove the dependency to ratqs from the variance of the vertical PDF CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs) WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs firstcall=.FALSE. ENDIF !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ecart-types des distributions !------------------------------------------------------------------------------- do ind1=1,ngrid if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 ath=1./(1.+(alth*Lv/cppd)) !al, p84 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviation !------------------------------------------------------------------------------ sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & & (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 ! sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*po(ind1) ELSE sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) ! tests ! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1) ! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) then ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs !------------------------------------------------------------------------------- deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s) xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s) xth1=(sth-deltasth)/(sqrt(2.)*sigma2s) xth2=(sth+deltasth)/(sqrt(2.)*sigma2s) coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1,ind2)=0.5*(1.+1.*erf(xth2)) cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ! Environment IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1)) IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! Thermal IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1)) IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2)) IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2)) IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1)) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ELSE IF (iflag_cloudth_vert >= 3) THEN IF (iflag_cloudth_vert < 5) THEN !------------------------------------------------------------------------------- ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) IF (iflag_cloudth_vert == 3) THEN deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha_th*sigma2s ELSE IF (iflag_cloudth_vert == 4) THEN IF (iflag_cloudth_vert_noratqs == 1) THEN deltasenv=vert_alpha*max(sigma1s_fraca,1e-10) deltasth=vert_alpha_th*sigma2s ELSE deltasenv=vert_alpha*sigma1s deltasth=vert_alpha_th*sigma2s ENDIF ENDIF xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) !CF_surfacique cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !CF_volumique & eau condense !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) if (deltasenv .lt. 1.e-10) then qlenv(ind1,ind2)=IntJ cenv_vol(ind1,ind2)=IntJ_CF else IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif !thermique IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) if (deltasth .lt. 1.e-10) then qlth(ind1,ind2)=IntJ cth_vol(ind1,ind2)=IntJ_CF else IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) ELSE IF (iflag_cloudth_vert == 5) THEN sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) !Environment sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) !Thermals !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) xth=sth/(sqrt(2.)*sigma2s) xenv=senv/(sqrt(2.)*sigma1s) !Volumique cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) !print *,'jeanjean_CV=',ctot_vol(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2)) qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !Surfacique !Neggers !beta=0.0044 !inverse_rho=1.+beta*dz(ind1,ind2) !print *,'jeanjean : beta=',beta !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !Brooks a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans shear !A_Maj_Brooks=0.17 !-- ARM LES !A_Maj_Brooks=0.18 !-- RICO LES !A_Maj_Brooks=0.19 !-- BOMEX LES Dx_Brooks=200000. f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) !print *,'jeanjean_f=',f_Brooks cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.)) cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.)) ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !print *,'JJ_ctot_1',ctot(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert<5) ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4) ! if (ctot(ind1,ind2).lt.1.e-10) then if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. ctot_vol(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif else ! gaussienne environnement seule zqenv(ind1)=po(ind1) Tbef=t(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) sth=0. sigma1s=ratqs(ind1,ind2)*zqenv(ind1) sigma2s=0. sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) if (ctot(ind1,ind2).lt.1.e-3) then ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else ! ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492 ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma1s cloudth_sigmath(ind1,ind2) = sigma2s enddo ! from the loop on ngrid l.333 return ! end END SUBROUTINE cloudth_vert_v3 ! SUBROUTINE cloudth_v6(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & & qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,T) USE ioipsl_getin_p_mod, ONLY : getin_p USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, & & cloudth_sigmath,cloudth_sigmaenv IMPLICIT NONE #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "nuage.h" !Domain variables INTEGER ngrid !indice Max lat-lon INTEGER klev !indice Max alt INTEGER ind1 !indice in [1:ngrid] INTEGER ind2 !indice in [1:klev] !thermal plume fraction REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox !temperatures REAL T(ngrid,klev) !temperature REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables) REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90) REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th) REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env) !pressure REAL paprs(ngrid,klev+1) !pressure at the interface of levels REAL pplay(ngrid,klev) !pressure at the middle of the level !humidity REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution REAL po(ngrid) !total water (qt) REAL zqenv(ngrid) !total water in the environment (qt_env) REAL zqta(ngrid,klev) !total water in the thermals (qt_th) REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th) REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env) REAL qlth(ngrid,klev) !condensed water in the thermals REAL qlenv(ngrid,klev) !condensed water in the environment REAL qltot(ngrid,klev) !condensed water in the gridbox !cloud fractions REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox !PDF of saturation deficit variables REAL rdd,cppd,Lv REAL Tbef,zdelta,qsatbef,zcor REAL alth,alenv,ath,aenv REAL sth,senv !saturation deficits in the thermals and environment REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF !cloud fraction variables REAL xth,xenv REAL inverse_rho,beta !Neggers et al. (2011) method REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method !Incloud total water variables REAL zqs(ngrid) !q_sat REAL qcloud(ngrid) !eau totale dans le nuage !Some arithmetic variables REAL erf,pi,sqrt2,sqrt2pi !Depth of the layer REAL dz(ngrid,klev) !epaisseur de la couche en metre REAL rhodz(ngrid,klev) REAL zrho(ngrid,klev) DO ind1 = 1, ngrid rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2] zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3] dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m] END DO !------------------------------------------------------------------------------ ! Initialization !------------------------------------------------------------------------------ qlth(:,:)=0. qlenv(:,:)=0. qltot(:,:)=0. cth_vol(:,:)=0. cenv_vol(:,:)=0. ctot_vol(:,:)=0. cth_surf(:,:)=0. cenv_surf(:,:)=0. ctot_surf(:,:)=0. qcloud(:)=0. rdd=287.04 cppd=1005.7 pi=3.14159 Lv=2.5e6 sqrt2=sqrt(2.) sqrt2pi=sqrt(2.*pi) DO ind1=1,ngrid !------------------------------------------------------------------------------- !Both thermal and environment in the gridbox !------------------------------------------------------------------------------- IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul de qsat_th !-------------------------------------------- Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatth(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_th !-------------------------------------------- alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 these Arnaud Jam ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations bi-Gaussian PDF !-------------------------------------------- sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) xth=sth/(sqrt2*sigma_th) xenv=senv/(sqrt2*sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth)) cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2)) qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- if (ctot_surf(ind1,ind2) .lt. 1.e-10) then ctot_vol(ind1,ind2)=0. ctot_surf(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1) endif !------------------------------------------------------------------------------- !Environment only in the gridbox !------------------------------------------------------------------------------- ELSE !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor zqsatenv(ind1,ind2)=qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations Gaussian PDF !-------------------------------------------- zqenv(ind1)=po(ind1) sigma_env=ratqs(ind1,ind2)*zqenv(ind1) xenv=senv/(sqrt2*sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2)) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks=0.6694 b_Brooks=0.1882 A_Maj_Brooks=0.1635 !-- sans dependence au shear Dx_Brooks=200000. f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks)) ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- if (ctot_surf(ind1,ind2) .lt. 1.e-8) then ctot_vol(ind1,ind2)=0. ctot_surf(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) else qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2) endif END IF ! From the separation (thermal/envrionnement) et (environnement only) ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma_env cloudth_sigmath(ind1,ind2) = sigma_th END DO ! From the loop on ngrid return END SUBROUTINE cloudth_v6 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points, & & temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl, & & ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS) ! Date: Adapted from cloudth_vert_v3 in 2021 ! Aim : computes qc and rneb in thermals with cold microphysical considerations ! + for mixed phase boundary layer clouds, calculate ql and qi from ! a stationary MPC model ! IMPORTANT NOTE: we assume iflag_clouth_vert=3 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE ioipsl_getin_p_mod, ONLY : getin_p USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF USE phys_local_var_mod, ONLY : qlth, qith, qsith IMPLICIT NONE #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "nuage.h" !------------------------------------------------------------------------------ ! Declaration !------------------------------------------------------------------------------ ! INPUT/OUTPUT INTEGER, INTENT(IN) :: klon,klev,ind2 REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! Temperature [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! Virtual potential temp [K] REAL, DIMENSION(klon), INTENT(IN) :: po ! specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: fraca ! Fraction of the mesh covered by thermals [0-1] REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! Pressure at layer interfaces [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! Pressure at the center of layers [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! Liquid temp [K] REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! Liquid potential temp [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib. REAL, DIMENSION(klon), INTENT(IN) :: zqs ! Saturation specific humidity in the mesh [kg/kg] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: snowflux ! snow flux at the interface of the layer [kg/m2/s] INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot ! Cloud fraction [0-1] REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg] REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1] ! LOCAL VARIABLES INTEGER itap,ind1,l,ig,iter,k INTEGER iflag_topthermals, niter REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon) REAL sigma1(klon,klev) REAL sigma2(klon,klev) REAL qcth(klon,klev) REAL qcenv(klon,klev) REAL qctot(klon,klev) REAL cth(klon,klev) REAL cenv(klon,klev) REAL cth_vol(klon,klev) REAL cenv_vol(klon,klev) REAL rneb(klon,klev) REAL zqenv(klon), zqth(klon), zqenvonly(klon) REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi REAL rdd,cppd,Lv REAL alth,alenv,ath,aenv REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL zdelta,qsatbef,zcor REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon) REAL qlbef REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon) REAL erf REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) REAL rhodz(klon,klev) REAL zrho(klon,klev) REAL dz(klon,klev) REAL qslenv(klon), qslth(klon) REAL alenvl, aenvl REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc REAL senvi, senvl, qbase, sbase, qliqth, qiceth REAL qmax, ttarget, stmp, cout, coutref, fraci REAL maxi, mini, pas, temp_lim REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1) ! Modifty the saturation deficit PDF in thermals ! in the presence of ice crystals REAL,SAVE :: C_mpc !$OMP THREADPRIVATE(C_mpc) REAL, SAVE :: Ni,C_cap,Ei,d_top !$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top) ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL,SAVE :: vert_alpha, vert_alpha_th !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th) REAL,SAVE :: sigma1s_factor=1.1 REAL,SAVE :: sigma1s_power=0.6 REAL,SAVE :: sigma2s_factor=0.09 REAL,SAVE :: sigma2s_power=0.5 REAL,SAVE :: cloudth_ratqsmin=-1. !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin) INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0 !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs) LOGICAL, SAVE :: firstcall = .TRUE. !$OMP THREADPRIVATE(firstcall) CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: routname = 'cloudth_mpc' !------------------------------------------------------------------------------ ! Initialisation !------------------------------------------------------------------------------ ! Few initial checksS IF (iflag_cloudth_vert.NE.3) THEN abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3' CALL abort_physic(routname,abort_message,1) ENDIF DO k = 1,klev DO ind1 = 1, klon rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2 zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3 dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre END DO END DO sigma1(:,:)=0. sigma2(:,:)=0. qcth(:,:)=0. qcenv(:,:)=0. qctot(:,:)=0. qlth(:,ind2)=0. qith(:,ind2)=0. rneb(:,:)=0. qcloud(:)=0. cth(:,:)=0. cenv(:,:)=0. ctot(:,:)=0. cth_vol(:,:)=0. cenv_vol(:,:)=0. ctot_vol(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. rdd=287.04 cppd=1005.7 pi=3.14159 sqrt2pi=sqrt(2.*pi) sqrt2=sqrt(2.) sqrtpi=sqrt(pi) icefrac(:,ind2)=0. althl=0. althi=0. athl=0. athi=0. senvl=0. senvi=0. sthi=0. sthl=0. sthil=0. IF (firstcall) THEN vert_alpha=0.5 CALL getin_p('cloudth_vert_alpha',vert_alpha) WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha ! The factor used for the thermal is equal to that of the environment ! if nothing is explicitly specified in the def file vert_alpha_th=vert_alpha CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th) WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th ! Factor used in the calculation of sigma1s CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor) WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor ! Power used in the calculation of sigma1s CALL getin_p('cloudth_sigma1s_power',sigma1s_power) WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power ! Factor used in the calculation of sigma2s CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor) WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor ! Power used in the calculation of sigma2s CALL getin_p('cloudth_sigma2s_power',sigma2s_power) WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power ! Minimum value for the environmental air subgrid water distrib CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin) WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin ! Remove the dependency to ratqs from the variance of the vertical PDF CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs) WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs ! Modifies the PDF in thermals when ice crystals are present C_mpc=1.e4 CALL getin_p('C_mpc',C_mpc) WRITE(*,*) 'C_mpc = ', C_mpc Ni=2.0e3 CALL getin_p('Ni', Ni) WRITE(*,*) 'Ni = ', Ni Ei=0.5 CALL getin_p('Ei', Ei) WRITE(*,*) 'Ei = ', Ei C_cap=0.5 CALL getin_p('C_cap', C_cap) WRITE(*,*) 'C_cap = ', C_cap d_top=1.2 CALL getin_p('d_top', d_top) WRITE(*,*) 'd_top = ', d_top firstcall=.FALSE. ENDIF !------------------------------------------------------------------------------- ! Identify grid points with potential mixed-phase conditions !------------------------------------------------------------------------------- temp_lim=RTT-40.0 DO ind1=1,klon IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) & .AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2) & .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN mpc_bl_points(ind1,ind2)=1 ELSE mpc_bl_points(ind1,ind2)=0 ENDIF ENDDO !------------------------------------------------------------------------------- ! Thermal fraction calculation and standard deviation of the distribution !------------------------------------------------------------------------------- ! calculation of temperature, humidity and saturation specific humidity Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2) Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2) Tbefenvonly(:)=temp(:,ind2) zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv zqth(:)=zqta(:,ind2) zqenvonly(:)=po(:) CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv) CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv) CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv) CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth) CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth) CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth) DO ind1=1,klon IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement ! Environment: IF (Tbefenv(ind1) .GE. RTT) THEN Lv=RLVTT ELSE Lv=RLSTT ENDIF alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 senv=aenv*(po(ind1)-zqsatenv(ind1)) !s, p84 ! For MPCs: IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2) aenvl=1./(1.+(alenv*Lv/cppd)) senvl=aenvl*(po(ind1)-qslenv(ind1)) senvi=senv ENDIF ! Thermals: IF (Tbefth(ind1) .GE. RTT) THEN Lv=RLVTT ELSE Lv=RLSTT ENDIF alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2) ath=1./(1.+(alth*Lv/cppd)) sth=ath*(zqta(ind1,ind2)-zqsatth(ind1)) ! For MPCs: IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN althi=alth althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2) athl=1./(1.+(alth*RLVTT/cppd)) athi=alth sthl=athl*(zqta(ind1,ind2)-qslth(ind1)) sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2)) sthil=athi*(zqta(ind1,ind2)-qslth(ind1)) ENDIF !------------------------------------------------------------------------------- ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s ! Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3 !------------------------------------------------------------------------------- IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC ! Standard deviation of the distributions sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & & (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*po(ind1) ELSE sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) deltasenv=aenv*vert_alpha*sigma1s deltasth=ath*vert_alpha_th*sigma2s xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) !surface CF cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !volume CF and condensed water !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) IF (deltasenv .LT. 1.e-10) THEN qcenv(ind1,ind2)=IntJ cenv_vol(ind1,ind2)=IntJ_CF ELSE IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF !thermals IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) IF (deltasth .LT. 1.e-10) THEN qcth(ind1,ind2)=IntJ cth_vol(ind1,ind2)=IntJ_CF ELSE IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN ctot(ind1,ind2)=0. ctot_vol(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1) qincloud(ind1)=0. ELSE qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2) ENDIF ELSE ! mpc_bl_points>0 ! Treat boundary layer mixed phase clouds ! thermals !========= ! ice phase !........... qiceth=0. deltazlev_mpc=dz(ind1,:) temp_mpc=ztla(ind1,:)*zpspsk(ind1,:) pres_mpc=pplay(ind1,:) fraca_mpc=fraca(ind1,:) snowf_mpc=snowflux(ind1,:) iflag_topthermals=0 IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN iflag_topthermals = 1 ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN iflag_topthermals = 2 ELSE iflag_topthermals = 0 ENDIF CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:),qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,snowf_mpc,fraca_mpc,qith(ind1,:)) ! qmax calculation sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) deltasth=athl*vert_alpha_th*sigma2s xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s) xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.) ! Liquid phase !................ ! We account for the effect of ice crystals in thermals on sthl ! and on the width of the distribution sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2)) & + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) deltasthc=athl*vert_alpha_th*sigma2sc xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc) xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1)) IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2) IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc) IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc) qliqth=IntJ+IntI1+IntI2+IntI3 qlth(ind1,ind2)=MAX(0.,qliqth) ! Condensed water qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2) ! consistency with subgrid distribution IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN fraci=qith(ind1,ind2)/qcth(ind1,ind2) qcth(ind1,ind2)=qmax qith(ind1,ind2)=fraci*qmax qlth(ind1,ind2)=(1.-fraci)*qmax ENDIF ! Cloud Fraction !............... ! calculation of qbase which is the value of the water vapor within mixed phase clouds ! such that the total water in cloud = qbase+qliqth+qiceth ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction ! sbase and qbase calculation (note that sbase is wrt liq so negative) ! look for an approximate solution with iteration ttarget=qcth(ind1,ind2) mini= athl*(qsith(ind1,ind2)-qslth(ind1)) maxi= 0. !athl*(3.*sqrt(sigma2s)) niter=20 pas=(maxi-mini)/niter stmp=mini sbase=stmp coutref=1.E6 DO iter=1,niter cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) & + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget) IF (cout .LT. coutref) THEN sbase=stmp coutref=cout ELSE stmp=stmp+pas ENDIF ENDDO qbase=MAX(0., sbase/athl+qslth(ind1)) ! surface cloud fraction in thermals cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s)) cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.) !volume cloud fraction in thermals !to be checked xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s) xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s) exp_xth1 = exp(-1.*xth1**2) exp_xth2 = exp(-1.*xth2**2) IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 IntJ_CF=0.5*(1.-1.*erf(xth2)) IF (deltasth .LT. 1.e-10) THEN cth_vol(ind1,ind2)=IntJ_CF ELSE IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.) ! Environment !============= ! In the environment/downdrafts, only liquid clouds ! See Shupe et al. 2008, JAS ! standard deviation of the distribution in the environment sth=sthl senv=senvl sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & & (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5 ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution ! in the environement sigma1s_ratqs=1E-10 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin*po(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs deltasenv=aenvl*vert_alpha*sigma1s xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s) xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s) exp_xenv1 = exp(-1.*xenv1**2) exp_xenv2 = exp(-1.*xenv2**2) !surface CF cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) !volume CF and condensed water IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 IntJ_CF=0.5*(1.-1.*erf(xenv2)) IF (deltasenv .LT. 1.e-10) THEN qcenv(ind1,ind2)=IntJ cenv_vol(ind1,ind2)=IntJ_CF ELSE IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF ENDIF qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.) cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.) ! Thermals + environment ! ======================= ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) IF (qcth(ind1,ind2) .GT. 0) THEN icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2)/(fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)) icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.) ELSE icefrac(ind1,ind2)=0. ENDIF IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN ctot(ind1,ind2)=0. ctot_vol(ind1,ind2)=0. qincloud(ind1)=0. qcloud(ind1)=zqsatenv(ind1) ELSE qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) & +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1)) qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) & +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.) ENDIF ENDIF ! mpc_bl_points ELSE ! gaussian for environment only IF (Tbefenvonly(ind1) .GE. RTT) THEN Lv=RLVTT ELSE Lv=RLSTT ENDIF zthl(ind1,ind2)=temp(ind1,ind2)*(101325./paprs(ind1,ind2))**(rdd/cppd) alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2) aenv=1./(1.+(alenv*Lv/cppd)) senv=aenv*(po(ind1)-zqsatenvonly(ind1)) sth=0. sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1) sigma2s=0. sqrt2pi=sqrt(2.*pi) xenv=senv/(sqrt(2.)*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) ctot_vol(ind1,ind2)=ctot(ind1,ind2) qctot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) IF (ctot(ind1,ind2).LT.1.e-3) THEN ctot(ind1,ind2)=0. qcloud(ind1)=zqsatenvonly(ind1) qincloud(ind1)=0. ELSE qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1) qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.) ENDIF ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492 ! Outputs used to check the PDFs cloudth_senv(ind1,ind2) = senv cloudth_sth(ind1,ind2) = sth cloudth_sigmaenv(ind1,ind2) = sigma1s cloudth_sigmath(ind1,ind2) = sigma2s ENDDO !loop on klon RETURN END SUBROUTINE cloudth_mpc !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,snowf,fraca,qith) ! parameterization of ice for boundary ! layer mixed-phase clouds assuming a stationary system ! ! Note that vapor deposition on ice crystals and riming of liquid droplets ! depend on the ice number concentration Ni ! One could assume that Ni depends on qi, e.g., Ni=beta*(qi*rho)**xi ! and use values from Hong et al. 2004, MWR for instance ! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962 ! One could also think of a more complex expression of Ni; ! function of qi, T, the concentration in aerosols or INP .. ! Here we prefer fixing Ni to a tuning parameter ! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard ! in Mioche et al. 2017 ! ! ! References: !------------ ! This parameterization is thoroughly described in Vignon et al. ! ! More specifically ! for the Water vapor deposition process: ! ! Rotstayn, L. D., 1997: A physically based scheme for the treat- ! ment of stratiform cloudfs and precipitation in large-scale ! models. I: Description and evaluation of the microphysical ! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282. ! ! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk ! stratiform cloud microphysics scheme in the NCAR Com- ! munity Atmosphere Model (CAM3). Part I: Description and ! numerical tests. J. Climate, 21, 3642–3659 ! ! for the Riming process: ! ! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro- ! scale structure and organization of clouds and precipitation in ! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’ ! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206 ! ! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit ! forecasts of winter precipitation using an improved bulk ! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit ! forecasts of winter precipitation using an improved bulk ! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542 ! ! For the formation of clouds by thermals: ! ! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of ! the Atmospheric Sciences, 65, 407–425. ! ! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a ! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3 ! ! ! ! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr !============================================================================= USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF USE ioipsl_getin_p_mod, ONLY : getin_p USE phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm IMPLICIT none INCLUDE "YOMCST.h" INCLUDE "nuage.h" INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ? REAL, INTENT(IN) :: Ni ! ice number concentration [m-3] REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency REAL, INTENT(IN) :: C_cap ! ice crystal capacitance REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa] REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg] REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg] REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg] REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m] REAL, DIMENSION(klev+1), INTENT(IN) :: snowf ! snow flux at the upper inferface REAL, DIMENSION(klev+1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg] INTEGER ind2p1,ind2p2 REAL rho(klev) REAL unsurtaudet, unsurtaustardep, unsurtaurim REAL qi, AA, BB, Ka, Dv, rhoi REAL p0, t0, fp1, fp2 REAL alpha, flux_term REAL det_term, precip_term, rim_term, dep_term ind2p1=ind2+1 ind2p2=ind2+2 rho=pres/temp/RD ! air density kg/m3 Ka=2.4e-2 ! thermal conductivity of the air, SI p0=101325.0 ! ref pressure T0=273.15 ! ref temp rhoi=500.0 ! cloud ice density following Reisner et al. 1998 alpha=700. ! fallvelocity param IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI ! Detrainment term: unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2) ! Deposition term AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.) BB=1./(rho(ind2)*Dv*qsith(ind2)) unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) & *4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33) ! Riming term neglected at this level !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4) qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12) qi=MAX(qi,0.)**(3./2.) ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2 Dv=0.0001*0.211*(p0/pres(ind2p1))*((temp(ind2p1)/T0)**1.94) ! water vapor diffusivity in air, SI ! Detrainment term: unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1) det_term=-unsurtaudet*qith(ind2p1)*rho(ind2p1) ! Deposition term AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.) BB=1./(rho(ind2p1)*Dv*qsith(ind2p1)) unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1))/qsith(ind2p1)*4.*RPI/(AA+BB)*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33) dep_term=rho(ind2p1)*fraca(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep ! Riming term unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4) rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim ! Precip term ! We assume that there is no solid precipitation outside thermals ! so the precipitation flux within thermals is equal to the precipitation flux ! at mesh-scale divided by thermals fraction fp2=0. fp1=0. IF (fraca(ind2p1) .GT. 0.) THEN fp2=-snowf(ind2p2)/fraca(ind2p1) ! flux defined positive upward fp1=-snowf(ind2p1)/fraca(ind2p1) ENDIF precip_term=-1./deltazlev(ind2p1)*(fp2-fp1) ! Calculation in a top-to-bottom loop IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN qi= 1./fm_therm(ind1,ind2p1)* & (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + & fm_therm(ind1,ind2p2)*(qith(ind2p1))) END IF ENDIF ! top thermals qith(ind2)=MAX(0.,qi) RETURN END SUBROUTINE ICE_MPC_BL_CLOUDS END MODULE cloudth_mod