Index: /LMDZ6/trunk/libf/phylmd/thermcell_main.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_main.F90 (revision 3450) +++ /LMDZ6/trunk/libf/phylmd/thermcell_main.F90 (revision 3451) @@ -440,12 +440,35 @@ ! if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out -!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & - -! Gestion temporaire de plusieurs appels à thermcell_plume au travers -! de la variable iflag_thermals - -! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed + +!===================================================================== +! Old version of thermcell_plume in thermcell_plume_6A.F90 +! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding +! to the 5B and 6A versions used for CMIP5 and CMIP6. +! The latest was previously named thermcellV1_plume. +! The new thermcell_plume is a clean version (removing obsolete +! options) of thermcell_plume_6A. +! The 3 versions are controled by +! flag_thermals_ed <= 9 thermcell_plume_6A +! <= 19 thermcell_plume_5B +! else thermcell_plume (default 20 for convergence with 6A) +! Fredho +!===================================================================== + if (iflag_thermals_ed<=9) then ! print*,'THERM NOUVELLE/NOUVELLE Arnaud' + CALL thermcell_plume_6A(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, & + & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & + & ,lev_out,lunout1,igout) + + elseif (iflag_thermals_ed<=19) then +! print*,'THERM RIO et al 2010, version d Arnaud' + CALL thermcell_plume_5B(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, & + & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & + & ,lev_out,lunout1,igout) + else CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & @@ -453,13 +476,4 @@ & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) - - 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,& - & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & - & lalim,f0,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) - endif Index: /LMDZ6/trunk/libf/phylmd/thermcell_plume.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_plume.F90 (revision 3450) +++ /LMDZ6/trunk/libf/phylmd/thermcell_plume.F90 (revision 3451) @@ -1,4 +1,4 @@ ! -! $Id$ +! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $ ! SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, & @@ -9,5 +9,15 @@ ! & ,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) !-------------------------------------------------------------------------- USE IOIPSL, ONLY : getin @@ -82,5 +92,5 @@ real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev) - real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis + real zdz2,zdz3,lmel,entrbis,zdzbis real d_temp(ngrid) real ztv1,ztv2,factinv,zinv,zlmel @@ -283,196 +293,33 @@ !-------------------------------------------------- - if (iflag_thermals_ed.lt.8) then -!-------------------------------------------------- -!AJ052014: J'ai remplac?? la boucle do par un do while -! afin de faire moins de calcul dans la boucle -!-------------------------------------------------- - do while (zlmelup.gt.zltup) - lt=lt+1 - zlt=zlev(ig,lt) - zdz3=zlev(ig,lt+1)-zlt - zltdwn=zlt-zdz3/2 - zltup=zlt+zdz3/2 - enddo -!-------------------------------------------------- -!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors -! on cherche o?? se trouve l'altitude d'inversion -! en calculant ztv1 (interpolation de la valeur de -! theta au niveau lt en utilisant les niveaux lt-1 et -! lt-2) et ztv2 (interpolation avec les niveaux lt+1 -! et lt+2). Si theta r??ellement calcul??e au niveau lt -! comprise entre ztv1 et ztv2, alors il y a inversion -! et on calcule son altitude zinv en supposant que ztv(lt) -! est une combinaison lineaire de ztv1 et ztv2. -! Ensuite, on calcule la flottabilite en comparant -! la temperature de la couche l a celle de l'air situe -! l+lmel plus haut, ce qui necessite de savoir quel fraction -! de cet air est au-dessus ou en-dessous de l'inversion -!-------------------------------------------------- - atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2)) - btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & - & /(zlev(ig,lt-1)-zlev(ig,lt-2)) - atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1)) - btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & - & /(zlev(ig,lt+2)-zlev(ig,lt+1)) - - ztv1=atv1*zlt+btv1 - ztv2=atv2*zlt+btv2 - - if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then - -!-------------------------------------------------- -!AJ052014: D??calage de zinv qui est entre le haut -! et le bas de la couche lt -!-------------------------------------------------- - factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1) - zinv=zltdwn+zdz3*factinv - - - if (zlmeldwn.ge.zinv) then - ztv_est(ig,l)=atv2*zlmel+btv2 - zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & - & +(1.-fact_shell)*zbuoy(ig,l) - elseif (zlmelup.ge.zinv) then - ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2 - ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1 - ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1 - - zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- & - & ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- & - & ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l) - - else - ztv_est(ig,l)=atv1*zlmel+btv1 - zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & - & +(1.-fact_shell)*zbuoy(ig,l) - endif - - else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then - - if (zlmeldwn.gt.zltdwn) then - zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- & - & ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l) - else - zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & - & ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & - & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l) - - endif - -! zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & -! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & -! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*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)) - endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then - - else ! if (iflag_thermals_ed.lt.8) then - lt=l+1 + lt=l+1 + zlt=zlev(ig,lt) + zdz2=zlev(ig,lt)-zlev(ig,l) + + do while (lmel.gt.zdz2) + lt=lt+1 zlt=zlev(ig,lt) - zdz2=zlev(ig,lt)-zlev(ig,l) - - do while (lmel.gt.zdz2) - lt=lt+1 - zlt=zlev(ig,lt) - zdz2=zlt-zlev(ig,l) - enddo - zdz3=zlev(ig,lt+1)-zlt - zltdwn=zlev(ig,lt)-zdz3/2 - 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 + zdz2=zlt-zlev(ig,l) + enddo + zdz3=zlev(ig,lt+1)-zlt + zltdwn=zlev(ig,lt)-zdz3/2 + 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) !------------------------------------------------ !AJAM:nouveau calcul de w? !------------------------------------------------ - zdz=zlev(ig,l+1)-zlev(ig,l) - 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) - 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) & -! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) -! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & -! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) -! 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)) - -!-------------------------------------------------- -!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l) -!-------------------------------------------------- - if (iflag_thermals_ed==8) then -! Ancienne version -! 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)) - - w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) - -! Nouvelle version Arnaud - else -! 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)) - - w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) - -! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* & -! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* & -! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis)) - - - -! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact)) - -! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & -! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) - -! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) - - endif - - - if (iflag_thermals_ed<6) then - zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) -! fact_epsilon=0.0005/(zalpha+0.025)**0.5 -! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) - fact_epsilon=0.0002/(zalpha+0.1) - 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,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) - -! w_est(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)) - - w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) - - - endif -!-------------------------------------------------- -!AJ052014: J'ai comment? ce if plus n?cessaire puisqu' -!on fait max(0.0001,.....) -!-------------------------------------------------- - -! if (w_est(ig,l+1).lt.0.) then -! w_est(ig,l+1)=zw2(ig,l) -! w_est(ig,l+1)=0.0001 -! endif - + zdz=zlev(ig,l+1)-zlev(ig,l) + 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) + zdw2=afact*zbuoy(ig,l)/fact_epsilon + zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon + lm=Max(1,l-2) + w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) endif enddo @@ -488,36 +335,12 @@ ! 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) - -! 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 ) - - - -! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) - !========================================================================= ! 4. Calcul de l'entrainement et du detrainement !========================================================================= - -! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & -! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ) -! entrbis=entr_star(ig,l) - - if (iflag_thermals_ed.lt.6) then - fact_epsilon=0.0002/(zalpha+0.1) - endif - - detr_star(ig,l)=f_star(ig,l)*zdz & @@ -526,29 +349,15 @@ & + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) -! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ & -! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1) - - 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)) - - -! entr_star(ig,l)=f_star(ig,l)*zdz* ( & -! & mix0 * 0.1 / (zalpha+0.001) & -! & + MAX(entr_min, & -! & zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon + & -! & detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) - - -! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ & -! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1) - -! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* & -! & afact*zbuoy(ig,l)/zw2m & -! & - 1.*fact_epsilon) - + 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)) + 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)) + endif ! En dessous de lalim, on prend le max de alim_star et entr_star pour @@ -558,11 +367,5 @@ entr_star(ig,l)=0. endif -! 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) & + f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & & -detr_star(ig,l) @@ -606,54 +409,13 @@ zdzbis=zlev(ig,l)-zlev(ig,l-1) zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) -!!!!!!! fact_epsilon=0.002 - 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) -! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & -! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) -! lm=Max(1,l-2) -! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & -! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) -! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) -! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & -! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) -! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) -! 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)) - 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)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* & -! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis)) - - - if (iflag_thermals_ed.lt.6) then - zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l) -! fact_epsilon=0.0005/(zalpha+0.025)**0.5 -! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) - fact_epsilon=0.0002/(zalpha+0.1)**1 - 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,(zdz/zdzbis)*(exp(-zw2fact)* & -! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & -! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) -! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) - zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) - - endif - - + 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) endif 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 ! !=========================================================================== @@ -732,455 +494,2 @@ return end - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - SUBROUTINE thermcellV1_plume(itap,ngrid,klev,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, & -& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & -& ,lev_out,lunout1,igout) -!& ,lev_out,lunout1,igout,zbuoy,zbuoyjam) - -!-------------------------------------------------------------------------- -!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance -! Version conforme a l'article de Rio et al. 2010. -! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin -!-------------------------------------------------------------------------- - - USE print_control_mod, ONLY: prt_level - IMPLICIT NONE - -#include "YOMCST.h" -#include "YOETHF.h" -#include "FCTTRE.h" -#include "thermcell.h" - - INTEGER itap - INTEGER lunout1,igout - INTEGER ngrid,klev - REAL ptimestep - REAL ztv(ngrid,klev) - REAL zthl(ngrid,klev) - REAL po(ngrid,klev) - REAL zl(ngrid,klev) - REAL rhobarz(ngrid,klev) - REAL zlev(ngrid,klev+1) - REAL pplev(ngrid,klev+1) - REAL pphi(ngrid,klev) - REAL zpspsk(ngrid,klev) - REAL alim_star(ngrid,klev) - REAL f0(ngrid) - INTEGER lalim(ngrid) - integer lev_out ! niveau pour les print - integer nbpb - - real alim_star_tot(ngrid) - - REAL ztva(ngrid,klev) - REAL ztla(ngrid,klev) - REAL zqla(ngrid,klev) - REAL zqta(ngrid,klev) - REAL zha(ngrid,klev) - - REAL detr_star(ngrid,klev) - REAL coefc - REAL entr_star(ngrid,klev) - REAL detr(ngrid,klev) - REAL entr(ngrid,klev) - - REAL csc(ngrid,klev) - - REAL zw2(ngrid,klev+1) - REAL w_est(ngrid,klev+1) - REAL f_star(ngrid,klev+1) - REAL wa_moy(ngrid,klev+1) - - REAL ztva_est(ngrid,klev) - REAL zqla_est(ngrid,klev) - REAL zqsatth(ngrid,klev) - REAL zta_est(ngrid,klev) - REAL zbuoyjam(ngrid,klev) - REAL ztemp(ngrid),zqsat(ngrid) - REAL zdw2 - REAL zw2modif - REAL zw2fact - REAL zeps(ngrid,klev) - - REAL linter(ngrid) - INTEGER lmix(ngrid) - INTEGER lmix_bis(ngrid) - REAL wmaxa(ngrid) - - INTEGER ig,l,k - - real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m - real zbuoybis - real zcor,zdelta,zcvm5,qlbef,zdz2 - 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,fact_epsilon2 - REAL c2(ngrid,klev) - Zsat=.false. -! Initialisation - - RLvCp = RLVTT/RCPD - fact_epsilon=0.002 - betalpha=0.9 - afact=2./3. - - zbetalpha=betalpha/(1.+betalpha) - - -! Initialisations des variables reeles -if (1==1) then - ztva(:,:)=ztv(:,:) - ztva_est(:,:)=ztva(:,:) - ztla(:,:)=zthl(:,:) - zqta(:,:)=po(:,:) - zha(:,:) = ztva(:,:) -else - ztva(:,:)=0. - ztva_est(:,:)=0. - ztla(:,:)=0. - zqta(:,:)=0. - zha(:,:) =0. -endif - - zqla_est(:,:)=0. - zqsatth(:,:)=0. - zqla(:,:)=0. - detr_star(:,:)=0. - entr_star(:,:)=0. - alim_star(:,:)=0. - alim_star_tot(:)=0. - csc(:,:)=0. - detr(:,:)=0. - 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. -! Initialisation des variables entieres - lmix(:)=1 - lmix_bis(:)=2 - wmaxa(:)=0. - lalim(:)=1 - - -!------------------------------------------------------------------------- -! On ne considere comme actif que les colonnes dont les deux premieres -! couches sont instables. -!------------------------------------------------------------------------- - active(:)=ztv(:,1)>ztv(:,2) - -!------------------------------------------------------------------------- -! Definition de l'alimentation a l'origine dans thermcell_init -!------------------------------------------------------------------------- - do l=1,klev-1 - do ig=1,ngrid - if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then - alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & - & *sqrt(zlev(ig,l+1)) - lalim(ig)=l+1 - alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) - endif - enddo - enddo - do l=1,klev - do ig=1,ngrid - if (alim_star_tot(ig) > 1.e-10 ) then - alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) - endif - enddo - enddo - alim_star_tot(:)=1. - - - -!------------------------------------------------------------------------------ -! 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 -! dans une couche l>1 -!------------------------------------------------------------------------------ -do ig=1,ngrid -! Le panache va prendre au debut les caracteristiques de l'air contenu -! dans cette couche. - if (active(ig)) then - ztla(ig,1)=zthl(ig,1) - zqta(ig,1)=po(ig,1) - zqla(ig,1)=zl(ig,1) -!cr: attention, prise en compte de f*(1)=1 - f_star(ig,2)=alim_star(ig,1) - zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & -& *(zlev(ig,2)-zlev(ig,1)) & -& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) - w_est(ig,2)=zw2(ig,2) - endif -enddo -! - -!============================================================================== -!boucle de calcul de la vitesse verticale dans le thermique -!============================================================================== -do l=2,klev-1 -!============================================================================== - - -! On decide si le thermique est encore actif ou non -! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test - do ig=1,ngrid - active(ig)=active(ig) & -& .and. zw2(ig,l)>1.e-10 & -& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 - enddo - - - -!--------------------------------------------------------------------------- -! 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 -!--------------------------------------------------------------------------- - - 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) - 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)) - -!------------------------------------------------ -!AJAM:nouveau calcul de w? -!------------------------------------------------ - zdz=zlev(ig,l+1)-zlev(ig,l) - zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) - - zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) - zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon) - w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) - - - if (w_est(ig,l+1).lt.0.) then - w_est(ig,l+1)=zw2(ig,l) - endif - endif - enddo - - -!------------------------------------------------- -!calcul des taux d'entrainement et de detrainement -!------------------------------------------------- - - 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) - 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) - - - entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*MAX(0., & - & afact*zbuoybis/zw2m - fact_epsilon ) - - - detr_star(ig,l)=f_star(ig,l)*zdz & - & *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m & - & + 0.012*(zdqt(ig,l)/zw2m)**0.5 ) - -! En dessous de lalim, on prend le max de alim_star et entr_star pour -! alim_star et 0 sinon - if (l.lt.lalim(ig)) then - alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) - entr_star(ig,l)=0. - endif - -! 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) - - endif - enddo - - -!---------------------------------------------------------------------------- -!calcul de la vitesse verticale en melangeant Tl et qt du thermique -!--------------------------------------------------------------------------- - activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 - do ig=1,ngrid - if (activetmp(ig)) then - Zsat=.false. - ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & - & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & - & /(f_star(ig,l+1)+detr_star(ig,l)) - zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & - & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & - & /(f_star(ig,l+1)+detr_star(ig,l)) - - endif - enddo - - 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) -!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)) - zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) - zdz=zlev(ig,l+1)-zlev(ig,l) - zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) - - zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) - zdw2=afact*zbuoy(ig,l)/(fact_epsilon) - zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) - endif - enddo - - if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l -! -!--------------------------------------------------------------------------- -!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max -!--------------------------------------------------------------------------- - - nbpb=0 - 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. - linter(ig)=l+1 - endif - - if (zw2(ig,l+1).lt.0.) then - linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & - & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) - zw2(ig,l+1)=0. - elseif (f_star(ig,l+1).lt.0.) then - linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & - & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) -! print*,"linter plume", linter(ig) - zw2(ig,l+1)=0. - endif - - wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) - - if (wa_moy(ig,l+1).gt.wmaxa(ig)) then -! lmix est le niveau de la couche ou w (wa_moy) est maximum -!on rajoute le calcul de lmix_bis - if (zqla(ig,l).lt.1.e-10) then - lmix_bis(ig)=l+1 - endif - lmix(ig)=l+1 - wmaxa(ig)=wa_moy(ig,l+1) - endif - enddo - - if (nbpb>0) then - print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' - endif - -!========================================================================= -! FIN DE LA BOUCLE VERTICALE - enddo -!========================================================================= - -!on recalcule alim_star_tot - do ig=1,ngrid - alim_star_tot(ig)=0. - enddo - do ig=1,ngrid - do l=1,lalim(ig)-1 - alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) - enddo - enddo - - - if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l - -#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 Index: /LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90 (revision 3451) +++ /LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90 (revision 3451) @@ -0,0 +1,1186 @@ +! +! $Id$ +! + SUBROUTINE thermcell_plume_6A(itap,ngrid,klev,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, & + & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & + & ,lev_out,lunout1,igout) +! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) +!-------------------------------------------------------------------------- +!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance +!-------------------------------------------------------------------------- +USE IOIPSL, ONLY : getin +USE ioipsl_getin_p_mod, ONLY : getin_p + + USE print_control_mod, ONLY: prt_level + IMPLICIT NONE + +#include "YOMCST.h" +#include "YOETHF.h" +#include "FCTTRE.h" +#include "thermcell.h" + + INTEGER itap + INTEGER lunout1,igout + INTEGER ngrid,klev + REAL ptimestep + REAL ztv(ngrid,klev) + REAL zthl(ngrid,klev) + REAL po(ngrid,klev) + REAL zl(ngrid,klev) + REAL rhobarz(ngrid,klev) + REAL zlev(ngrid,klev+1) + REAL pplev(ngrid,klev+1) + REAL pphi(ngrid,klev) + REAL zpspsk(ngrid,klev) + REAL alim_star(ngrid,klev) + REAL f0(ngrid) + INTEGER lalim(ngrid) + integer lev_out ! niveau pour les print + integer nbpb + + real alim_star_tot(ngrid) + + REAL ztva(ngrid,klev) + REAL ztla(ngrid,klev) + REAL zqla(ngrid,klev) + REAL zqta(ngrid,klev) + REAL zha(ngrid,klev) + + REAL detr_star(ngrid,klev) + REAL coefc + REAL entr_star(ngrid,klev) + REAL detr(ngrid,klev) + REAL entr(ngrid,klev) + + REAL csc(ngrid,klev) + + REAL zw2(ngrid,klev+1) + REAL w_est(ngrid,klev+1) + REAL f_star(ngrid,klev+1) + REAL wa_moy(ngrid,klev+1) + + 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 ztemp(ngrid),zqsat(ngrid) + REAL zdw2,zdw2bis + REAL zw2modif + REAL zw2fact,zw2factbis + REAL zeps(ngrid,klev) + + REAL linter(ngrid) + INTEGER lmix(ngrid) + INTEGER lmix_bis(ngrid) + REAL wmaxa(ngrid) + + INTEGER ig,l,k,lt,it,lm + + real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m + real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev) + real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis + real d_temp(ngrid) + real ztv1,ztv2,factinv,zinv,zlmel + real zlmelup,zlmeldwn,zlt,zltdwn,zltup + real atv1,atv2,btv1,btv2 + real ztv_est1,ztv_est2 + real zcor,zdelta,zcvm5,qlbef + real zbetalpha, coefzlmel + real eps + REAL REPS,RLvCp,DDT0 + PARAMETER (DDT0=.01) + logical Zsat + LOGICAL active(ngrid),activetmp(ngrid) + REAL fact_gamma,fact_gamma2,fact_epsilon2 + + REAL, SAVE :: fact_epsilon=0.002 + REAL, SAVE :: betalpha=0.9 + REAL, SAVE :: afact=2./3. + REAL, SAVE :: fact_shell=1. + REAL,SAVE :: detr_min=1.e-5 + REAL,SAVE :: entr_min=1.e-5 + REAL,SAVE :: detr_q_coef=0.012 + REAL,SAVE :: detr_q_power=0.5 + REAL,SAVE :: mix0=0. + INTEGER,SAVE :: thermals_flag_alim=0 + +!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell) +!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power) +!$OMP THREADPRIVATE( mix0, thermals_flag_alim) + + LOGICAL, SAVE :: first=.true. + !$OMP THREADPRIVATE(first) + + + REAL c2(ngrid,klev) + + if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11' + Zsat=.false. +! Initialisation + + RLvCp = RLVTT/RCPD + IF (first) THEN + + CALL getin_p('thermals_fact_epsilon',fact_epsilon) + CALL getin_p('thermals_betalpha',betalpha) + CALL getin_p('thermals_afact',afact) + CALL getin_p('thermals_fact_shell',fact_shell) + CALL getin_p('thermals_detr_min',detr_min) + CALL getin_p('thermals_entr_min',entr_min) + CALL getin_p('thermals_detr_q_coef',detr_q_coef) + CALL getin_p('thermals_detr_q_power',detr_q_power) + CALL getin_p('thermals_mix0',mix0) + CALL getin_p('thermals_flag_alim',thermals_flag_alim) + + + first=.false. + ENDIF + + 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. + zqta(:,:)=0. + zha(:,:) =0. +endif + + zqla_est(:,:)=0. + zqsatth(:,:)=0. + zqla(:,:)=0. + detr_star(:,:)=0. + entr_star(:,:)=0. + alim_star(:,:)=0. + alim_star_tot(:)=0. + csc(:,:)=0. + detr(:,:)=0. + 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. +! Initialisation des variables entieres + lmix(:)=1 + lmix_bis(:)=2 + wmaxa(:)=0. + + +!------------------------------------------------------------------------- +! On ne considere comme actif que les colonnes dont les deux premieres +! couches sont instables. +!------------------------------------------------------------------------- + + active(:)=ztv(:,1)>ztv(:,2) + d_temp(:)=0. ! Pour activer un contraste de temperature a la base + ! du panache +! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main + CALL thermcell_alim(thermals_flag_alim,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim) + +!------------------------------------------------------------------------------ +! 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 +! dans une couche l>1 +!------------------------------------------------------------------------------ +do ig=1,ngrid +! Le panache va prendre au debut les caracteristiques de l'air contenu +! dans cette couche. + if (active(ig)) then + ztla(ig,1)=zthl(ig,1) + zqta(ig,1)=po(ig,1) + zqla(ig,1)=zl(ig,1) +!cr: attention, prise en compte de f*(1)=1 + f_star(ig,2)=alim_star(ig,1) + zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & +& *(zlev(ig,2)-zlev(ig,1)) & +& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) + w_est(ig,2)=zw2(ig,2) + endif +enddo +! + +!============================================================================== +!boucle de calcul de la vitesse verticale dans le thermique +!============================================================================== +do l=2,klev-1 +!============================================================================== + + +! On decide si le thermique est encore actif ou non +! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test + do ig=1,ngrid + active(ig)=active(ig) & +& .and. zw2(ig,l)>1.e-10 & +& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 + enddo + + + +!--------------------------------------------------------------------------- +! 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 +!--------------------------------------------------------------------------- + + 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) + 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)) + + +!Modif AJAM + + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + zdz=zlev(ig,l+1)-zlev(ig,l) + lmel=fact_thermals_ed_dz*zlev(ig,l) +! lmel=0.09*zlev(ig,l) + zlmel=zlev(ig,l)+lmel + zlmelup=zlmel+(zdz/2) + zlmeldwn=zlmel-(zdz/2) + + lt=l+1 + zlt=zlev(ig,lt) + zdz3=zlev(ig,lt+1)-zlt + zltdwn=zlt-zdz3/2 + zltup=zlt+zdz3/2 + +!========================================================================= +! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus +!========================================================================= + +!-------------------------------------------------- + if (iflag_thermals_ed.lt.8) then +!-------------------------------------------------- +!AJ052014: J'ai remplac?? la boucle do par un do while +! afin de faire moins de calcul dans la boucle +!-------------------------------------------------- + do while (zlmelup.gt.zltup) + lt=lt+1 + zlt=zlev(ig,lt) + zdz3=zlev(ig,lt+1)-zlt + zltdwn=zlt-zdz3/2 + zltup=zlt+zdz3/2 + enddo +!-------------------------------------------------- +!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors +! on cherche o?? se trouve l'altitude d'inversion +! en calculant ztv1 (interpolation de la valeur de +! theta au niveau lt en utilisant les niveaux lt-1 et +! lt-2) et ztv2 (interpolation avec les niveaux lt+1 +! et lt+2). Si theta r??ellement calcul??e au niveau lt +! comprise entre ztv1 et ztv2, alors il y a inversion +! et on calcule son altitude zinv en supposant que ztv(lt) +! est une combinaison lineaire de ztv1 et ztv2. +! Ensuite, on calcule la flottabilite en comparant +! la temperature de la couche l a celle de l'air situe +! l+lmel plus haut, ce qui necessite de savoir quel fraction +! de cet air est au-dessus ou en-dessous de l'inversion +!-------------------------------------------------- + atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2)) + btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & + & /(zlev(ig,lt-1)-zlev(ig,lt-2)) + atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1)) + btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & + & /(zlev(ig,lt+2)-zlev(ig,lt+1)) + + ztv1=atv1*zlt+btv1 + ztv2=atv2*zlt+btv2 + + if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then + +!-------------------------------------------------- +!AJ052014: D??calage de zinv qui est entre le haut +! et le bas de la couche lt +!-------------------------------------------------- + factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1) + zinv=zltdwn+zdz3*factinv + + + if (zlmeldwn.ge.zinv) then + ztv_est(ig,l)=atv2*zlmel+btv2 + zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & + & +(1.-fact_shell)*zbuoy(ig,l) + elseif (zlmelup.ge.zinv) then + ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2 + ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1 + ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1 + + zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- & + & ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- & + & ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l) + + else + ztv_est(ig,l)=atv1*zlmel+btv1 + zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & + & +(1.-fact_shell)*zbuoy(ig,l) + endif + + else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then + + if (zlmeldwn.gt.zltdwn) then + zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- & + & ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l) + else + zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & + & ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & + & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l) + + endif + +! zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & +! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & +! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*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)) + endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then + + else ! if (iflag_thermals_ed.lt.8) then + lt=l+1 + zlt=zlev(ig,lt) + zdz2=zlev(ig,lt)-zlev(ig,l) + + do while (lmel.gt.zdz2) + lt=lt+1 + zlt=zlev(ig,lt) + zdz2=zlt-zlev(ig,l) + enddo + zdz3=zlev(ig,lt+1)-zlt + zltdwn=zlev(ig,lt)-zdz3/2 + 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 + +!------------------------------------------------ +!AJAM:nouveau calcul de w? +!------------------------------------------------ + zdz=zlev(ig,l+1)-zlev(ig,l) + 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) + 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) & +! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) +! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & +! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) +! 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)) + +!-------------------------------------------------- +!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l) +!-------------------------------------------------- + if (iflag_thermals_ed==8) then +! Ancienne version +! 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)) + + w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) + +! Nouvelle version Arnaud + else +! 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)) + + w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) + +! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* & +! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* & +! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis)) + + + +! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact)) + +! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & +! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) + +! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) + + endif + + + if (iflag_thermals_ed<6) then + zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) +! fact_epsilon=0.0005/(zalpha+0.025)**0.5 +! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) + fact_epsilon=0.0002/(zalpha+0.1) + 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,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) + +! w_est(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)) + + w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) + + + endif +!-------------------------------------------------- +!AJ052014: J'ai comment? ce if plus n?cessaire puisqu' +!on fait max(0.0001,.....) +!-------------------------------------------------- + +! if (w_est(ig,l+1).lt.0.) then +! 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 +!------------------------------------------------- + + 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(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 ) + + + +! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + +!========================================================================= +! 4. Calcul de l'entrainement et du detrainement +!========================================================================= + +! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & +! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ) +! entrbis=entr_star(ig,l) + + if (iflag_thermals_ed.lt.6) then + fact_epsilon=0.0002/(zalpha+0.1) + endif + + + + 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)) + +! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ & +! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1) + + 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)) + + +! entr_star(ig,l)=f_star(ig,l)*zdz* ( & +! & mix0 * 0.1 / (zalpha+0.001) & +! & + MAX(entr_min, & +! & zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon + & +! & detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) + + +! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ & +! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1) + +! 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 + if (l.lt.lalim(ig)) then + alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) + entr_star(ig,l)=0. + endif +! 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) + + endif + enddo + + +!============================================================================ +! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique +!=========================================================================== + + activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 + do ig=1,ngrid + if (activetmp(ig)) then + Zsat=.false. + ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & + & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & + & /(f_star(ig,l+1)+detr_star(ig,l)) + zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & + & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & + & /(f_star(ig,l+1)+detr_star(ig,l)) + + endif + enddo + + 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) +!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)) + zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) + zdz=zlev(ig,l+1)-zlev(ig,l) + zdzbis=zlev(ig,l)-zlev(ig,l-1) + zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) +!!!!!!! fact_epsilon=0.002 + 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) +! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & +! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) +! lm=Max(1,l-2) +! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & +! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) +! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) +! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & +! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) +! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) +! 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)) + 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)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* & +! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis)) + + + if (iflag_thermals_ed.lt.6) then + zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l) +! fact_epsilon=0.0005/(zalpha+0.025)**0.5 +! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) + fact_epsilon=0.0002/(zalpha+0.1)**1 + 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,(zdz/zdzbis)*(exp(-zw2fact)* & +! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & +! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) +! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) + zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) + + endif + + + endif + enddo + + if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l +! +!=========================================================================== +! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max +!=========================================================================== + + nbpb=0 + 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. + linter(ig)=l+1 + endif + + if (zw2(ig,l+1).lt.0.) then + linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & + & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) + zw2(ig,l+1)=0. +!+CR:04/05/12:correction calcul linter pour calcul de zmax continu + elseif (f_star(ig,l+1).lt.0.) then + linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & + & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) + zw2(ig,l+1)=0. +!fin CR:04/05/12 + endif + + wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) + + if (wa_moy(ig,l+1).gt.wmaxa(ig)) then +! lmix est le niveau de la couche ou w (wa_moy) est maximum +!on rajoute le calcul de lmix_bis + if (zqla(ig,l).lt.1.e-10) then + lmix_bis(ig)=l+1 + endif + lmix(ig)=l+1 + wmaxa(ig)=wa_moy(ig,l+1) + endif + enddo + + if (nbpb>0) then + print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' + endif + +!========================================================================= +! FIN DE LA BOUCLE VERTICALE + enddo +!========================================================================= + +!on recalcule alim_star_tot + do ig=1,ngrid + alim_star_tot(ig)=0. + enddo + do ig=1,ngrid + do l=1,lalim(ig)-1 + alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) + enddo + enddo + + + if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l + +#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 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE thermcell_plume_5B(itap,ngrid,klev,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, & +& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & +& ,lev_out,lunout1,igout) +!& ,lev_out,lunout1,igout,zbuoy,zbuoyjam) + +!-------------------------------------------------------------------------- +!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance +! Version conforme a l'article de Rio et al. 2010. +! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin +!-------------------------------------------------------------------------- + + USE print_control_mod, ONLY: prt_level + IMPLICIT NONE + +#include "YOMCST.h" +#include "YOETHF.h" +#include "FCTTRE.h" +#include "thermcell.h" + + INTEGER itap + INTEGER lunout1,igout + INTEGER ngrid,klev + REAL ptimestep + REAL ztv(ngrid,klev) + REAL zthl(ngrid,klev) + REAL po(ngrid,klev) + REAL zl(ngrid,klev) + REAL rhobarz(ngrid,klev) + REAL zlev(ngrid,klev+1) + REAL pplev(ngrid,klev+1) + REAL pphi(ngrid,klev) + REAL zpspsk(ngrid,klev) + REAL alim_star(ngrid,klev) + REAL f0(ngrid) + INTEGER lalim(ngrid) + integer lev_out ! niveau pour les print + integer nbpb + + real alim_star_tot(ngrid) + + REAL ztva(ngrid,klev) + REAL ztla(ngrid,klev) + REAL zqla(ngrid,klev) + REAL zqta(ngrid,klev) + REAL zha(ngrid,klev) + + REAL detr_star(ngrid,klev) + REAL coefc + REAL entr_star(ngrid,klev) + REAL detr(ngrid,klev) + REAL entr(ngrid,klev) + + REAL csc(ngrid,klev) + + REAL zw2(ngrid,klev+1) + REAL w_est(ngrid,klev+1) + REAL f_star(ngrid,klev+1) + REAL wa_moy(ngrid,klev+1) + + REAL ztva_est(ngrid,klev) + REAL zqla_est(ngrid,klev) + REAL zqsatth(ngrid,klev) + REAL zta_est(ngrid,klev) + REAL zbuoyjam(ngrid,klev) + REAL ztemp(ngrid),zqsat(ngrid) + REAL zdw2 + REAL zw2modif + REAL zw2fact + REAL zeps(ngrid,klev) + + REAL linter(ngrid) + INTEGER lmix(ngrid) + INTEGER lmix_bis(ngrid) + REAL wmaxa(ngrid) + + INTEGER ig,l,k + + real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m + real zbuoybis + real zcor,zdelta,zcvm5,qlbef,zdz2 + 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,fact_epsilon2 + REAL c2(ngrid,klev) + Zsat=.false. +! Initialisation + + RLvCp = RLVTT/RCPD + fact_epsilon=0.002 + betalpha=0.9 + afact=2./3. + + zbetalpha=betalpha/(1.+betalpha) + + +! Initialisations des variables reeles +if (1==1) then + ztva(:,:)=ztv(:,:) + ztva_est(:,:)=ztva(:,:) + ztla(:,:)=zthl(:,:) + zqta(:,:)=po(:,:) + zha(:,:) = ztva(:,:) +else + ztva(:,:)=0. + ztva_est(:,:)=0. + ztla(:,:)=0. + zqta(:,:)=0. + zha(:,:) =0. +endif + + zqla_est(:,:)=0. + zqsatth(:,:)=0. + zqla(:,:)=0. + detr_star(:,:)=0. + entr_star(:,:)=0. + alim_star(:,:)=0. + alim_star_tot(:)=0. + csc(:,:)=0. + detr(:,:)=0. + 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. +! Initialisation des variables entieres + lmix(:)=1 + lmix_bis(:)=2 + wmaxa(:)=0. + lalim(:)=1 + + +!------------------------------------------------------------------------- +! On ne considere comme actif que les colonnes dont les deux premieres +! couches sont instables. +!------------------------------------------------------------------------- + active(:)=ztv(:,1)>ztv(:,2) + +!------------------------------------------------------------------------- +! Definition de l'alimentation a l'origine dans thermcell_init +!------------------------------------------------------------------------- + do l=1,klev-1 + do ig=1,ngrid + if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then + alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & + & *sqrt(zlev(ig,l+1)) + lalim(ig)=l+1 + alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) + endif + enddo + enddo + do l=1,klev + do ig=1,ngrid + if (alim_star_tot(ig) > 1.e-10 ) then + alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) + endif + enddo + enddo + alim_star_tot(:)=1. + + + +!------------------------------------------------------------------------------ +! 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 +! dans une couche l>1 +!------------------------------------------------------------------------------ +do ig=1,ngrid +! Le panache va prendre au debut les caracteristiques de l'air contenu +! dans cette couche. + if (active(ig)) then + ztla(ig,1)=zthl(ig,1) + zqta(ig,1)=po(ig,1) + zqla(ig,1)=zl(ig,1) +!cr: attention, prise en compte de f*(1)=1 + f_star(ig,2)=alim_star(ig,1) + zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & +& *(zlev(ig,2)-zlev(ig,1)) & +& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) + w_est(ig,2)=zw2(ig,2) + endif +enddo +! + +!============================================================================== +!boucle de calcul de la vitesse verticale dans le thermique +!============================================================================== +do l=2,klev-1 +!============================================================================== + + +! On decide si le thermique est encore actif ou non +! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test + do ig=1,ngrid + active(ig)=active(ig) & +& .and. zw2(ig,l)>1.e-10 & +& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 + enddo + + + +!--------------------------------------------------------------------------- +! 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 +!--------------------------------------------------------------------------- + + 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) + 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)) + +!------------------------------------------------ +!AJAM:nouveau calcul de w? +!------------------------------------------------ + zdz=zlev(ig,l+1)-zlev(ig,l) + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + + zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) + zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon) + w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) + + + if (w_est(ig,l+1).lt.0.) then + w_est(ig,l+1)=zw2(ig,l) + endif + endif + enddo + + +!------------------------------------------------- +!calcul des taux d'entrainement et de detrainement +!------------------------------------------------- + + 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) + 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) + + + entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*MAX(0., & + & afact*zbuoybis/zw2m - fact_epsilon ) + + + detr_star(ig,l)=f_star(ig,l)*zdz & + & *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m & + & + 0.012*(zdqt(ig,l)/zw2m)**0.5 ) + +! En dessous de lalim, on prend le max de alim_star et entr_star pour +! alim_star et 0 sinon + if (l.lt.lalim(ig)) then + alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) + entr_star(ig,l)=0. + endif + +! 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) + + endif + enddo + + +!---------------------------------------------------------------------------- +!calcul de la vitesse verticale en melangeant Tl et qt du thermique +!--------------------------------------------------------------------------- + activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 + do ig=1,ngrid + if (activetmp(ig)) then + Zsat=.false. + ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & + & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & + & /(f_star(ig,l+1)+detr_star(ig,l)) + zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & + & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & + & /(f_star(ig,l+1)+detr_star(ig,l)) + + endif + enddo + + 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) +!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)) + zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) + zdz=zlev(ig,l+1)-zlev(ig,l) + zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) + + zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) + zdw2=afact*zbuoy(ig,l)/(fact_epsilon) + zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) + endif + enddo + + if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l +! +!--------------------------------------------------------------------------- +!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max +!--------------------------------------------------------------------------- + + nbpb=0 + 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. + linter(ig)=l+1 + endif + + if (zw2(ig,l+1).lt.0.) then + linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & + & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) + zw2(ig,l+1)=0. + elseif (f_star(ig,l+1).lt.0.) then + linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & + & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) +! print*,"linter plume", linter(ig) + zw2(ig,l+1)=0. + endif + + wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) + + if (wa_moy(ig,l+1).gt.wmaxa(ig)) then +! lmix est le niveau de la couche ou w (wa_moy) est maximum +!on rajoute le calcul de lmix_bis + if (zqla(ig,l).lt.1.e-10) then + lmix_bis(ig)=l+1 + endif + lmix(ig)=l+1 + wmaxa(ig)=wa_moy(ig,l+1) + endif + enddo + + if (nbpb>0) then + print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' + endif + +!========================================================================= +! FIN DE LA BOUCLE VERTICALE + enddo +!========================================================================= + +!on recalcule alim_star_tot + do ig=1,ngrid + alim_star_tot(ig)=0. + enddo + do ig=1,ngrid + do l=1,lalim(ig)-1 + alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) + enddo + enddo + + + if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l + +#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