Index: /LMDZ6/trunk/libf/phylmd/lmdz_thermcell_main.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/lmdz_thermcell_main.F90 (revision 5894) +++ /LMDZ6/trunk/libf/phylmd/lmdz_thermcell_main.F90 (revision 5895) @@ -429,5 +429,5 @@ CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & - & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & + & lalim,f0,zmax0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) Index: /LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume.f90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume.f90 (revision 5894) +++ /LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume.f90 (revision 5895) @@ -1,5 +1,5 @@ MODULE lmdz_thermcell_plume ! -! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $ +! $Id: lmdz_thermcell_plume.f90 5854 2025-10-10 14:00:56Z fhourdin $ ! CONTAINS @@ -7,20 +7,10 @@ SUBROUTINE thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & - & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & + & lalim,f0,zmax0,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) !-------------------------------------------------------------------------- -! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam -! !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance -! This versions starts from a cleaning of thermcell_plume_6A (2019/01/20) -! thermcell_plume_6A is activate for flag_thermas_ed < 10 -! thermcell_plume_5B for flag_thermas_ed < 20 -! thermcell_plume for flag_thermals_ed>= 20 -! Various options are controled by the flag_thermals_ed parameter -! = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8 -! = 21 : the Jam strato-cumulus modif is not activated in detrainment -! = 29 : an other way to compute the modified buoyancy (to be tested) !-------------------------------------------------------------------------- @@ -46,5 +36,5 @@ real,intent(in),dimension(ngrid,nlay) :: pphi real,intent(in),dimension(ngrid,nlay) :: zpspsk - real,intent(in),dimension(ngrid) :: f0 + real,intent(in),dimension(ngrid) :: f0,zmax0 integer,intent(out) :: lalim(ngrid) @@ -64,15 +54,8 @@ real,intent(out),dimension(ngrid,nlay) :: ztva_est real,intent(out),dimension(ngrid,nlay) :: zqsatth - integer,intent(out),dimension(ngrid) :: lmix(ngrid) - integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid) - real,intent(out),dimension(ngrid) :: linter(ngrid) - - - REAL,dimension(ngrid,nlay+1) :: wa_moy - REAL,dimension(ngrid,nlay) :: entr,detr - REAL,dimension(ngrid,nlay) :: ztv_est - REAL,dimension(ngrid,nlay) :: zqla_est - REAL,dimension(ngrid,nlay) :: zta_est - REAL,dimension(ngrid) :: ztemp,zqsat + integer,intent(out),dimension(ngrid) :: lmix + integer,intent(out),dimension(ngrid) :: lmix_bis + real,intent(out),dimension(ngrid) :: linter + REAL zdw2,zdw2bis REAL zw2modif @@ -80,13 +63,22 @@ REAL,dimension(ngrid,nlay) :: zeps - REAL,dimension(ngrid) :: wmaxa - - INTEGER ig,l,k,lt,it,lm,nbpb + REAL, dimension(ngrid) :: wmaxa + + INTEGER ig,l,k,lt,it,lm + integer nbpb + + real,dimension(ngrid,nlay) :: detr + real,dimension(ngrid,nlay) :: entr + real,dimension(ngrid,nlay+1) :: wa_moy + real,dimension(ngrid,nlay) :: ztv_est + real,dimension(ngrid) :: ztemp,zqsat + real,dimension(ngrid,nlay) :: zqla_est + real,dimension(ngrid,nlay) :: zta_est real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt real zdz,zalpha,zw2m real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam - real zdz2,zdz3,lmel,entrbis,zdzbis - real,dimension(ngrid) :: d_temp + real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis + real, dimension(ngrid) :: d_temp real ztv1,ztv2,factinv,zinv,zlmel real zlmelup,zlmeldwn,zlt,zltdwn,zltup @@ -99,12 +91,26 @@ LOGICAL,dimension(ngrid) :: active,activetmp REAL fact_gamma,fact_gamma2,fact_epsilon2 - - + REAL coefc REAL,dimension(ngrid,nlay) :: c2 - if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11' - Zsat=.false. -! Initialisation - + real, dimension(ngrid,nlay) :: z_detr_q,z_cld_radius + real, dimension(ngrid) :: z_cld_radius_base,z_cld_base_height,z_alpha_base + + integer choice_ed_main,choice_ed_dq + + z_cld_base_height(:)=0. + + if (ngrid==1) print*,'THERMCELL PLUME Modifie 2025/11/11' + +! --------------------------------------------------------------- +! Controling entrainment and detrainment +! iflag_thermals_ed=1XY +! choice_ed_main=X +! choice_ed_dq=Y +! X=0 et Y=0 <=> thermcell_plume_6A + choice_ed_main=(iflag_thermals_ed-100)/10 + choice_ed_dq=(iflag_thermals_ed-100)-10*choice_ed_main +! --------------------------------------------------------------- + zbetalpha=betalpha/(1.+betalpha) @@ -112,4 +118,5 @@ ! Initialisations des variables r?elles +Zsat=.false. if (1==1) then ztva(:,:)=ztv(:,:) @@ -154,4 +161,11 @@ wmaxa(:)=0. +! Initialisation a 0 en cas de sortie dans replay + zqsat(:)=0. + zta_est(:,:)=0. + zdqt(:,:)=0. + zdqtjam(:,:)=0. + c2(:,:)=0. + !------------------------------------------------------------------------- @@ -221,11 +235,10 @@ do ig=1,ngrid ! print*,'active',active(ig),ig,l - if(active(ig)) then + 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) ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) - ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & - & -zqla_est(ig,l))-zqla_est(ig,l)) + ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)-zqla_est(ig,l))) @@ -250,5 +263,4 @@ !========================================================================= -!-------------------------------------------------- lt=l+1 zlt=zlev(ig,lt) @@ -256,7 +268,7 @@ do while (lmel.gt.zdz2) - lt=lt+1 - zlt=zlev(ig,lt) - zdz2=zlt-zlev(ig,l) + lt=lt+1 + zlt=zlev(ig,lt) + zdz2=zlt-zlev(ig,l) enddo zdz3=zlev(ig,lt+1)-zlt @@ -274,4 +286,5 @@ zdzbis=zlev(ig,l)-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) @@ -280,4 +293,6 @@ lm=Max(1,l-2) w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) + + endif enddo @@ -285,5 +300,5 @@ !------------------------------------------------- -!calcul des taux d'entrainement et de detrainement +! 4. calcul des taux d'entrainement et de detrainement !------------------------------------------------- @@ -293,29 +308,91 @@ ! 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(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) -!========================================================================= -! 4. Calcul de l'entrainement et du detrainement -!========================================================================= - - detr_star(ig,l)=f_star(ig,l)*zdz & - & *( mix0 * 0.1 / (zalpha+0.001) & - & + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & - & + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) - - if ( iflag_thermals_ed == 20 ) then - entr_star(ig,l)=f_star(ig,l)*zdz* ( & - & mix0 * 0.1 / (zalpha+0.001) & - & + zbetalpha*MAX(entr_min, & - & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) + + +! detr_q_coef = thermals_detr_q_coef lu dans les .def +! vrai flux de masse : f0*fstar +! equation de conservation de la masse s'écrit : f*(k+1) - f*(k) = e*(k) - d*(k) +! e=f0 e* / dz + + if ( choice_ed_dq == 0 ) then + + z_detr_q(ig,l)=detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power + + elseif ( choice_ed_dq == 1 ) then + + ! Tant que la couche du dessous n'est pas condensée, on ne tient pas + ! compte du détrainement en q + ! FH V1 : if ( zqla(ig,l-1) > 0. ) then + if ( zqla_est(ig,l) > 0. ) then + if ( z_cld_base_height(ig) == 0. ) then + ! Cloud base : height and fraction + z_cld_base_height(ig)=zlev(ig,l) ! z at cloud base + z_alpha_base(ig)=zalpha + endif + ! Cloud radius. At cloud base = cloud_height / 2. Cloud_height=zmax0-z_cld_base_height + ! With min value dz of the layer + ! FH V1 : z_cld_radius(ig,l)=sqrt(zalpha/z_alpha_base(ig))*0.5*(max(zmax0(ig)-z_cld_base_height(ig),zlev(ig,l+1)-zlev(ig,l))) + z_cld_radius(ig,l)=0.5*(max(zmax0(ig)-z_cld_base_height(ig),zlev(ig,l+1)-zlev(ig,l))) + z_detr_q(ig,l)=detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power*2./z_cld_radius(ig,l) + else + z_detr_q(ig,l)=0. + endif + else - entr_star(ig,l)=f_star(ig,l)*zdz* ( & - & mix0 * 0.1 / (zalpha+0.001) & - & + zbetalpha*MAX(entr_min, & - & afact*zbuoy(ig,l)/zw2m - fact_epsilon)) + + print*,'choice_ed_dq=',choice_ed_dq,' non prevu' + stop + endif + +! zbetalpha = B1 / ( 1 + B1 ) +! afact = A1 +! zbuoyjam = B +! zw2m = w2 +! z_detr_q = terme de détrainement de Delta q + + if ( choice_ed_main == 0 ) then + + detr_star(ig,l)=f_star(ig,l)*zdz & + & *( mix0 * 0.1 / (zalpha+0.001) & + & + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & + & + z_detr_q(ig,l)) ) + + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + entr_star(ig,l)=f_star(ig,l)*zdz* ( & + & mix0 * 0.1 / (zalpha+0.001) & + & + zbetalpha*MAX(entr_min, & + & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) + + elseif ( choice_ed_main == 1 ) then + + detr_star(ig,l)=f_star(ig,l)*zdz & + & *( mix0 * 0.1 / (zalpha+0.001) & + & + detr_min & + & + MAX(-afact*zbetalpha*zbuoyjam(ig,l)/zw2m,0.) & + & + z_detr_q(ig,l) ) + + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + entr_star(ig,l)=f_star(ig,l)*zdz* ( & + & mix0 * 0.1 / (zalpha+0.001) & + & + entr_min & + & + zbetalpha*MAX(0., & + & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) + + + else + print*,'choice_ed_dq=',choice_ed_dq,' non prevu' + stop + endif + + ! En dessous de lalim, on prend le max de alim_star et entr_star pour @@ -325,5 +402,11 @@ entr_star(ig,l)=0. endif - f_star(ig,l+1)=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',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),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) & & -detr_star(ig,l) @@ -361,6 +444,5 @@ !on rajoute le calcul de zha pour diagnostiques (temp potentielle) zha(ig,l) = ztva(ig,l) - ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & - & -zqla(ig,l))-zqla(ig,l)) + ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)-zqla(ig,l))) zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) zdz=zlev(ig,l+1)-zlev(ig,l) @@ -375,5 +457,5 @@ enddo - if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l + if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l ! !=========================================================================== @@ -384,6 +466,4 @@ do ig=1,ngrid 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' -! print*,'On tombe sur le cas particulier de thermcell_plume' nbpb=nbpb+1 zw2(ig,l+1)=0. @@ -438,5 +518,4 @@ if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l - RETURN END SUBROUTINE thermcell_plume