source: LMDZ6/branches/contrails/libf/phylmd/lmdz_thermcell_plume_6A.f90 @ 5450

Last change on this file since 5450 was 5450, checked in by aborella, 7 hours ago

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