Index: /LMDZ5/trunk/libf/phylmd/cloudth.F90 =================================================================== --- /LMDZ5/trunk/libf/phylmd/cloudth.F90 (revision 2266) +++ /LMDZ5/trunk/libf/phylmd/cloudth.F90 (revision 2267) @@ -5,4 +5,5 @@ + USE IOIPSL, ONLY : getin IMPLICIT NONE @@ -39,5 +40,5 @@ - REAL sigma1(ngrid,klev) + REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) @@ -48,5 +49,5 @@ REAL ctot(ngrid,klev) REAL rneb(ngrid,klev) - REAL t(ngrid,klev) + REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi REAL rdd,cppd,Lv @@ -62,18 +63,39 @@ REAL erf - - - - -! print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), & -! & fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), & -! & 'verif' - - -! LOGICAL active(ngrid) - -!----------------------------------------------------------------------------------------------------------------- + REAL, SAVE :: iflag_cloudth_vert, iflag_cloudth_vert_omp=0 + + + LOGICAL, SAVE :: first=.true. + + + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Astuce pour gérer deux versions de cloudth en attendant +! de converger sur une version nouvelle +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IF (first) THEN + !$OMP MASTER + CALL getin('iflag_cloudth_vert',iflag_cloudth_vert_omp) + !$OMP END MASTER + !$OMP BARRIER + iflag_cloudth_vert=iflag_cloudth_vert_omp + first=.false. + ENDIF + IF (iflag_cloudth_vert==1) THEN + CALL cloudth_vert(ngrid,klev,ind2, & + & ztv,po,zqta,fraca, & + & qcloud,ctot,zpspsk,paprs,ztla,zthl, & + & ratqs,zqs,t) + RETURN + ENDIF +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + +!------------------------------------------------------------------------------- ! Initialisation des variables r?elles -!----------------------------------------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- sigma1(:,:)=0. sigma2(:,:)=0. @@ -96,7 +118,7 @@ -!------------------------------------------------------------------------------------------------------------------ +!------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions -!------------------------------------------------------------------------------------------------------------------ +!------------------------------------------------------------------------------- do ind1=1,ngrid @@ -139,7 +161,7 @@ -!----------------------------------------------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! 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) @@ -155,7 +177,7 @@ ! 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) @@ -176,5 +198,5 @@ ! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha' -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (ctot(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. @@ -242,9 +264,285 @@ - - - - - - - +!=========================================================================== + SUBROUTINE cloudth_vert(ngrid,klev,ind2, & + & ztv,po,zqta,fraca, & + & qcloud,ctot,zpspsk,paprs,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 "iniprint.h" +#include "thermcell.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 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 alpha,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 + + + + + +!------------------------------------------------------------------------------ +! 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) + + + +!------------------------------------------------------------------------------- +! 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) +! ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*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) +! qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2) + + +! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha' + + +!------------------------------------------------------------------------------- +! 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) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + 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 Index: /LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 =================================================================== --- /LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 (revision 2266) +++ /LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 (revision 2267) @@ -6,6 +6,6 @@ & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & - & ,lev_out,lunout1,igout) -! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) + & ,lev_out,lunout1,igout) +! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) !-------------------------------------------------------------------------- !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance @@ -87,5 +87,5 @@ real ztv_est1,ztv_est2 real zcor,zdelta,zcvm5,qlbef - real zbetalpha + real zbetalpha, coefzlmel real eps REAL REPS,RLvCp,DDT0 @@ -396,4 +396,5 @@ else ! if (iflag_thermals_ed.lt.8) then lt=l+1 + zlt=zlev(ig,lt) zdz2=zlev(ig,lt)-zlev(ig,l) @@ -405,7 +406,8 @@ zdz3=zlev(ig,lt+1)-zlt zltdwn=zlev(ig,lt)-zdz3/2 - - zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- & - & ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- & + zlmelup=zlmel+(zdz/2) + coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz) + zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & + & ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- & & ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) endif ! if (iflag_thermals_ed.lt.8) then @@ -422,4 +424,5 @@ zdw2=afact*zbuoy(ig,l)/fact_epsilon zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon +! zdw2bis=0.5*(zdw2+zdw2bis) lm=Max(1,l-2) ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & @@ -442,5 +445,5 @@ ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) - w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) + w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) ! Nouvelle version Arnaud @@ -556,5 +559,5 @@ & mix0 * 0.1 / (zalpha+0.001) & & + zbetalpha*MAX(entr_min, & - & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )) + & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) @@ -645,7 +648,11 @@ ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) + if (iflag_thermals_ed==8) then + zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) + else zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) + endif ! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* & -! & (zw2(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdz+zdzbis))* & +! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* & ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))