Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 (revision 1337) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90 (revision 1338) @@ -130,4 +130,5 @@ real wmax(klon) + real wmax_tmp(klon) real wmax_sec(klon) real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev) @@ -191,6 +192,5 @@ -! #define wrgrads_thermcell -#undef wrgrads_thermcell +#define wrgrads_thermcell #ifdef wrgrads_thermcell ! Initialisation des sorties grads pour les thermiques. @@ -208,4 +208,6 @@ fm=0. ; entr=0. ; detr=0. + + print*,'THERMCELL MAIN OPT7' icount=icount+1 @@ -395,6 +397,5 @@ elseif (iflag_thermals_ed<=19) then -! Version d'Arnaud Jam -! print*,'THERM RIO et al 2010' +! print*,'THERM RIO et al 2010, version d Arnaud' CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & @@ -426,4 +427,11 @@ CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, & & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out) +! Attention, w2 est transforme en sa racine carree dans cette routine +! Le probleme vient du fait que linter et lmix sont souvent égaux à 1. + wmax_tmp=0. + do l=1,nlay + wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l)) + enddo +! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax @@ -665,56 +673,25 @@ enddo enddo -!calcul de ale_bl et alp_bl -!pour le calcul d'une valeur integree entre la surface et lmax - do ig=1,ngrid - alp_int(ig)=0. - ale_int(ig)=0. - n_int(ig)=0 - enddo -! - do l=1,nlay - do ig=1,ngrid - if(l.LE.lmax(ig)) THEN - alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l) - ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2 - n_int(ig)=n_int(ig)+1 - endif - enddo - enddo +! ! print*,'avant calcul ale et alp' !calcul de ALE et ALP pour la convection + Alp_bl(:)=0. + Ale_bl(:)=0. +! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)' + do l=1,nlay do ig=1,ngrid -! Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig)) -! Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig)) -! Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig)) -! & *0.1 -!valeur integree de alp_bl * 0.5: - if (n_int(ig).gt.0) then - Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig) -! if (Alp_bl(ig).lt.0.) then -! Alp_bl(ig)=0. - endif -! endif -! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)), -! s wth3(ig,nivcon(ig)),Alp_bl(ig) -! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig) -! Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2 -! if (nivcon(ig).eq.1) then -! Ale_bl(ig)=0. -! else -!valeur max de ale_bl: - Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2 -! & /2. -! & *0.1 -! Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2 -! if (n_int(ig).gt.0) then -! Ale_bl(ig)=ale_int(ig)/n_int(ig) -! Ale_bl(ig)=4. -! endif -! endif -! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig)) -! Ale_bl(ig)=wth2(ig,nivcon(ig)) -! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig) - enddo + Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) ) + Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2) +! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig) + enddo + enddo + +! print*,'AAAAAAA ',Alp_bl,Ale_bl,lmix + + +! TEST. IL FAUT REECRIRE LES ALE et ALP +! Ale_bl(:)=0.5*wmax(:)*wmax(:) +! Alp_bl(:)=0.1*wmax(:)*wmax(:)*wmax(:) + !test:calcul de la ponderation des couches pour KE !initialisations Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 (revision 1337) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 (revision 1338) @@ -64,4 +64,5 @@ REAL zqsatth(ngrid,klev) REAL zta_est(ngrid,klev) + REAL ztemp(ngrid),zqsat(ngrid) REAL zdw2 REAL zw2modif @@ -77,6 +78,4 @@ real zdz,zfact,zbuoy,zalpha real zcor,zdelta,zcvm5,qlbef - real Tbef,qsatbef - real dqsat_dT,DT,num,denom REAL REPS,RLvCp,DDT0 PARAMETER (DDT0=.01) @@ -116,4 +115,5 @@ ! write(lunout,*)'THERM 31H ' + print*,'THERMCELL_PLUME OPTIMISE V0 ' ! Initialisations des variables reeles if (1==0) then @@ -239,13 +239,17 @@ ! 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 !--------------------------------------------------------------------------- - call thermcell_condens(ngrid,active, & -& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l)) + ztemp(:)=zpspsk(:,l)*ztla(:,l-1) + call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat) + 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) @@ -366,6 +370,6 @@ enddo - call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l)) - + ztemp(:)=zpspsk(:,l)*ztla(:,l) + call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) do ig=1,ngrid @@ -373,4 +377,5 @@ 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) @@ -382,5 +387,4 @@ !on ecrit zqsat - zqsatth(ig,l)=qsatbef !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -554,4 +558,5 @@ REAL zqsatth(ngrid,klev) REAL zta_est(ngrid,klev) + REAL ztemp(ngrid),zqsat(ngrid) REAL zdw2 REAL zw2modif @@ -570,6 +575,5 @@ real zcor,zdelta,zcvm5,qlbef,zdz2 real betalpha,zbetalpha - real Tbef,qsatbef,b1,eps, afact - real dqsat_dT,DT,num,denom,m + real eps, afact REAL REPS,RLvCp,DDT0 PARAMETER (DDT0=.01) @@ -589,4 +593,5 @@ ! print*,'THERM 31B' + print*,'THERMCELL_PLUME OPTIMISE V1 CCC ' ! Initialisations des variables reeles @@ -624,5 +629,4 @@ linter(:)=1. ! linter(:)=1. - b1=2. ! Initialisation des variables entieres lmix(:)=1 @@ -631,5 +635,5 @@ lalim(:)=1 -! print*,'THERMCELL PLUME ARNAUD DEDANS' + print*,'THERMCELL PLUME QSAT2 NDDDDN' !------------------------------------------------------------------------- @@ -708,15 +712,17 @@ ! 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 !--------------------------------------------------------------------------- - call thermcell_condens(ngrid,active, & -& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l)) - + ztemp(:)=zpspsk(:,l)*ztla(:,l-1) + call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) do ig=1,ngrid ! print*,'active',active(ig),ig,l 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) @@ -747,4 +753,5 @@ !------------------------------------------------- + print*,'THERM V1 SANS DQ' do ig=1,ngrid if (active(ig)) then @@ -758,4 +765,5 @@ zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) + zdqt(ig,l)=0. @@ -800,6 +808,6 @@ enddo - call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l)) - + ztemp(:)=zpspsk(:,l)*ztla(:,l) + call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) do ig=1,ngrid @@ -808,4 +816,5 @@ ! on ecrit de maniere conservative (sat ou non) ! T = Tl +Lv/Cp ql + zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l)) ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) @@ -814,8 +823,4 @@ ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & & -zqla(ig,l))-zqla(ig,l)) - -!on ecrit zqsat - zqsatth(ig,l)=qsatbef - zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) zdz=zlev(ig,l+1)-zlev(ig,l) Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_qsat.F90 =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_qsat.F90 (revision 1337) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_qsat.F90 (revision 1338) @@ -1,3 +1,3 @@ -subroutine thermcell_qsat(klon,mask,pplev,ptemp,pqt,pqsat) +subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat) implicit none @@ -5,4 +5,5 @@ #include "YOETHF.h" #include "FCTTRE.h" + !==================================================================== @@ -12,12 +13,12 @@ ! Arguments INTEGER klon -REAL pplev(klon) -REAL ptemp(klon),pqt(klon),pqsat(klon) -LOGICAL mask(klon) +REAL zpspsk(klon),pplev(klon) +REAL ztemp(klon),zqta(klon),zqsat(klon) +LOGICAL active(klon) ! Variables locales INTEGER ig,iter REAL Tbef(klon),DT(klon) -REAL tdelta,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT +REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT logical Zsat REAL RLvCp @@ -31,50 +32,56 @@ RLvCp = RLVTT/RCPD tout_converge=.false. -afaire(:)=mask(:) +afaire(:)=.false. DT(:)=0. !==================================================================== -! Routine a vectoriser en copiant mask dans converge et en mettant +! Routine a vectoriser en copiant active dans converge et en mettant ! la boucle sur les iterations a l'exterieur est en mettant ! converge= false des que la convergence est atteinte. -! Pr Tprec=Tl calcul de qsat -! Si qsat>qT T=Tl, q=qT -! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt) -! On cherche DDT < DDT0 !==================================================================== do ig=1,klon - if (mask(ig)) then - Tbef(ig)=ptemp(ig) + if (active(ig)) then + Tbef(ig)=ztemp(ig) zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig))) - pqsat(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig) - pqsat(ig)=MIN(0.5,pqsat(ig)) - zcor=1./(1.-retv*pqsat(ig)) - pqsat(ig)=pqsat(ig)*zcor - qlbef=max(0.,pqt(ig)-pqsat(ig)) - afaire(ig)=qlbef>1.e-10 + qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig) + qsatbef=MIN(0.5,qsatbef) + zcor=1./(1.-retv*qsatbef) + qsatbef=qsatbef*zcor + qlbef=max(0.,zqta(ig)-qsatbef) DT(ig) = 0.5*RLvCp*qlbef + zqsat(ig)=qsatbef endif enddo +! Traitement du cas ou il y a condensation mais faible +! On ne condense pas mais on dit que le qsat est le qta +do ig=1,klon + if (active(ig)) then + if (0.DDT0).and.afaire(:) + afaire(:)=abs(DT(:)).gt.DDT0 do ig=1,klon if (afaire(ig)) then Tbef(ig)=Tbef(ig)+DT(ig) zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig))) - pqsat(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig) - pqsat(ig)=MIN(0.5,pqsat(ig)) - zcor=1./(1.-retv*pqsat(ig)) - pqsat(ig)=pqsat(ig)*zcor - qlbef=pqt(ig)-pqsat(ig) + qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig) + qsatbef=MIN(0.5,qsatbef) + zcor=1./(1.-retv*qsatbef) + qsatbef=qsatbef*zcor + qlbef=zqta(ig)-qsatbef zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig))) zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta - zcor=1./(1.-retv*pqsat(ig)) - dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,pqsat(ig),zcor) - num=-Tbef(ig)+ptemp(ig)+RLvCp*qlbef + zcor=1./(1.-retv*qsatbef) + dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor) + num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef denom=1.+RLvCp*dqsat_dT + zqsat(ig) = qsatbef DT(ig)=num/denom endif