source: LMDZ6/trunk/libf/phylmd/thermcell_plume_6A.F90 @ 3761

Last change on this file since 3761 was 3451, checked in by fhourdin, 6 years ago

Saving 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

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 40.7 KB
Line 
1!
2! $Id: thermcell_plume_6A.F90 3451 2019-01-27 11:07:30Z adurocher $
3!
4      SUBROUTINE thermcell_plume_6A(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
5     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
6     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
7     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8    &           ,lev_out,lunout1,igout)
9!     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
10!--------------------------------------------------------------------------
11!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12!--------------------------------------------------------------------------
13USE IOIPSL, ONLY : getin
14USE ioipsl_getin_p_mod, ONLY : getin_p
15
16       USE print_control_mod, ONLY: prt_level
17       IMPLICIT NONE
18
19#include "YOMCST.h"
20#include "YOETHF.h"
21#include "FCTTRE.h"
22#include "thermcell.h"
23
24      INTEGER itap
25      INTEGER lunout1,igout
26      INTEGER ngrid,klev
27      REAL ptimestep
28      REAL ztv(ngrid,klev)
29      REAL zthl(ngrid,klev)
30      REAL po(ngrid,klev)
31      REAL zl(ngrid,klev)
32      REAL rhobarz(ngrid,klev)
33      REAL zlev(ngrid,klev+1)
34      REAL pplev(ngrid,klev+1)
35      REAL pphi(ngrid,klev)
36      REAL zpspsk(ngrid,klev)
37      REAL alim_star(ngrid,klev)
38      REAL f0(ngrid)
39      INTEGER lalim(ngrid)
40      integer lev_out                           ! niveau pour les print
41      integer nbpb
42   
43      real alim_star_tot(ngrid)
44
45      REAL ztva(ngrid,klev)
46      REAL ztla(ngrid,klev)
47      REAL zqla(ngrid,klev)
48      REAL zqta(ngrid,klev)
49      REAL zha(ngrid,klev)
50
51      REAL detr_star(ngrid,klev)
52      REAL coefc
53      REAL entr_star(ngrid,klev)
54      REAL detr(ngrid,klev)
55      REAL entr(ngrid,klev)
56
57      REAL csc(ngrid,klev)
58
59      REAL zw2(ngrid,klev+1)
60      REAL w_est(ngrid,klev+1)
61      REAL f_star(ngrid,klev+1)
62      REAL wa_moy(ngrid,klev+1)
63
64      REAL ztva_est(ngrid,klev)
65      REAL ztv_est(ngrid,klev)
66      REAL zqla_est(ngrid,klev)
67      REAL zqsatth(ngrid,klev)
68      REAL zta_est(ngrid,klev)
69      REAL ztemp(ngrid),zqsat(ngrid)
70      REAL zdw2,zdw2bis
71      REAL zw2modif
72      REAL zw2fact,zw2factbis
73      REAL zeps(ngrid,klev)
74
75      REAL linter(ngrid)
76      INTEGER lmix(ngrid)
77      INTEGER lmix_bis(ngrid)
78      REAL    wmaxa(ngrid)
79
80      INTEGER ig,l,k,lt,it,lm
81
82      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
83      real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev)
84      real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
85      real d_temp(ngrid)
86      real ztv1,ztv2,factinv,zinv,zlmel
87      real zlmelup,zlmeldwn,zlt,zltdwn,zltup
88      real atv1,atv2,btv1,btv2
89      real ztv_est1,ztv_est2
90      real zcor,zdelta,zcvm5,qlbef
91      real zbetalpha, coefzlmel
92      real eps
93      REAL REPS,RLvCp,DDT0
94      PARAMETER (DDT0=.01)
95      logical Zsat
96      LOGICAL active(ngrid),activetmp(ngrid)
97      REAL fact_gamma,fact_gamma2,fact_epsilon2
98
99      REAL, SAVE :: fact_epsilon=0.002
100      REAL, SAVE :: betalpha=0.9
101      REAL, SAVE :: afact=2./3.
102      REAL, SAVE :: fact_shell=1.
103      REAL,SAVE :: detr_min=1.e-5
104      REAL,SAVE :: entr_min=1.e-5
105      REAL,SAVE :: detr_q_coef=0.012
106      REAL,SAVE :: detr_q_power=0.5
107      REAL,SAVE :: mix0=0.
108      INTEGER,SAVE :: thermals_flag_alim=0
109
110!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell)
111!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power)
112!$OMP THREADPRIVATE( mix0, thermals_flag_alim)
113
114      LOGICAL, SAVE :: first=.true.
115  !$OMP THREADPRIVATE(first)
116
117
118      REAL c2(ngrid,klev)
119
120      if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
121      Zsat=.false.
122! Initialisation
123
124      RLvCp = RLVTT/RCPD
125      IF (first) THEN
126
127     CALL getin_p('thermals_fact_epsilon',fact_epsilon)
128     CALL getin_p('thermals_betalpha',betalpha)
129     CALL getin_p('thermals_afact',afact)
130     CALL getin_p('thermals_fact_shell',fact_shell)
131     CALL getin_p('thermals_detr_min',detr_min)
132     CALL getin_p('thermals_entr_min',entr_min)
133     CALL getin_p('thermals_detr_q_coef',detr_q_coef)
134     CALL getin_p('thermals_detr_q_power',detr_q_power)
135     CALL getin_p('thermals_mix0',mix0)
136     CALL getin_p('thermals_flag_alim',thermals_flag_alim)
137
138
139      first=.false.
140      ENDIF
141
142      zbetalpha=betalpha/(1.+betalpha)
143
144
145! Initialisations des variables r?elles
146if (1==1) then
147      ztva(:,:)=ztv(:,:)
148      ztva_est(:,:)=ztva(:,:)
149      ztv_est(:,:)=ztv(:,:)
150      ztla(:,:)=zthl(:,:)
151      zqta(:,:)=po(:,:)
152      zqla(:,:)=0.
153      zha(:,:) = ztva(:,:)
154else
155      ztva(:,:)=0.
156      ztv_est(:,:)=0.
157      ztva_est(:,:)=0.
158      ztla(:,:)=0.
159      zqta(:,:)=0.
160      zha(:,:) =0.
161endif
162
163      zqla_est(:,:)=0.
164      zqsatth(:,:)=0.
165      zqla(:,:)=0.
166      detr_star(:,:)=0.
167      entr_star(:,:)=0.
168      alim_star(:,:)=0.
169      alim_star_tot(:)=0.
170      csc(:,:)=0.
171      detr(:,:)=0.
172      entr(:,:)=0.
173      zw2(:,:)=0.
174      zbuoy(:,:)=0.
175      zbuoyjam(:,:)=0.
176      gamma(:,:)=0.
177      zeps(:,:)=0.
178      w_est(:,:)=0.
179      f_star(:,:)=0.
180      wa_moy(:,:)=0.
181      linter(:)=1.
182!     linter(:)=1.
183! Initialisation des variables entieres
184      lmix(:)=1
185      lmix_bis(:)=2
186      wmaxa(:)=0.
187
188
189!-------------------------------------------------------------------------
190! On ne considere comme actif que les colonnes dont les deux premieres
191! couches sont instables.
192!-------------------------------------------------------------------------
193
194      active(:)=ztv(:,1)>ztv(:,2)
195      d_temp(:)=0. ! Pour activer un contraste de temperature a la base
196                   ! du panache
197!  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
198      CALL thermcell_alim(thermals_flag_alim,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim)
199
200!------------------------------------------------------------------------------
201! Calcul dans la premiere couche
202! On decide dans cette version que le thermique n'est actif que si la premiere
203! couche est instable.
204! Pourrait etre change si on veut que le thermiques puisse se d??clencher
205! dans une couche l>1
206!------------------------------------------------------------------------------
207do ig=1,ngrid
208! Le panache va prendre au debut les caracteristiques de l'air contenu
209! dans cette couche.
210    if (active(ig)) then
211    ztla(ig,1)=zthl(ig,1)
212    zqta(ig,1)=po(ig,1)
213    zqla(ig,1)=zl(ig,1)
214!cr: attention, prise en compte de f*(1)=1
215    f_star(ig,2)=alim_star(ig,1)
216    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
217&                     *(zlev(ig,2)-zlev(ig,1))  &
218&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
219    w_est(ig,2)=zw2(ig,2)
220    endif
221enddo
222!
223
224!==============================================================================
225!boucle de calcul de la vitesse verticale dans le thermique
226!==============================================================================
227do l=2,klev-1
228!==============================================================================
229
230
231! On decide si le thermique est encore actif ou non
232! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
233    do ig=1,ngrid
234       active(ig)=active(ig) &
235&                 .and. zw2(ig,l)>1.e-10 &
236&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
237    enddo
238
239
240
241!---------------------------------------------------------------------------
242! calcul des proprietes thermodynamiques et de la vitesse de la couche l
243! sans tenir compte du detrainement et de l'entrainement dans cette
244! couche
245! C'est a dire qu'on suppose
246! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
247! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
248! avant) a l'alimentation pour avoir un calcul plus propre
249!---------------------------------------------------------------------------
250
251   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
252   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
253    do ig=1,ngrid
254!       print*,'active',active(ig),ig,l
255        if(active(ig)) then
256        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
257        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
258        zta_est(ig,l)=ztva_est(ig,l)
259        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
260        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
261     &      -zqla_est(ig,l))-zqla_est(ig,l))
262 
263
264!Modif AJAM
265
266        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
267        zdz=zlev(ig,l+1)-zlev(ig,l)         
268        lmel=fact_thermals_ed_dz*zlev(ig,l)
269!        lmel=0.09*zlev(ig,l)
270        zlmel=zlev(ig,l)+lmel
271        zlmelup=zlmel+(zdz/2)
272        zlmeldwn=zlmel-(zdz/2)
273
274        lt=l+1
275        zlt=zlev(ig,lt)
276        zdz3=zlev(ig,lt+1)-zlt
277        zltdwn=zlt-zdz3/2
278        zltup=zlt+zdz3/2
279         
280!=========================================================================
281! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
282!=========================================================================
283
284!--------------------------------------------------
285        if (iflag_thermals_ed.lt.8) then
286!--------------------------------------------------
287!AJ052014: J'ai remplac?? la boucle do par un do while
288! afin de faire moins de calcul dans la boucle
289!--------------------------------------------------
290            do while (zlmelup.gt.zltup)
291               lt=lt+1
292               zlt=zlev(ig,lt)
293               zdz3=zlev(ig,lt+1)-zlt
294               zltdwn=zlt-zdz3/2
295               zltup=zlt+zdz3/2       
296            enddo
297!--------------------------------------------------
298!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
299! on cherche o?? se trouve l'altitude d'inversion
300! en calculant ztv1 (interpolation de la valeur de
301! theta au niveau lt en utilisant les niveaux lt-1 et
302! lt-2) et ztv2 (interpolation avec les niveaux lt+1
303! et lt+2). Si theta r??ellement calcul??e au niveau lt
304! comprise entre ztv1 et ztv2, alors il y a inversion
305! et on calcule son altitude zinv en supposant que ztv(lt)
306! est une combinaison lineaire de ztv1 et ztv2.
307! Ensuite, on calcule la flottabilite en comparant
308! la temperature de la couche l a celle de l'air situe
309! l+lmel plus haut, ce qui necessite de savoir quel fraction
310! de cet air est au-dessus ou en-dessous de l'inversion   
311!--------------------------------------------------
312            atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
313            btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
314    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
315            atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
316            btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
317    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
318
319             ztv1=atv1*zlt+btv1
320             ztv2=atv2*zlt+btv2
321
322             if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 
323
324!--------------------------------------------------
325!AJ052014: D??calage de zinv qui est entre le haut
326!          et le bas de la couche lt
327!--------------------------------------------------
328                factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
329                zinv=zltdwn+zdz3*factinv
330
331         
332                if (zlmeldwn.ge.zinv) then
333                   ztv_est(ig,l)=atv2*zlmel+btv2
334                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
335    &                    +(1.-fact_shell)*zbuoy(ig,l)
336                elseif (zlmelup.ge.zinv) then
337                 ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
338                   ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1
339                   ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1
340
341                   zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
342    &            ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
343    &            ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
344
345                else
346                   ztv_est(ig,l)=atv1*zlmel+btv1
347                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
348    &                           +(1.-fact_shell)*zbuoy(ig,l)
349                endif
350
351             else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
352
353                if (zlmeldwn.gt.zltdwn) then
354                   zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- &
355    &                ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
356                else
357                   zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
358    &                ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
359    &                ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
360
361                endif
362
363!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
364!    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
365!    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
366!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
367!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
368!     &          po(ig,lt-1))/po(ig,lt-1))
369          endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
370
371        else  !   if (iflag_thermals_ed.lt.8) then
372           lt=l+1
373           zlt=zlev(ig,lt)
374           zdz2=zlev(ig,lt)-zlev(ig,l)
375
376           do while (lmel.gt.zdz2)
377             lt=lt+1
378             zlt=zlev(ig,lt)
379             zdz2=zlt-zlev(ig,l)
380           enddo
381           zdz3=zlev(ig,lt+1)-zlt
382           zltdwn=zlev(ig,lt)-zdz3/2
383           zlmelup=zlmel+(zdz/2)
384           coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
385           zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- &
386    &          ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
387    &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
388        endif !   if (iflag_thermals_ed.lt.8) then
389
390!------------------------------------------------
391!AJAM:nouveau calcul de w? 
392!------------------------------------------------
393              zdz=zlev(ig,l+1)-zlev(ig,l)
394              zdzbis=zlev(ig,l)-zlev(ig,l-1)
395              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
396
397              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
398              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
399              zdw2=afact*zbuoy(ig,l)/fact_epsilon
400              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
401!              zdw2bis=0.5*(zdw2+zdw2bis)
402              lm=Max(1,l-2)
403!              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
404!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
405!              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
406!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
407!             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
408!             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
409!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
410!    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
411!              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
412
413!--------------------------------------------------
414!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
415!--------------------------------------------------
416         if (iflag_thermals_ed==8) then
417! Ancienne version
418!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
419!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
420!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
421
422            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
423
424! Nouvelle version Arnaud
425         else
426!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
427!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
428!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
429
430            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
431
432!             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
433!    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
434!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
435
436
437
438!            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
439
440!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
441!    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
442
443!             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
444
445         endif
446
447
448         if (iflag_thermals_ed<6) then
449             zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
450!              fact_epsilon=0.0005/(zalpha+0.025)**0.5
451!              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
452              fact_epsilon=0.0002/(zalpha+0.1)
453              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
454              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
455              zdw2=afact*zbuoy(ig,l)/fact_epsilon
456              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
457!              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
458
459!              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
460!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
461!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
462
463            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
464
465
466         endif
467!--------------------------------------------------
468!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
469!on fait max(0.0001,.....)
470!--------------------------------------------------         
471
472!             if (w_est(ig,l+1).lt.0.) then
473!               w_est(ig,l+1)=zw2(ig,l)
474!                w_est(ig,l+1)=0.0001
475!             endif
476
477       endif
478    enddo
479
480
481!-------------------------------------------------
482!calcul des taux d'entrainement et de detrainement
483!-------------------------------------------------
484
485     do ig=1,ngrid
486        if (active(ig)) then
487
488!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
489          zw2m=w_est(ig,l+1)
490!          zw2m=zw2(ig,l)
491          zdz=zlev(ig,l+1)-zlev(ig,l)
492          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
493!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
494          zbuoybis=zbuoy(ig,l)
495          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
496          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
497
498         
499!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
500!    &     afact*zbuoybis/zw2m - fact_epsilon )
501
502!          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
503!    &     afact*zbuoybis/zw2m - fact_epsilon )
504
505
506
507!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
508
509!=========================================================================
510! 4. Calcul de l'entrainement et du detrainement
511!=========================================================================
512
513!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
514!    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
515!          entrbis=entr_star(ig,l)
516
517          if (iflag_thermals_ed.lt.6) then
518          fact_epsilon=0.0002/(zalpha+0.1)
519          endif
520         
521
522
523          detr_star(ig,l)=f_star(ig,l)*zdz             &
524    &     *( mix0 * 0.1 / (zalpha+0.001)               &
525    &     + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
526    &     + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
527
528!         detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
529!    &                    ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
530
531          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
532
533          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
534    &       mix0 * 0.1 / (zalpha+0.001)               &
535    &     + zbetalpha*MAX(entr_min,                   &
536    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
537
538
539!          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
540!    &       mix0 * 0.1 / (zalpha+0.001)               &
541!    &     + MAX(entr_min,                   &
542!    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
543!    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
544
545
546!         entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
547!    &                    ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
548
549!          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
550!    &     afact*zbuoy(ig,l)/zw2m &
551!    &     - 1.*fact_epsilon)
552
553         
554! En dessous de lalim, on prend le max de alim_star et entr_star pour
555! alim_star et 0 sinon
556        if (l.lt.lalim(ig)) then
557          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
558          entr_star(ig,l)=0.
559        endif
560!        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
561!          alim_star(ig,l)=entrbis
562!        endif
563
564!        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
565! Calcul du flux montant normalise
566      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
567     &              -detr_star(ig,l)
568
569      endif
570   enddo
571
572
573!============================================================================
574! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
575!===========================================================================
576
577   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
578   do ig=1,ngrid
579       if (activetmp(ig)) then
580           Zsat=.false.
581           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
582     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
583     &            /(f_star(ig,l+1)+detr_star(ig,l))
584           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
585     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
586     &            /(f_star(ig,l+1)+detr_star(ig,l))
587
588        endif
589    enddo
590
591   ztemp(:)=zpspsk(:,l)*ztla(:,l)
592   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
593   do ig=1,ngrid
594      if (activetmp(ig)) then
595! on ecrit de maniere conservative (sat ou non)
596!          T = Tl +Lv/Cp ql
597           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
598           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
599           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
600!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
601           zha(ig,l) = ztva(ig,l)
602           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
603     &              -zqla(ig,l))-zqla(ig,l))
604           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
605           zdz=zlev(ig,l+1)-zlev(ig,l)
606           zdzbis=zlev(ig,l)-zlev(ig,l-1)
607           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
608!!!!!!!          fact_epsilon=0.002
609            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
610            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
611            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
612            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
613!              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
614!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
615!              lm=Max(1,l-2)
616!              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
617!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
618!            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
619!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
620!     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
621!            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
622!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
623!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
624!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
625            if (iflag_thermals_ed==8) then
626            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
627            else
628            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
629            endif
630!            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
631!    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
632!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
633
634
635           if (iflag_thermals_ed.lt.6) then
636           zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
637!           fact_epsilon=0.0005/(zalpha+0.025)**0.5
638!           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
639           fact_epsilon=0.0002/(zalpha+0.1)**1
640            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
641            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
642            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
643            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
644
645!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
646!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
647!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
648!            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
649            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
650
651           endif
652
653
654      endif
655   enddo
656
657        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
658!
659!===========================================================================
660! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
661!===========================================================================
662
663   nbpb=0
664   do ig=1,ngrid
665            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
666!               stop'On tombe sur le cas particulier de thermcell_dry'
667!               print*,'On tombe sur le cas particulier de thermcell_plume'
668                nbpb=nbpb+1
669                zw2(ig,l+1)=0.
670                linter(ig)=l+1
671            endif
672
673        if (zw2(ig,l+1).lt.0.) then
674           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
675     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
676           zw2(ig,l+1)=0.
677!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
678        elseif (f_star(ig,l+1).lt.0.) then
679           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
680     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
681           zw2(ig,l+1)=0.
682!fin CR:04/05/12
683        endif
684
685           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
686
687        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
688!   lmix est le niveau de la couche ou w (wa_moy) est maximum
689!on rajoute le calcul de lmix_bis
690            if (zqla(ig,l).lt.1.e-10) then
691               lmix_bis(ig)=l+1
692            endif
693            lmix(ig)=l+1
694            wmaxa(ig)=wa_moy(ig,l+1)
695        endif
696   enddo
697
698   if (nbpb>0) then
699   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
700   endif
701
702!=========================================================================
703! FIN DE LA BOUCLE VERTICALE
704      enddo
705!=========================================================================
706
707!on recalcule alim_star_tot
708       do ig=1,ngrid
709          alim_star_tot(ig)=0.
710       enddo
711       do ig=1,ngrid
712          do l=1,lalim(ig)-1
713          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
714          enddo
715       enddo
716       
717
718        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
719
720#undef wrgrads_thermcell
721#ifdef wrgrads_thermcell
722         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
723         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
724         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
725         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
726         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
727         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
728         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
729#endif
730
731
732     return
733     end
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
770!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
771!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
772!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
773!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
774!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
775 SUBROUTINE thermcell_plume_5B(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
776&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
777&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
778&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
779&           ,lev_out,lunout1,igout)
780!&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
781
782!--------------------------------------------------------------------------
783!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
784! Version conforme a l'article de Rio et al. 2010.
785! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
786!--------------------------------------------------------------------------
787
788      USE print_control_mod, ONLY: prt_level
789      IMPLICIT NONE
790
791#include "YOMCST.h"
792#include "YOETHF.h"
793#include "FCTTRE.h"
794#include "thermcell.h"
795
796      INTEGER itap
797      INTEGER lunout1,igout
798      INTEGER ngrid,klev
799      REAL ptimestep
800      REAL ztv(ngrid,klev)
801      REAL zthl(ngrid,klev)
802      REAL po(ngrid,klev)
803      REAL zl(ngrid,klev)
804      REAL rhobarz(ngrid,klev)
805      REAL zlev(ngrid,klev+1)
806      REAL pplev(ngrid,klev+1)
807      REAL pphi(ngrid,klev)
808      REAL zpspsk(ngrid,klev)
809      REAL alim_star(ngrid,klev)
810      REAL f0(ngrid)
811      INTEGER lalim(ngrid)
812      integer lev_out                           ! niveau pour les print
813      integer nbpb
814   
815      real alim_star_tot(ngrid)
816
817      REAL ztva(ngrid,klev)
818      REAL ztla(ngrid,klev)
819      REAL zqla(ngrid,klev)
820      REAL zqta(ngrid,klev)
821      REAL zha(ngrid,klev)
822
823      REAL detr_star(ngrid,klev)
824      REAL coefc
825      REAL entr_star(ngrid,klev)
826      REAL detr(ngrid,klev)
827      REAL entr(ngrid,klev)
828
829      REAL csc(ngrid,klev)
830
831      REAL zw2(ngrid,klev+1)
832      REAL w_est(ngrid,klev+1)
833      REAL f_star(ngrid,klev+1)
834      REAL wa_moy(ngrid,klev+1)
835
836      REAL ztva_est(ngrid,klev)
837      REAL zqla_est(ngrid,klev)
838      REAL zqsatth(ngrid,klev)
839      REAL zta_est(ngrid,klev)
840      REAL zbuoyjam(ngrid,klev)
841      REAL ztemp(ngrid),zqsat(ngrid)
842      REAL zdw2
843      REAL zw2modif
844      REAL zw2fact
845      REAL zeps(ngrid,klev)
846
847      REAL linter(ngrid)
848      INTEGER lmix(ngrid)
849      INTEGER lmix_bis(ngrid)
850      REAL    wmaxa(ngrid)
851
852      INTEGER ig,l,k
853
854      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
855      real zbuoybis
856      real zcor,zdelta,zcvm5,qlbef,zdz2
857      real betalpha,zbetalpha
858      real eps, afact
859      REAL REPS,RLvCp,DDT0
860      PARAMETER (DDT0=.01)
861      logical Zsat
862      LOGICAL active(ngrid),activetmp(ngrid)
863      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
864      REAL c2(ngrid,klev)
865      Zsat=.false.
866! Initialisation
867
868      RLvCp = RLVTT/RCPD
869      fact_epsilon=0.002
870      betalpha=0.9
871      afact=2./3.           
872
873      zbetalpha=betalpha/(1.+betalpha)
874
875
876! Initialisations des variables reeles
877if (1==1) then
878      ztva(:,:)=ztv(:,:)
879      ztva_est(:,:)=ztva(:,:)
880      ztla(:,:)=zthl(:,:)
881      zqta(:,:)=po(:,:)
882      zha(:,:) = ztva(:,:)
883else
884      ztva(:,:)=0.
885      ztva_est(:,:)=0.
886      ztla(:,:)=0.
887      zqta(:,:)=0.
888      zha(:,:) =0.
889endif
890
891      zqla_est(:,:)=0.
892      zqsatth(:,:)=0.
893      zqla(:,:)=0.
894      detr_star(:,:)=0.
895      entr_star(:,:)=0.
896      alim_star(:,:)=0.
897      alim_star_tot(:)=0.
898      csc(:,:)=0.
899      detr(:,:)=0.
900      entr(:,:)=0.
901      zw2(:,:)=0.
902      zbuoy(:,:)=0.
903      zbuoyjam(:,:)=0.
904      gamma(:,:)=0.
905      zeps(:,:)=0.
906      w_est(:,:)=0.
907      f_star(:,:)=0.
908      wa_moy(:,:)=0.
909      linter(:)=1.
910!     linter(:)=1.
911! Initialisation des variables entieres
912      lmix(:)=1
913      lmix_bis(:)=2
914      wmaxa(:)=0.
915      lalim(:)=1
916
917
918!-------------------------------------------------------------------------
919! On ne considere comme actif que les colonnes dont les deux premieres
920! couches sont instables.
921!-------------------------------------------------------------------------
922      active(:)=ztv(:,1)>ztv(:,2)
923
924!-------------------------------------------------------------------------
925! Definition de l'alimentation a l'origine dans thermcell_init
926!-------------------------------------------------------------------------
927      do l=1,klev-1
928         do ig=1,ngrid
929            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
930               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
931     &                       *sqrt(zlev(ig,l+1))
932               lalim(ig)=l+1
933               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
934            endif
935         enddo
936      enddo
937      do l=1,klev
938         do ig=1,ngrid
939            if (alim_star_tot(ig) > 1.e-10 ) then
940               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
941            endif
942         enddo
943      enddo
944      alim_star_tot(:)=1.
945
946
947
948!------------------------------------------------------------------------------
949! Calcul dans la premiere couche
950! On decide dans cette version que le thermique n'est actif que si la premiere
951! couche est instable.
952! Pourrait etre change si on veut que le thermiques puisse se d??clencher
953! dans une couche l>1
954!------------------------------------------------------------------------------
955do ig=1,ngrid
956! Le panache va prendre au debut les caracteristiques de l'air contenu
957! dans cette couche.
958    if (active(ig)) then
959    ztla(ig,1)=zthl(ig,1)
960    zqta(ig,1)=po(ig,1)
961    zqla(ig,1)=zl(ig,1)
962!cr: attention, prise en compte de f*(1)=1
963    f_star(ig,2)=alim_star(ig,1)
964    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
965&                     *(zlev(ig,2)-zlev(ig,1))  &
966&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
967    w_est(ig,2)=zw2(ig,2)
968    endif
969enddo
970!
971
972!==============================================================================
973!boucle de calcul de la vitesse verticale dans le thermique
974!==============================================================================
975do l=2,klev-1
976!==============================================================================
977
978
979! On decide si le thermique est encore actif ou non
980! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
981    do ig=1,ngrid
982       active(ig)=active(ig) &
983&                 .and. zw2(ig,l)>1.e-10 &
984&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
985    enddo
986
987
988
989!---------------------------------------------------------------------------
990! calcul des proprietes thermodynamiques et de la vitesse de la couche l
991! sans tenir compte du detrainement et de l'entrainement dans cette
992! couche
993! C'est a dire qu'on suppose
994! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
995! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
996! avant) a l'alimentation pour avoir un calcul plus propre
997!---------------------------------------------------------------------------
998
999   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
1000   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
1001
1002    do ig=1,ngrid
1003!       print*,'active',active(ig),ig,l
1004        if(active(ig)) then
1005        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
1006        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
1007        zta_est(ig,l)=ztva_est(ig,l)
1008        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
1009        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
1010     &      -zqla_est(ig,l))-zqla_est(ig,l))
1011
1012!------------------------------------------------
1013!AJAM:nouveau calcul de w? 
1014!------------------------------------------------
1015              zdz=zlev(ig,l+1)-zlev(ig,l)
1016              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1017
1018              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1019              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
1020              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
1021 
1022
1023             if (w_est(ig,l+1).lt.0.) then
1024                w_est(ig,l+1)=zw2(ig,l)
1025             endif
1026       endif
1027    enddo
1028
1029
1030!-------------------------------------------------
1031!calcul des taux d'entrainement et de detrainement
1032!-------------------------------------------------
1033
1034     do ig=1,ngrid
1035        if (active(ig)) then
1036
1037          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
1038          zw2m=w_est(ig,l+1)
1039          zdz=zlev(ig,l+1)-zlev(ig,l)
1040          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1041!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
1042          zbuoybis=zbuoy(ig,l)
1043          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
1044          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
1045
1046         
1047          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
1048    &     afact*zbuoybis/zw2m - fact_epsilon )
1049
1050
1051          detr_star(ig,l)=f_star(ig,l)*zdz                        &
1052    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
1053    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
1054         
1055! En dessous de lalim, on prend le max de alim_star et entr_star pour
1056! alim_star et 0 sinon
1057        if (l.lt.lalim(ig)) then
1058          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
1059          entr_star(ig,l)=0.
1060        endif
1061
1062! Calcul du flux montant normalise
1063      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
1064     &              -detr_star(ig,l)
1065
1066      endif
1067   enddo
1068
1069
1070!----------------------------------------------------------------------------
1071!calcul de la vitesse verticale en melangeant Tl et qt du thermique
1072!---------------------------------------------------------------------------
1073   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
1074   do ig=1,ngrid
1075       if (activetmp(ig)) then
1076           Zsat=.false.
1077           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
1078     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
1079     &            /(f_star(ig,l+1)+detr_star(ig,l))
1080           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
1081     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
1082     &            /(f_star(ig,l+1)+detr_star(ig,l))
1083
1084        endif
1085    enddo
1086
1087   ztemp(:)=zpspsk(:,l)*ztla(:,l)
1088   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
1089
1090   do ig=1,ngrid
1091      if (activetmp(ig)) then
1092! on ecrit de maniere conservative (sat ou non)
1093!          T = Tl +Lv/Cp ql
1094           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
1095           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1096           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1097!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1098           zha(ig,l) = ztva(ig,l)
1099           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
1100     &              -zqla(ig,l))-zqla(ig,l))
1101           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1102           zdz=zlev(ig,l+1)-zlev(ig,l)
1103           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
1104
1105            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1106            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
1107            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
1108      endif
1109   enddo
1110
1111        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1112!
1113!---------------------------------------------------------------------------
1114!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1115!---------------------------------------------------------------------------
1116
1117   nbpb=0
1118   do ig=1,ngrid
1119            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1120!               stop'On tombe sur le cas particulier de thermcell_dry'
1121!               print*,'On tombe sur le cas particulier de thermcell_plume'
1122                nbpb=nbpb+1
1123                zw2(ig,l+1)=0.
1124                linter(ig)=l+1
1125            endif
1126
1127        if (zw2(ig,l+1).lt.0.) then
1128           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
1129     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1130           zw2(ig,l+1)=0.
1131        elseif (f_star(ig,l+1).lt.0.) then
1132           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
1133     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
1134!           print*,"linter plume", linter(ig)
1135           zw2(ig,l+1)=0.
1136        endif
1137
1138           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1139
1140        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1141!   lmix est le niveau de la couche ou w (wa_moy) est maximum
1142!on rajoute le calcul de lmix_bis
1143            if (zqla(ig,l).lt.1.e-10) then
1144               lmix_bis(ig)=l+1
1145            endif
1146            lmix(ig)=l+1
1147            wmaxa(ig)=wa_moy(ig,l+1)
1148        endif
1149   enddo
1150
1151   if (nbpb>0) then
1152   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
1153   endif
1154
1155!=========================================================================
1156! FIN DE LA BOUCLE VERTICALE
1157      enddo
1158!=========================================================================
1159
1160!on recalcule alim_star_tot
1161       do ig=1,ngrid
1162          alim_star_tot(ig)=0.
1163       enddo
1164       do ig=1,ngrid
1165          do l=1,lalim(ig)-1
1166          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1167          enddo
1168       enddo
1169       
1170
1171        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1172
1173#undef wrgrads_thermcell
1174#ifdef wrgrads_thermcell
1175         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
1176         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
1177         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
1178         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
1179         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
1180         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
1181         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
1182#endif
1183
1184
1185     return
1186     end
Note: See TracBrowser for help on using the repository browser.