source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/lmdz_thermcell_plume_6A.F90 @ 5097

Last change on this file since 5097 was 4727, checked in by idelkadi, 13 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch
  • 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 
1MODULE lmdz_thermcell_plume_6A
2!
3! $Id: lmdz_thermcell_plume_6A.F90 4727 2023-10-19 14:02:57Z snguyen $
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(ngrid)
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#undef wrgrads_thermcell
687#ifdef wrgrads_thermcell
688         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
689         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
690         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
691         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
692         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
693         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
694         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
695#endif
696
697
698 RETURN
699     end
700
701
702
703
704
705
706
707!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
708!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
709!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
710!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
711!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
712!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
713 SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
714&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
715&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
716&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
717&           ,lev_out,lunout1,igout)
718!&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
719
720!--------------------------------------------------------------------------
721!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
722! Version conforme a l'article de Rio et al. 2010.
723! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
724!--------------------------------------------------------------------------
725
726      USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
727       USE lmdz_thermcell_qsat, ONLY : thermcell_qsat
728      IMPLICIT NONE
729
730      INTEGER itap
731      INTEGER lunout1,igout
732      INTEGER ngrid,nlay
733      REAL ptimestep
734      REAL ztv(ngrid,nlay)
735      REAL zthl(ngrid,nlay)
736      REAL, INTENT(IN) :: po(ngrid,nlay)
737      REAL zl(ngrid,nlay)
738      REAL rhobarz(ngrid,nlay)
739      REAL zlev(ngrid,nlay+1)
740      REAL pplev(ngrid,nlay+1)
741      REAL pphi(ngrid,nlay)
742      REAL zpspsk(ngrid,nlay)
743      REAL alim_star(ngrid,nlay)
744      REAL f0(ngrid)
745      INTEGER lalim(ngrid)
746      integer lev_out                           ! niveau pour les print
747      integer nbpb
748   
749      real alim_star_tot(ngrid)
750
751      REAL ztva(ngrid,nlay)
752      REAL ztla(ngrid,nlay)
753      REAL zqla(ngrid,nlay)
754      REAL zqta(ngrid,nlay)
755      REAL zha(ngrid,nlay)
756
757      REAL detr_star(ngrid,nlay)
758      REAL coefc
759      REAL entr_star(ngrid,nlay)
760      REAL detr(ngrid,nlay)
761      REAL entr(ngrid,nlay)
762
763      REAL csc(ngrid,nlay)
764
765      REAL zw2(ngrid,nlay+1)
766      REAL w_est(ngrid,nlay+1)
767      REAL f_star(ngrid,nlay+1)
768      REAL wa_moy(ngrid,nlay+1)
769
770      REAL ztva_est(ngrid,nlay)
771      REAL zqla_est(ngrid,nlay)
772      REAL zqsatth(ngrid,nlay)
773      REAL zta_est(ngrid,nlay)
774      REAL zbuoyjam(ngrid,nlay)
775      REAL ztemp(ngrid),zqsat(ngrid)
776      REAL zdw2
777      REAL zw2modif
778      REAL zw2fact
779      REAL zeps(ngrid,nlay)
780
781      REAL linter(ngrid)
782      INTEGER lmix(ngrid)
783      INTEGER lmix_bis(ngrid)
784      REAL    wmaxa(ngrid)
785
786      INTEGER ig,l,k
787
788      real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
789      real zbuoybis
790      real zcor,zdelta,zcvm5,qlbef,zdz2
791      real betalpha,zbetalpha
792      real eps, afact
793      logical Zsat
794      LOGICAL active(ngrid),activetmp(ngrid)
795      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
796      REAL c2(ngrid,nlay)
797      Zsat=.false.
798! Initialisation
799
800      fact_epsilon=0.002
801      betalpha=0.9
802      afact=2./3.           
803
804      zbetalpha=betalpha/(1.+betalpha)
805
806
807! Initialisations des variables reeles
808if (1==1) then
809      ztva(:,:)=ztv(:,:)
810      ztva_est(:,:)=ztva(:,:)
811      ztla(:,:)=zthl(:,:)
812      zqta(:,:)=po(:,:)
813      zha(:,:) = ztva(:,:)
814else
815      ztva(:,:)=0.
816      ztva_est(:,:)=0.
817      ztla(:,:)=0.
818      zqta(:,:)=0.
819      zha(:,:) =0.
820endif
821
822      zqla_est(:,:)=0.
823      zqsatth(:,:)=0.
824      zqla(:,:)=0.
825      detr_star(:,:)=0.
826      entr_star(:,:)=0.
827      alim_star(:,:)=0.
828      alim_star_tot(:)=0.
829      csc(:,:)=0.
830      detr(:,:)=0.
831      entr(:,:)=0.
832      zw2(:,:)=0.
833      zbuoy(:,:)=0.
834      zbuoyjam(:,:)=0.
835      gamma(:,:)=0.
836      zeps(:,:)=0.
837      w_est(:,:)=0.
838      f_star(:,:)=0.
839      wa_moy(:,:)=0.
840      linter(:)=1.
841!     linter(:)=1.
842! Initialisation des variables entieres
843      lmix(:)=1
844      lmix_bis(:)=2
845      wmaxa(:)=0.
846      lalim(:)=1
847
848
849!-------------------------------------------------------------------------
850! On ne considere comme actif que les colonnes dont les deux premieres
851! couches sont instables.
852!-------------------------------------------------------------------------
853      active(:)=ztv(:,1)>ztv(:,2)
854
855!-------------------------------------------------------------------------
856! Definition de l'alimentation
857!-------------------------------------------------------------------------
858      do l=1,nlay-1
859         do ig=1,ngrid
860            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
861               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
862     &                       *sqrt(zlev(ig,l+1))
863               lalim(ig)=l+1
864               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
865            endif
866         enddo
867      enddo
868      do l=1,nlay
869         do ig=1,ngrid
870            if (alim_star_tot(ig) > 1.e-10 ) then
871               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
872            endif
873         enddo
874      enddo
875      alim_star_tot(:)=1.
876
877
878
879!------------------------------------------------------------------------------
880! Calcul dans la premiere couche
881! On decide dans cette version que le thermique n'est actif que si la premiere
882! couche est instable.
883! Pourrait etre change si on veut que le thermiques puisse se d??clencher
884! dans une couche l>1
885!------------------------------------------------------------------------------
886do ig=1,ngrid
887! Le panache va prendre au debut les caracteristiques de l'air contenu
888! dans cette couche.
889    if (active(ig)) then
890    ztla(ig,1)=zthl(ig,1)
891    zqta(ig,1)=po(ig,1)
892    zqla(ig,1)=zl(ig,1)
893!cr: attention, prise en compte de f*(1)=1
894    f_star(ig,2)=alim_star(ig,1)
895    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
896&                     *(zlev(ig,2)-zlev(ig,1))  &
897&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
898    w_est(ig,2)=zw2(ig,2)
899    endif
900enddo
901!
902
903!==============================================================================
904!boucle de calcul de la vitesse verticale dans le thermique
905!==============================================================================
906do l=2,nlay-1
907!==============================================================================
908
909
910! On decide si le thermique est encore actif ou non
911! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
912    do ig=1,ngrid
913       active(ig)=active(ig) &
914&                 .and. zw2(ig,l)>1.e-10 &
915&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
916    enddo
917
918
919
920!---------------------------------------------------------------------------
921! calcul des proprietes thermodynamiques et de la vitesse de la couche l
922! sans tenir compte du detrainement et de l'entrainement dans cette
923! couche
924! C'est a dire qu'on suppose
925! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
926! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
927! avant) a l'alimentation pour avoir un calcul plus propre
928!---------------------------------------------------------------------------
929
930   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
931   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
932
933    do ig=1,ngrid
934!       print*,'active',active(ig),ig,l
935        if(active(ig)) then
936        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
937        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
938        zta_est(ig,l)=ztva_est(ig,l)
939        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
940        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
941     &      -zqla_est(ig,l))-zqla_est(ig,l))
942
943!------------------------------------------------
944!AJAM:nouveau calcul de w? 
945!------------------------------------------------
946              zdz=zlev(ig,l+1)-zlev(ig,l)
947              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
948
949              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
950              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
951              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
952 
953
954             if (w_est(ig,l+1).lt.0.) then
955                w_est(ig,l+1)=zw2(ig,l)
956             endif
957       endif
958    enddo
959
960
961!-------------------------------------------------
962!calcul des taux d'entrainement et de detrainement
963!-------------------------------------------------
964
965     do ig=1,ngrid
966        if (active(ig)) then
967
968          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
969          zw2m=w_est(ig,l+1)
970          zdz=zlev(ig,l+1)-zlev(ig,l)
971          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
972!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
973          zbuoybis=zbuoy(ig,l)
974          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
975          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
976
977         
978          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
979    &     afact*zbuoybis/zw2m - fact_epsilon )
980
981
982          detr_star(ig,l)=f_star(ig,l)*zdz                        &
983    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
984    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
985         
986! En dessous de lalim, on prend le max de alim_star et entr_star pour
987! alim_star et 0 sinon
988        if (l.lt.lalim(ig)) then
989          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
990          entr_star(ig,l)=0.
991        endif
992
993! Calcul du flux montant normalise
994      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
995     &              -detr_star(ig,l)
996
997      endif
998   enddo
999
1000
1001!----------------------------------------------------------------------------
1002!calcul de la vitesse verticale en melangeant Tl et qt du thermique
1003!---------------------------------------------------------------------------
1004   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
1005   do ig=1,ngrid
1006       if (activetmp(ig)) then
1007           Zsat=.false.
1008           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
1009     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
1010     &            /(f_star(ig,l+1)+detr_star(ig,l))
1011           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
1012     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
1013     &            /(f_star(ig,l+1)+detr_star(ig,l))
1014
1015        endif
1016    enddo
1017
1018   ztemp(:)=zpspsk(:,l)*ztla(:,l)
1019   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
1020
1021   do ig=1,ngrid
1022      if (activetmp(ig)) then
1023! on ecrit de maniere conservative (sat ou non)
1024!          T = Tl +Lv/Cp ql
1025           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
1026           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1027           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1028!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1029           zha(ig,l) = ztva(ig,l)
1030           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
1031     &              -zqla(ig,l))-zqla(ig,l))
1032           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1033           zdz=zlev(ig,l+1)-zlev(ig,l)
1034           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
1035
1036            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1037            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
1038            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
1039      endif
1040   enddo
1041
1042        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1043!
1044!---------------------------------------------------------------------------
1045!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1046!---------------------------------------------------------------------------
1047
1048   nbpb=0
1049   do ig=1,ngrid
1050            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1051!               stop'On tombe sur le cas particulier de thermcell_dry'
1052!               print*,'On tombe sur le cas particulier de thermcell_plume'
1053                nbpb=nbpb+1
1054                zw2(ig,l+1)=0.
1055                linter(ig)=l+1
1056            endif
1057
1058        if (zw2(ig,l+1).lt.0.) then
1059           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
1060     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1061           zw2(ig,l+1)=0.
1062        elseif (f_star(ig,l+1).lt.0.) then
1063           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
1064     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
1065!           print*,"linter plume", linter(ig)
1066           zw2(ig,l+1)=0.
1067        endif
1068
1069           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1070
1071        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1072!   lmix est le niveau de la couche ou w (wa_moy) est maximum
1073!on rajoute le calcul de lmix_bis
1074            if (zqla(ig,l).lt.1.e-10) then
1075               lmix_bis(ig)=l+1
1076            endif
1077            lmix(ig)=l+1
1078            wmaxa(ig)=wa_moy(ig,l+1)
1079        endif
1080   enddo
1081
1082   if (nbpb>0) then
1083   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
1084   endif
1085
1086!=========================================================================
1087! FIN DE LA BOUCLE VERTICALE
1088      enddo
1089!=========================================================================
1090
1091!on recalcule alim_star_tot
1092       do ig=1,ngrid
1093          alim_star_tot(ig)=0.
1094       enddo
1095       do ig=1,ngrid
1096          do l=1,lalim(ig)-1
1097          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1098          enddo
1099       enddo
1100       
1101
1102        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1103
1104#undef wrgrads_thermcell
1105#ifdef wrgrads_thermcell
1106         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
1107         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
1108         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
1109         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
1110         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
1111         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
1112         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
1113#endif
1114
1115
1116     return
1117     end
1118END MODULE lmdz_thermcell_plume_6A
Note: See TracBrowser for help on using the repository browser.