source: LMDZ5/trunk/libf/phylmd/thermcell_plume.F90 @ 2989

Last change on this file since 2989 was 2406, checked in by fhourdin, 9 years ago

Small corrections

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