source: LMDZ5/branches/testing/libf/phylmd/thermcell_plume.F90 @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago

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