source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/thermcell_plume.F90

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

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

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