Index: LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 (revision 1967) +++ LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 (revision 1968) @@ -9,8 +9,10 @@ !-------------------------------------------------------------------------- -!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance +! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance +! Last modified : Arnaud Jam 2014/02/11 +! Better representation of stratocumulus !-------------------------------------------------------------------------- - IMPLICIT NONE + IMPLICIT NONE #include "YOMCST.h" @@ -38,5 +40,4 @@ integer lev_out ! niveau pour les print integer nbpb - real zcon2(ngrid) real alim_star_tot(ngrid) @@ -62,10 +63,13 @@ REAL ztva_est(ngrid,klev) + REAL ztv_est(ngrid,klev) REAL zqla_est(ngrid,klev) REAL zqsatth(ngrid,klev) REAL zta_est(ngrid,klev) - REAL zdw2 + REAL ztemp(ngrid),zqsat(ngrid) + REAL zdw2,zdw2bis REAL zw2modif - REAL zeps + REAL zw2fact,zw2factbis + REAL zeps(ngrid,klev) REAL linter(ngrid) @@ -74,41 +78,41 @@ REAL wmaxa(ngrid) - INTEGER ig,l,k - - real zdz,zfact,zbuoy,zalpha,zdrag + INTEGER ig,l,k,lt,it + + real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m + real zbuoyjam(ngrid,klev) + real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis real zcor,zdelta,zcvm5,qlbef - real Tbef,qsatbef - real dqsat_dT,DT,num,denom + real betalpha,zbetalpha + real eps, afact REAL REPS,RLvCp,DDT0 PARAMETER (DDT0=.01) logical Zsat LOGICAL active(ngrid),activetmp(ngrid) - REAL fact_gamma,fact_epsilon,fact_gamma2 + REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2 REAL c2(ngrid,klev) - REAL a1,m - - REAL zw2fact,expa Zsat=.false. ! Initialisation + RLvCp = RLVTT/RCPD - - - fact_epsilon=0.002 - a1=2./3. - fact_gamma=0.9 - zfact=fact_gamma/(1+fact_gamma) - fact_gamma2=zfact - expa=0. - - -! Initialisations des variables reeles + fact_epsilon=0.002 + betalpha=0.9 + afact=2./3. + + zbetalpha=betalpha/(1.+betalpha) + + +! Initialisations des variables r?elles if (1==1) then ztva(:,:)=ztv(:,:) ztva_est(:,:)=ztva(:,:) + ztv_est(:,:)=ztv(:,:) ztla(:,:)=zthl(:,:) zqta(:,:)=po(:,:) + zqla(:,:)=0. zha(:,:) = ztva(:,:) else ztva(:,:)=0. + ztv_est(:,:)=0. ztva_est(:,:)=0. ztla(:,:)=0. @@ -128,10 +132,13 @@ entr(:,:)=0. zw2(:,:)=0. + zbuoy(:,:)=0. + zbuoyjam(:,:)=0. + gamma(:,:)=0. + zeps(:,:)=0. w_est(:,:)=0. f_star(:,:)=0. wa_moy(:,:)=0. linter(:)=1. - linter(:)=1. - +! linter(:)=1. ! Initialisation des variables entieres lmix(:)=1 @@ -139,4 +146,5 @@ wmaxa(:)=0. lalim(:)=1 + !------------------------------------------------------------------------- @@ -156,4 +164,5 @@ lalim(ig)=l+1 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) + print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l) endif enddo @@ -169,9 +178,11 @@ + + !------------------------------------------------------------------------------ ! Calcul dans la premiere couche ! 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 !------------------------------------------------------------------------------ @@ -210,28 +221,20 @@ -! Premier calcul de la vitesse verticale a partir de la temperature -! potentielle virtuelle -! if (1.eq.1) then -! w_est(ig,3)=zw2(ig,2)* & -! & ((f_star(ig,2))**2) & -! & /(f_star(ig,2)+alim_star(ig,2))**2+ & -! & 2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) & -! & *(zlev(ig,3)-zlev(ig,2)) -! endif - - !--------------------------------------------------------------------------- ! calcul des proprietes thermodynamiques et de la vitesse de la couche l ! 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)) - - do ig=1,ngrid - if(active(ig)) then + 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) @@ -239,30 +242,33 @@ ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & & -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) & - & /(f_star(ig,l)+alim_star(ig,l))**2+ & - & 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 - + + +!------------------------------------------------ +!AJAM:nouveau calcul de w? +!------------------------------------------------ + zdz=zlev(ig,l+1)-zlev(ig,l) + zdzbis=zlev(ig,l+1)-zlev(ig,l-1) + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + + zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) + zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) + zdw2=afact*zbuoy(ig,l)/fact_epsilon + zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon +! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) +! w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* & +! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & +! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) +! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact)) + w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & + & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & + & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) if (w_est(ig,l+1).lt.0.) then - w_est(ig,l+1)=zw2(ig,l) +! w_est(ig,l+1)=zw2(ig,l) + w_est(ig,l+1)=0.0001 endif endif enddo + !------------------------------------------------- !calcul des taux d'entrainement et de detrainement @@ -271,37 +277,67 @@ do ig=1,ngrid if (active(ig)) then + +! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) + zw2m=w_est(ig,l+1) +! zw2m=zw2(ig,l) zdz=zlev(ig,l+1)-zlev(ig,l) - zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) - -! 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 -! Calcul du taux d'entrainement entr_star (epsilon) - entr_star(ig,l)=f_star(ig,l)*zdz * ( zfact * MAX(0., & - & 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 * ( & -! & 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 * ( & -! & 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. )) - - + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) +! zbuoybis=zbuoy(ig,l)+RG*0.1/300. + zbuoybis=zbuoy(ig,l) + 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) + + +! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & +! & afact*zbuoybis/zw2m - fact_epsilon ) + +! entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha* & +! & afact*zbuoybis/zw2m - fact_epsilon ) + +!Modif AJAM + + lmel=0.1*zlev(ig,l) +! lmel=2.5*(zlev(ig,l)-zlev(ig,l-1)) + lt=l+1 + do it=1,klev-(l+1) + zdz2=zlev(ig,lt)-zlev(ig,l) + if (zdz2.gt.lmel) then + zdz3=zlev(ig,lt)-zlev(ig,lt-1) +! ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l) +! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) + + 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)- & + & ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) + +! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & +! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- & +! & po(ig,lt-1))/po(ig,lt-1)) + + else + lt=lt+1 + endif + enddo + +! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + + entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & + & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ) + + entrbis=entr_star(ig,l) + + + detr_star(ig,l)=f_star(ig,l)*zdz & + & *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & + & + 0.012*(zdqt(ig,l)/zw2m)**0.5 ) + + +! zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) +! +! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* & +! & afact*zbuoy(ig,l)/zw2m & +! & - 1.*fact_epsilon) + + ! En dessous de lalim, on prend le max de alim_star et entr_star pour ! alim_star et 0 sinon @@ -310,9 +346,9 @@ entr_star(ig,l)=0. 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) +! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then +! alim_star(ig,l)=entrbis ! endif + + print*,'alim0',l,lalim(ig),alim_star(ig,l),entrbis,f_star(ig,l) ! Calcul du flux montant normalise f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & @@ -321,4 +357,5 @@ endif enddo + !---------------------------------------------------------------------------- @@ -339,11 +376,11 @@ 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 if (activetmp(ig)) then ! 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) @@ -352,33 +389,17 @@ ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & & -zqla(ig,l))-zqla(ig,l)) - -!on ecrit zqsat - zqsatth(ig,l)=qsatbef - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! zw2(ig,l+1)=& -! & zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) & -! & +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) & -! & *1./(1.+fact_gamma) & -! & *(zlev(ig,l+1)-zlev(ig,l)) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! La meme en plus modulaire : - zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) + zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) zdz=zlev(ig,l+1)-zlev(ig,l) - - - zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) - -!if (1==0) then -! zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz) -! zdw2=2.*zbuoy/(1.+fact_gamma)*zdz -! zw2(ig,l+1)=zw2modif+zdw2 -!else - 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 -!endif - + zdzbis=zlev(ig,l+1)-zlev(ig,l-1) + zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) + + zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) + zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) + zdw2= afact*zbuoy(ig,l)/(fact_epsilon) + zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) +! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) + zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & + & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & + & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) endif enddo @@ -441,7 +462,19 @@ if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l - - return - end +#undef wrgrads_thermcell +#ifdef wrgrads_thermcell + call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ') + call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ') + call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ') + call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ') + call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ') + call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ') + call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ') +#endif + + + return + end + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -625,5 +658,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 !------------------------------------------------------------------------------ @@ -686,5 +719,5 @@ !------------------------------------------------ -!AJAM:nouveau calcul de w² +!AJAM:nouveau calcul de w? !------------------------------------------------ zdz=zlev(ig,l+1)-zlev(ig,l) @@ -855,3 +888,2 @@ return end -