Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90 (revision 1398) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90 (revision 1399) @@ -8,5 +8,5 @@ & ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth, & & ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, & - & zmax0,f0,zw2,fraca) + & zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl) USE dimphy @@ -49,4 +49,8 @@ real zqla(klon,klev) real zqta(klon,klev) + real ztv(klon,klev) + real zpspsk(klon,klev) + real ztla(klon,klev) + real zthl(klon,klev) real wmax_sec(klon) real zmax_sec(klon) @@ -214,5 +218,6 @@ & ,r_aspect_thermals,l_mix_thermals & & ,tau_thermals,iflag_thermals_ed,Ale,Alp,lalim_conv,wght_th & - & ,zmax0,f0,zw2,fraca) + & ,zmax0,f0,zw2,fraca,ztv,zpspsk & + & ,ztla,zthl) if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK' else Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90 (revision 1399) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90 (revision 1399) @@ -0,0 +1,263 @@ + SUBROUTINE cloudth(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 + REAL alth,alenv,ath,aenv + REAL sth,senv,sigma1s,sigma2s,xth,xenv + 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 + + + + + +! 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) + +!----------------------------------------------------------------------------------------------------------------- +! 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 +!----------------------------------------------------------------------------------------------------------------- +! alpha=0.5*(fraca(ind1,ind2)+fraca(ind1,ind2-1)) +! sigma1s=(Max(1.2*fraca(ind1,ind2)*(sth-senv)**2,(senv/100)**2))**0.5 +! sigma2s=((0.021/(fraca(ind1,ind2)+0.0)**0.5)*(sth-senv)**2+(sth/100)**2)**0.5 +! sigma1s=(1.5**0.5)*(fraca(ind1,ind2)**0.65)*(sth-senv)+0.000045 +! sigma2s=0.1265*(sth-senv)/(fraca(ind1,ind2)+0.0)**0.35+0.000045 + + sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.00003 + sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 + +! sigma1s=(1.5**0.5)*(alpha**0.65)*(sth-senv)+0.000045 +! sigma2s=0.1265*(sth-senv)/alpha**0.35+0.000045 + +! sigma1s=(2.8**0.5)*(0.1**0.7)*(sth-senv)+0.00002 +! sigma2s=((0.126/(0.1)**0.3)*(sth-senv))+0.00002 + + + +!----------------------------------------------------------------------------------------------------------------- +! 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' + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + 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)) + +! if (zqenv(ind1).gt.0.005) then + sigma1s=0.005*zqenv(ind1) +! else +! sigma1s=0.0005 +! endif + +! sigma1s=ratqs(ind1,ind2)*zqenv(ind1) + + + +! sigma1s=0.00003 + 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: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F (revision 1398) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F (revision 1399) @@ -7,5 +7,6 @@ s pfrac_impa, pfrac_nucl, pfrac_1nucl, s frac_impa, frac_nucl, - s prfl, psfl, rhcl) + s prfl, psfl, rhcl, zqta, fraca, + s ztv, zpspsk, ztla, zthl, iflag_cldcon) c @@ -41,5 +42,12 @@ REAL snow(klon) ! neige (mm/s) REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) - REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) + REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) + REAL ztv(klon,klev) + REAL zqta(klon,klev),fraca(klon,klev) + REAL sigma1(klon,klev),sigma2(klon,klev) + REAL qltot(klon,klev),ctot(klon,klev) + REAL zpspsk(klon,klev),ztla(klon,klev) + REAL zthl(klon,klev) + cAA c Coeffients de fraction lessivee : pour OFF-LINE @@ -63,5 +71,6 @@ INTEGER ninter ! sous-intervals pour la precipitation - INTEGER ncoreczq + INTEGER ncoreczq + INTEGER iflag_cldcon PARAMETER (ninter=5) LOGICAL evap_prec ! evaporation de la pluie @@ -72,5 +81,6 @@ real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) - real erf + real erf + REAL qcloud(klon) c LOGICAL cpartiel ! condensation partielle @@ -82,5 +92,5 @@ c INTEGER i, k, n, kk - REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 + REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 REAL zrfl(klon), zrfln(klon), zqev, zqevt REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq @@ -130,4 +140,5 @@ zdelq=0.0 + print*,'CLOUDTH4 A. JAM' IF (appel1er) THEN c @@ -322,4 +333,5 @@ c de l'eau condensee: c + IF (cpartiel) THEN @@ -351,6 +363,21 @@ zq(i)=1.e-15 endif - enddo - do i=1,klon + enddo + + if (iflag_cldcon.eq.5) then + + call cloudth(klon,klev,k,ztv, + . zq,zqta,fraca, + . qcloud,ctot,zpspsk,paprs,ztla,zthl, + . ratqs,zqs,t) + + do i=1,klon + rneb(i,k)=ctot(i,k) + zqn(i)=qcloud(i) + enddo + + else + + do i=1,klon zpdf_sig(i)=ratqs(i,k)*zq(i) zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) @@ -372,5 +399,7 @@ endif - enddo + enddo + + endif ! iflag_cldcon endif ! iflag_pdf Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F (revision 1398) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F (revision 1399) @@ -669,5 +669,9 @@ REAL zw2(klon,klev+1) - REAL fraca(klon,klev+1) + REAL fraca(klon,klev+1) + REAL ztv(klon,klev) + REAL zpspsk(klon,klev) + REAL ztla(klon,klev) + REAL zthl(klon,klev) c Variables locales pour la couche limite (al1): @@ -2440,5 +2444,6 @@ s ,ratqsdiff,zqsatth con rajoute ale et alp, et les caracteristiques de la couche alim - s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca) + s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca + s ,ztv,zpspsk,ztla,zthl) ! ---------------------------------------------------------------------- @@ -2702,5 +2707,6 @@ . pfrac_impa, pfrac_nucl, pfrac_1nucl, . frac_impa, frac_nucl, - . prfl, psfl, rhcl) + . prfl, psfl, rhcl, + . zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon ) WHERE (rain_lsc < 0) rain_lsc = 0. Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 (revision 1398) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 (revision 1399) @@ -10,5 +10,6 @@ & ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed & & ,Ale_bl,Alp_bl,lalim_conv,wght_th & - & ,zmax0, f0,zw2,fraca) + & ,zmax0, f0,zw2,fraca,ztv & + & ,zpspsk,ztla,zthl) USE dimphy @@ -389,5 +390,5 @@ ! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed if (iflag_thermals_ed<=9) then -! print*,'THERM NOVUELLE/NOUVELLE/ANCIENNE' +! print*,'THERM NOUVELLE/NOUVELLE Arnaud' CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & @@ -396,5 +397,5 @@ & ,lev_out,lunout1,igout) - elseif (iflag_thermals_ed<=19) then + elseif (iflag_thermals_ed>9) then ! print*,'THERM RIO et al 2010, version d Arnaud' CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 (revision 1398) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 (revision 1399) @@ -64,5 +64,4 @@ REAL zqsatth(ngrid,klev) REAL zta_est(ngrid,klev) - REAL ztemp(ngrid),zqsat(ngrid) REAL zdw2 REAL zw2modif @@ -76,6 +75,8 @@ INTEGER ig,l,k - real zdz,zfact,zbuoy,zalpha + real zdz,zfact,zbuoy,zalpha,zdrag real zcor,zdelta,zcvm5,qlbef + real Tbef,qsatbef + real dqsat_dT,DT,num,denom REAL REPS,RLvCp,DDT0 PARAMETER (DDT0=.01) @@ -84,4 +85,5 @@ REAL fact_gamma,fact_epsilon,fact_gamma2 REAL c2(ngrid,klev) + REAL a1,m REAL zw2fact,expa @@ -90,32 +92,16 @@ RLvCp = RLVTT/RCPD - if (iflag_thermals_ed==0) then - fact_gamma=1. - fact_epsilon=1. - else if (iflag_thermals_ed==1) then -! Valeurs au moment de la premiere soumission des papiers - fact_gamma=1. + fact_epsilon=0.002 - fact_gamma2=0.6 -! Suggestions des Fleurs, Septembre 2009 - fact_epsilon=0.015 -!test cr -! fact_epsilon=0.002 + a1=2./3. fact_gamma=0.9 - fact_gamma2=0.7 - - else if (iflag_thermals_ed==2) then - fact_gamma=1. - fact_epsilon=2. - else if (iflag_thermals_ed==3) then - fact_gamma=3./4. - fact_epsilon=3. - endif - -! write(lunout,*)'THERM 31H ' - - print*,'THERMCELL_PLUME OPTIMISE V0 ' + zfact=fact_gamma/(1+fact_gamma) + fact_gamma2=zfact + expa=0. + + print*,'THERM 31H ' + ! Initialisations des variables reeles -if (1==0) then +if (1==1) then ztva(:,:)=ztv(:,:) ztva_est(:,:)=ztva(:,:) @@ -168,5 +154,5 @@ alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & & *sqrt(zlev(ig,l+1)) - lalim(ig)=l+1 + lalim(:)=l+1 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) endif @@ -187,5 +173,5 @@ ! On decide dans cette version que le thermique n'est actif que si la premiere ! couche est instable. -! Pourrait etre change si on veut que le thermiques puisse se déclencher +! Pourrait etre change si on veut que le thermiques puisse se déclencher ! dans une couche l>1 !------------------------------------------------------------------------------ @@ -239,17 +225,13 @@ ! sans tenir compte du detrainement et de l'entrainement dans cette ! couche -! C'est a dire qu'on suppose -! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer ! avant) a l'alimentation pour avoir un calcul plus propre !--------------------------------------------------------------------------- - ztemp(:)=zpspsk(:,l)*ztla(:,l-1) - call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat) - + call thermcell_condens(ngrid,active, & +& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l)) do ig=1,ngrid if(active(ig)) then - zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig)) ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) zta_est(ig,l)=ztva_est(ig,l) @@ -258,4 +240,6 @@ & -zqla_est(ig,l))-zqla_est(ig,l)) + if (1.eq.0) then +!calcul de w_est sans prendre en compte le drag w_est(ig,l+1)=zw2(ig,l)* & & ((f_star(ig,l))**2) & @@ -263,4 +247,16 @@ & 2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) & & *(zlev(ig,l+1)-zlev(ig,l)) + else + + zdz=zlev(ig,l+1)-zlev(ig,l) + zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l) + zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + zdrag=fact_epsilon/(zalpha**expa) + zw2fact=zbuoy/zdrag*a1 + w_est(ig,l+1)=(w_est(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) & + & +zw2fact + + endif + if (w_est(ig,l+1).lt.0.) then w_est(ig,l+1)=zw2(ig,l) @@ -277,65 +273,34 @@ zdz=zlev(ig,l+1)-zlev(ig,l) zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) - zfact=fact_gamma/(1.+fact_gamma) ! estimation de la fraction couverte par les thermiques zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l) -!calcul de la soumission papier - if (1.eq.1) then - fact_epsilon=0.0007 -! zfact=0.9/(1.+0.9) - zfact=0.3 - fact_gamma=0.7 - fact_gamma2=0.6 - expa=0.25 +!calcul de la soumission papier ! Calcul du taux d'entrainement entr_star (epsilon) entr_star(ig,l)=f_star(ig,l)*zdz * ( zfact * MAX(0., & - & zbuoy/w_est(ig,l+1) )& -!- fact_epsilon/zalpha**0.25 ) & - & +0.000 ) - -! entr_star(ig,l)=f_star(ig,l)*zdz * ( 1./3 * MAX(0., & -! & zbuoy/w_est(ig,l+1) - 1./zalpha**0.25 ) & -! & +0.000 ) -! Calcul du taux de detrainement detr_star (delta) -! if (zqla_est(ig,l).lt.1.e-10) then + & a1*zbuoy/w_est(ig,l+1) & + & - fact_epsilon/zalpha**expa ) & + & +0. ) + +!calcul du taux de detrainment (delta) ! detr_star(ig,l)=f_star(ig,l)*zdz * ( & -! & fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) ) & -! & +0.0006 ) -! else +! & MAX(1.e-3, & +! & -fact_gamma2*a1*zbuoy/w_est(ig,l+1) & +! & +0.01*(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))**0.5 & +! & +0. )) + + m=0.5 + + detr_star(ig,l)=1.*f_star(ig,l)*zdz * & + & MAX(5.e-4,-fact_gamma2*a1*(1./w_est(ig,l+1))*((1.-(1.-m)/(1.+70*zqta(ig,l-1)))*zbuoy & + & -40*(1.-m)*(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) ) ) + ! detr_star(ig,l)=f_star(ig,l)*zdz * ( & -! & fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) ) & -! & +0.002 ) -! endif - detr_star(ig,l)=f_star(ig,l)*zdz * ( & - & fact_gamma2 * MAX(0., & -!fact_epsilon/zalpha**0.25 - & -zbuoy/w_est(ig,l+1) ) & -! & + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6 & -!test - & + 0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6 & - & +0.0000 ) - else - -! Calcul du taux d'entrainement entr_star (epsilon) - entr_star(ig,l)=f_star(ig,l)*zdz * ( zfact * MAX(0., & - & zbuoy/w_est(ig,l+1) - fact_epsilon ) & - & +0.0000 ) - -! Calcul du taux de detrainement detr_star (delta) - detr_star(ig,l)=f_star(ig,l)*zdz * ( & - & fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) ) & - & + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6 & - & +0.0000 ) - - endif -!endif choix du calcul de E* et D* - -!cr test -! entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l) - -! Prise en compte de la fraction -! detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5)) +! & MAX(0.0, & +! & -fact_gamma2*a1*zbuoy/w_est(ig,l+1) & +! & +20*(max(zqta(ig,l-1)-po(ig,l),0.))**1*(zalpha/w_est(ig,l+1))**0.5 & +! & +0. )) + ! En dessous de lalim, on prend le max de alim_star et entr_star pour @@ -346,4 +311,8 @@ endif +!attention test +! if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then +! detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) +! endif ! Calcul du flux montant normalise f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & @@ -370,6 +339,6 @@ enddo - ztemp(:)=zpspsk(:,l)*ztla(:,l) - call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) + call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l)) + do ig=1,ngrid @@ -377,5 +346,4 @@ if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l ! on ecrit de maniere conservative (sat ou non) - zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l)) ! T = Tl +Lv/Cp ql ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) @@ -387,4 +355,5 @@ !on ecrit zqsat + zqsatth(ig,l)=qsatbef !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -407,25 +376,12 @@ zw2(ig,l+1)=zw2modif+zdw2 else -! Tentative de reecriture de l'equation de w2. A reprendre ... -! zdw2=2*zdz*zbuoy -! zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon)) -!!!!! zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l)) -!formulation Arnaud -! zw2fact=zbuoy*zalpha**expa/fact_epsilon -! zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) & -! & +zw2fact - if (zbuoy.gt.1.e-10) then - zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact) - zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) & + zdrag=fact_epsilon/(zalpha**expa) + zw2fact=zbuoy/zdrag*a1 + zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) & & +zw2fact - else - zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma) - zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) & - & +zw2fact - - endif + endif -! zw2(ig,l+1)=zw2modif+zdw2 + endif enddo @@ -440,6 +396,5 @@ if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then ! stop'On tombe sur le cas particulier de thermcell_dry' - write(lunout,*) & - & 'On tombe sur le cas particulier de thermcell_plume' + print*,'On tombe sur le cas particulier de thermcell_plume' zw2(ig,l+1)=0. linter(ig)=l+1 @@ -487,5 +442,4 @@ return end - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!