source: lmdz_wrf/trunk/WRFV3/lmdz/thermcell_plume.F90 @ 1554

Last change on this file since 1554 was 179, checked in by lfita, 10 years ago

No Ale, Alp, ... values because pphi is at the intermediate level!!!
Removing raw prints of checkig for Ale, Alp

File size: 27.7 KB
Line 
1!
2! $Id: thermcell_plume.F90 1503 2011-03-23 11:57:52Z idelkadi $
3!
4      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
5     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
6     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
7     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8     &           ,lev_out,lunout1,igout)
9
10!--------------------------------------------------------------------------
11!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12!--------------------------------------------------------------------------
13
14      IMPLICIT NONE
15
16#include "YOMCST.h"
17#include "YOETHF.h"
18#include "FCTTRE.h"
19#include "iniprint.h"
20#include "thermcell.h"
21
22      INTEGER itap
23      INTEGER lunout1,igout
24      INTEGER ngrid,klev
25      REAL ptimestep
26      REAL ztv(ngrid,klev)
27      REAL zthl(ngrid,klev)
28      REAL po(ngrid,klev)
29      REAL zl(ngrid,klev)
30      REAL rhobarz(ngrid,klev)
31      REAL zlev(ngrid,klev+1)
32      REAL pplev(ngrid,klev+1)
33      REAL pphi(ngrid,klev)
34      REAL zpspsk(ngrid,klev)
35      REAL alim_star(ngrid,klev)
36      REAL f0(ngrid)
37      INTEGER lalim(ngrid)
38      integer lev_out                           ! niveau pour les print
39      integer nbpb
40      real zcon2(ngrid)
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 zqla_est(ngrid,klev)
65      REAL zqsatth(ngrid,klev)
66      REAL zta_est(ngrid,klev)
67      REAL zdw2
68      REAL zw2modif
69      REAL zeps
70
71      REAL linter(ngrid)
72      INTEGER lmix(ngrid)
73      INTEGER lmix_bis(ngrid)
74      REAL    wmaxa(ngrid)
75
76      INTEGER ig,l,k
77
78      real zdz,zfact,zbuoy,zalpha,zdrag
79      real zcor,zdelta,zcvm5,qlbef
80      real Tbef,qsatbef
81      real dqsat_dT,DT,num,denom
82      REAL REPS,RLvCp,DDT0
83      PARAMETER (DDT0=.01)
84      logical Zsat
85      LOGICAL active(ngrid),activetmp(ngrid)
86      REAL fact_gamma,fact_epsilon,fact_gamma2
87      REAL c2(ngrid,klev)
88      REAL a1,m
89
90      REAL zw2fact,expa
91      Zsat=.false.
92
93! Initialisation
94      RLvCp = RLVTT/RCPD
95     
96     
97         fact_epsilon=0.002
98         a1=2./3.
99         fact_gamma=0.9
100         zfact=fact_gamma/(1+fact_gamma)
101         fact_gamma2=zfact
102         expa=0.
103     
104
105! Initialisations des variables reeles
106if (1==1) then
107      ztva(:,:)=ztv(:,:)
108      ztva_est(:,:)=ztva(:,:)
109      ztla(:,:)=zthl(:,:)
110      zqta(:,:)=po(:,:)
111      zha(:,:) = ztva(:,:)
112else
113      ztva(:,:)=0.
114      ztva_est(:,:)=0.
115      ztla(:,:)=0.
116      zqta(:,:)=0.
117      zha(:,:) =0.
118endif
119
120      zqla_est(:,:)=0.
121      zqsatth(:,:)=0.
122      zqla(:,:)=0.
123      detr_star(:,:)=0.
124      entr_star(:,:)=0.
125      alim_star(:,:)=0.
126      alim_star_tot(:)=0.
127      csc(:,:)=0.
128      detr(:,:)=0.
129      entr(:,:)=0.
130      zw2(:,:)=0.
131      w_est(:,:)=0.
132      f_star(:,:)=0.
133      wa_moy(:,:)=0.
134      linter(:)=1.
135      linter(:)=1.
136
137! Initialisation des variables entieres
138      lmix(:)=1
139      lmix_bis(:)=2
140      wmaxa(:)=0.
141      lalim(:)=1
142
143!-------------------------------------------------------------------------
144! On ne considere comme actif que les colonnes dont les deux premieres
145! couches sont instables.
146!-------------------------------------------------------------------------
147      active(:)=ztv(:,1)>ztv(:,2)
148
149!-------------------------------------------------------------------------
150! Definition de l'alimentation a l'origine dans thermcell_init
151!-------------------------------------------------------------------------
152      do l=1,klev-1
153         do ig=1,ngrid
154            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
155               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
156     &                       *sqrt(zlev(ig,l+1))
157               lalim(ig)=l+1
158               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
159            endif
160         enddo
161      enddo
162      do l=1,klev
163         do ig=1,ngrid
164            if (alim_star_tot(ig) > 1.e-10 ) then
165               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
166            endif
167         enddo
168      enddo
169      alim_star_tot(:)=1.
170
171!------------------------------------------------------------------------------
172! Calcul dans la premiere couche
173! On decide dans cette version que le thermique n'est actif que si la premiere
174! couche est instable.
175! Pourrait etre change si on veut que le thermiques puisse se déclencher
176! dans une couche l>1
177!------------------------------------------------------------------------------
178do ig=1,ngrid
179! Le panache va prendre au debut les caracteristiques de l'air contenu
180! dans cette couche.
181    if (active(ig)) then
182    ztla(ig,1)=zthl(ig,1)
183    zqta(ig,1)=po(ig,1)
184    zqla(ig,1)=zl(ig,1)
185!cr: attention, prise en compte de f*(1)=1
186    f_star(ig,2)=alim_star(ig,1)
187    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
188&                     *(zlev(ig,2)-zlev(ig,1))  &
189&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
190    w_est(ig,2)=zw2(ig,2)
191    endif
192enddo
193!
194!==============================================================================
195!boucle de calcul de la vitesse verticale dans le thermique
196!==============================================================================
197do l=2,klev-1
198!==============================================================================
199
200
201! On decide si le thermique est encore actif ou non
202! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
203    do ig=1,ngrid
204       active(ig)=active(ig) &
205&                 .and. zw2(ig,l)>1.e-10 &
206&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
207    enddo
208
209
210
211! Premier calcul de la vitesse verticale a partir de la temperature
212! potentielle virtuelle
213!     if (1.eq.1) then
214!         w_est(ig,3)=zw2(ig,2)* &
215!    &      ((f_star(ig,2))**2) &
216!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
217!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
218!    &      *(zlev(ig,3)-zlev(ig,2))
219!      endif
220
221
222!---------------------------------------------------------------------------
223! calcul des proprietes thermodynamiques et de la vitesse de la couche l
224! sans tenir compte du detrainement et de l'entrainement dans cette
225! couche
226! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
227! avant) a l'alimentation pour avoir un calcul plus propre
228!---------------------------------------------------------------------------
229
230     call thermcell_condens(ngrid,active, &
231&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
232
233    do ig=1,ngrid
234        if(active(ig)) then
235        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
236        zta_est(ig,l)=ztva_est(ig,l)
237        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
238        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
239     &      -zqla_est(ig,l))-zqla_est(ig,l))
240
241         if (1.eq.0) then
242!calcul de w_est sans prendre en compte le drag
243             w_est(ig,l+1)=zw2(ig,l)*  &
244     &                   ((f_star(ig,l))**2)  &
245     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
246     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
247     &                   *(zlev(ig,l+1)-zlev(ig,l))
248         else
249
250           zdz=zlev(ig,l+1)-zlev(ig,l)
251           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
252           zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
253           zdrag=fact_epsilon/(zalpha**expa)
254           zw2fact=zbuoy/zdrag*a1
255           w_est(ig,l+1)=(w_est(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
256      &    +zw2fact
257
258         endif
259
260             if (w_est(ig,l+1).lt.0.) then
261                w_est(ig,l+1)=zw2(ig,l)
262             endif
263       endif
264    enddo
265
266!-------------------------------------------------
267!calcul des taux d'entrainement et de detrainement
268!-------------------------------------------------
269
270     do ig=1,ngrid
271        if (active(ig)) then
272          zdz=zlev(ig,l+1)-zlev(ig,l)
273          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
274 
275! estimation de la fraction couverte par les thermiques
276          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
277
278!calcul de la soumission papier
279! Calcul  du taux d'entrainement entr_star (epsilon)
280           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
281     &     a1*zbuoy/w_est(ig,l+1) &
282     &     - fact_epsilon/zalpha**expa  ) &
283     &     +0. )
284           
285!calcul du taux de detrainment (delta)
286!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
287!     &      MAX(1.e-3, &
288!     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
289!     &      +0.01*(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))**0.5    &   
290!     &     +0. ))
291
292          m=0.5
293
294          detr_star(ig,l)=1.*f_star(ig,l)*zdz *                    &
295    &     MAX(5.e-4,-fact_gamma2*a1*(1./w_est(ig,l+1))*((1.-(1.-m)/(1.+70*zqta(ig,l-1)))*zbuoy        &
296    &     -40*(1.-m)*(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) )   )
297
298!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
299!     &      MAX(0.0, &
300!     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
301!     &      +20*(max(zqta(ig,l-1)-po(ig,l),0.))**1*(zalpha/w_est(ig,l+1))**0.5    &   
302!     &     +0. ))
303
304
305! En dessous de lalim, on prend le max de alim_star et entr_star pour
306! alim_star et 0 sinon
307        if (l.lt.lalim(ig)) then
308          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
309          entr_star(ig,l)=0.
310        endif
311
312!attention test
313!        if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then       
314!            detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
315!        endif
316! Calcul du flux montant normalise
317      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
318     &              -detr_star(ig,l)
319
320      endif
321   enddo
322
323!----------------------------------------------------------------------------
324!calcul de la vitesse verticale en melangeant Tl et qt du thermique
325!---------------------------------------------------------------------------
326   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
327   do ig=1,ngrid
328       if (activetmp(ig)) then
329           Zsat=.false.
330           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
331     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
332     &            /(f_star(ig,l+1)+detr_star(ig,l))
333           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
334     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
335     &            /(f_star(ig,l+1)+detr_star(ig,l))
336
337        endif
338    enddo
339
340   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
341
342
343   do ig=1,ngrid
344      if (activetmp(ig)) then
345! on ecrit de maniere conservative (sat ou non)
346!          T = Tl +Lv/Cp ql
347           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
348           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
349!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
350           zha(ig,l) = ztva(ig,l)
351           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
352     &              -zqla(ig,l))-zqla(ig,l))
353
354!on ecrit zqsat
355           zqsatth(ig,l)=qsatbef 
356
357!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
358!          zw2(ig,l+1)=&
359!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
360!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
361!     &                 *1./(1.+fact_gamma) &
362!     &                 *(zlev(ig,l+1)-zlev(ig,l))
363!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
364! La meme en plus modulaire :
365           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
366           zdz=zlev(ig,l+1)-zlev(ig,l)
367
368
369           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
370
371!if (1==0) then
372!           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
373!           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
374!           zw2(ig,l+1)=zw2modif+zdw2
375!else
376           zdrag=fact_epsilon/(zalpha**expa)
377           zw2fact=zbuoy/zdrag*a1
378           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
379      &    +zw2fact
380!endif
381
382      endif
383   enddo
384
385        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
386!
387!---------------------------------------------------------------------------
388!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
389!---------------------------------------------------------------------------
390
391   nbpb=0
392   do ig=1,ngrid
393            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
394!               stop'On tombe sur le cas particulier de thermcell_dry'
395!               print*,'On tombe sur le cas particulier de thermcell_plume'
396                nbpb=nbpb+1
397                zw2(ig,l+1)=0.
398                linter(ig)=l+1
399            endif
400
401        if (zw2(ig,l+1).lt.0.) then
402           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
403     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
404           zw2(ig,l+1)=0.
405        endif
406
407           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
408
409        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
410!   lmix est le niveau de la couche ou w (wa_moy) est maximum
411!on rajoute le calcul de lmix_bis
412            if (zqla(ig,l).lt.1.e-10) then
413               lmix_bis(ig)=l+1
414            endif
415            lmix(ig)=l+1
416            wmaxa(ig)=wa_moy(ig,l+1)
417        endif
418   enddo
419
420   if (nbpb>0) then
421   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
422   endif
423
424!=========================================================================
425! FIN DE LA BOUCLE VERTICALE
426      enddo
427!=========================================================================
428
429!on recalcule alim_star_tot
430       do ig=1,ngrid
431          alim_star_tot(ig)=0.
432       enddo
433       do ig=1,ngrid
434          do l=1,lalim(ig)-1
435          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
436          enddo
437       enddo
438       
439
440        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
441
442
443      return
444      end
445
446!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452 SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
453&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
454&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
455&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
456&           ,lev_out,lunout1,igout)
457
458!--------------------------------------------------------------------------
459!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
460! Version conforme a l'article de Rio et al. 2010.
461! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
462!--------------------------------------------------------------------------
463
464      IMPLICIT NONE
465
466#include "YOMCST.h"
467#include "YOETHF.h"
468#include "FCTTRE.h"
469#include "iniprint.h"
470#include "thermcell.h"
471
472      INTEGER itap
473      INTEGER lunout1,igout
474      INTEGER ngrid,klev
475      REAL ptimestep
476      REAL ztv(ngrid,klev)
477      REAL zthl(ngrid,klev)
478      REAL po(ngrid,klev)
479      REAL zl(ngrid,klev)
480      REAL rhobarz(ngrid,klev)
481      REAL zlev(ngrid,klev+1)
482      REAL pplev(ngrid,klev+1)
483      REAL pphi(ngrid,klev)
484      REAL zpspsk(ngrid,klev)
485      REAL alim_star(ngrid,klev)
486      REAL f0(ngrid)
487      INTEGER lalim(ngrid)
488      integer lev_out                           ! niveau pour les print
489      integer nbpb
490   
491      real alim_star_tot(ngrid)
492
493      REAL ztva(ngrid,klev)
494      REAL ztla(ngrid,klev)
495      REAL zqla(ngrid,klev)
496      REAL zqta(ngrid,klev)
497      REAL zha(ngrid,klev)
498
499      REAL detr_star(ngrid,klev)
500      REAL coefc
501      REAL entr_star(ngrid,klev)
502      REAL detr(ngrid,klev)
503      REAL entr(ngrid,klev)
504
505      REAL csc(ngrid,klev)
506
507      REAL zw2(ngrid,klev+1)
508      REAL w_est(ngrid,klev+1)
509      REAL f_star(ngrid,klev+1)
510      REAL wa_moy(ngrid,klev+1)
511
512      REAL ztva_est(ngrid,klev)
513      REAL zqla_est(ngrid,klev)
514      REAL zqsatth(ngrid,klev)
515      REAL zta_est(ngrid,klev)
516      REAL ztemp(ngrid),zqsat(ngrid)
517      REAL zdw2
518      REAL zw2modif
519      REAL zw2fact
520      REAL zeps(ngrid,klev)
521
522      REAL linter(ngrid)
523      INTEGER lmix(ngrid)
524      INTEGER lmix_bis(ngrid)
525      REAL    wmaxa(ngrid)
526
527      INTEGER ig,l,k
528
529      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
530      real zbuoybis
531      real zcor,zdelta,zcvm5,qlbef,zdz2
532      real betalpha,zbetalpha
533      real eps, afact
534      REAL REPS,RLvCp,DDT0
535      PARAMETER (DDT0=.01)
536      logical Zsat
537      LOGICAL active(ngrid),activetmp(ngrid)
538      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
539      REAL c2(ngrid,klev)
540
541      Zsat=.false.
542! Initialisation
543
544      RLvCp = RLVTT/RCPD
545      fact_epsilon=0.002
546      betalpha=0.9
547      afact=2./3.           
548
549      zbetalpha=betalpha/(1.+betalpha)
550
551
552! Initialisations des variables reeles
553if (1==0) then
554      ztva(:,:)=ztv(:,:)
555      ztva_est(:,:)=ztva(:,:)
556      ztla(:,:)=zthl(:,:)
557      zqta(:,:)=po(:,:)
558      zha(:,:) = ztva(:,:)
559else
560      ztva(:,:)=0.
561      ztva_est(:,:)=0.
562      ztla(:,:)=0.
563      zqta(:,:)=0.
564      zha(:,:) =0.
565endif
566
567      zqla_est(:,:)=0.
568      zqsatth(:,:)=0.
569      zqla(:,:)=0.
570      detr_star(:,:)=0.
571      entr_star(:,:)=0.
572      alim_star(:,:)=0.
573      alim_star_tot(:)=0.
574      csc(:,:)=0.
575      detr(:,:)=0.
576      entr(:,:)=0.
577      zw2(:,:)=0.
578      zbuoy(:,:)=0.
579      gamma(:,:)=0.
580      zeps(:,:)=0.
581      w_est(:,:)=0.
582      f_star(:,:)=0.
583      wa_moy(:,:)=0.
584      linter(:)=1.
585!     linter(:)=1.
586! Initialisation des variables entieres
587      lmix(:)=1
588      lmix_bis(:)=2
589      wmaxa(:)=0.
590      lalim(:)=1
591
592
593!-------------------------------------------------------------------------
594! On ne considere comme actif que les colonnes dont les deux premieres
595! couches sont instables.
596!-------------------------------------------------------------------------
597      active(:)=ztv(:,1)>ztv(:,2)
598
599!-------------------------------------------------------------------------
600! Definition de l'alimentation a l'origine dans thermcell_init
601!-------------------------------------------------------------------------
602      do l=1,klev-1
603         do ig=1,ngrid
604            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
605! L. Fita, LMD July 2014. Adding IF for zlev == 0. NOT NEEDED!!
606!!               IF ( zlev(ig,l+1) < 0.) THEN
607!!                 alim_star(ig,l) = 0.
608!!               ELSE
609                 alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
610     &                       *sqrt(zlev(ig,l+1))
611                 lalim(ig)=l+1
612!!               END IF
613               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
614            endif
615         enddo
616      enddo
617      do l=1,klev
618         do ig=1,ngrid
619            if (alim_star_tot(ig) > 1.e-10 ) then
620               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
621            endif
622         enddo
623      enddo
624      alim_star_tot(:)=1.
625!------------------------------------------------------------------------------
626! Calcul dans la premiere couche
627! On decide dans cette version que le thermique n'est actif que si la premiere
628! couche est instable.
629! Pourrait etre change si on veut que le thermiques puisse se déclencher
630! dans une couche l>1
631!------------------------------------------------------------------------------
632do ig=1,ngrid
633! Le panache va prendre au debut les caracteristiques de l'air contenu
634! dans cette couche.
635    if (active(ig)) then
636    ztla(ig,1)=zthl(ig,1)
637    zqta(ig,1)=po(ig,1)
638    zqla(ig,1)=zl(ig,1)
639!cr: attention, prise en compte de f*(1)=1
640    f_star(ig,2)=alim_star(ig,1)
641    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
642&                     *(zlev(ig,2)-zlev(ig,1))  &
643&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
644    w_est(ig,2)=zw2(ig,2)
645    endif
646enddo
647!
648
649!==============================================================================
650!boucle de calcul de la vitesse verticale dans le thermique
651!==============================================================================
652do l=2,klev-1
653!==============================================================================
654
655
656! On decide si le thermique est encore actif ou non
657! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
658    do ig=1,ngrid
659       active(ig)=active(ig) &
660&                 .and. zw2(ig,l)>1.e-10 &
661&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
662    enddo
663
664
665
666!---------------------------------------------------------------------------
667! calcul des proprietes thermodynamiques et de la vitesse de la couche l
668! sans tenir compte du detrainement et de l'entrainement dans cette
669! couche
670! C'est a dire qu'on suppose
671! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
672! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
673! avant) a l'alimentation pour avoir un calcul plus propre
674!---------------------------------------------------------------------------
675
676   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
677   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
678
679    do ig=1,ngrid
680!       print*,'active',active(ig),ig,l
681        if(active(ig)) then
682        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
683        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
684        zta_est(ig,l)=ztva_est(ig,l)
685        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
686        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
687     &      -zqla_est(ig,l))-zqla_est(ig,l))
688
689!------------------------------------------------
690!AJAM:nouveau calcul de w\B2 
691!------------------------------------------------
692              zdz=zlev(ig,l+1)-zlev(ig,l)
693              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
694
695              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
696              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
697              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
698 
699
700             if (w_est(ig,l+1).lt.0.) then
701                w_est(ig,l+1)=zw2(ig,l)
702             endif
703       endif
704    enddo
705
706
707!-------------------------------------------------
708!calcul des taux d'entrainement et de detrainement
709!-------------------------------------------------
710
711     do ig=1,ngrid
712        if (active(ig)) then
713
714          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
715          zw2m=w_est(ig,l+1)
716          zdz=zlev(ig,l+1)-zlev(ig,l)
717          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
718!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
719          zbuoybis=zbuoy(ig,l)
720          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
721          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
722
723         
724          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
725    &     afact*zbuoybis/zw2m - fact_epsilon )
726
727
728          detr_star(ig,l)=f_star(ig,l)*zdz                        &
729    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
730    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
731         
732! En dessous de lalim, on prend le max de alim_star et entr_star pour
733! alim_star et 0 sinon
734        if (l.lt.lalim(ig)) then
735          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
736          entr_star(ig,l)=0.
737        endif
738
739! Calcul du flux montant normalise
740      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
741     &              -detr_star(ig,l)
742
743      endif
744   enddo
745
746!----------------------------------------------------------------------------
747!calcul de la vitesse verticale en melangeant Tl et qt du thermique
748!---------------------------------------------------------------------------
749   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
750   do ig=1,ngrid
751       if (activetmp(ig)) then
752           Zsat=.false.
753           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
754     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
755     &            /(f_star(ig,l+1)+detr_star(ig,l))
756           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
757     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
758     &            /(f_star(ig,l+1)+detr_star(ig,l))
759
760        endif
761    enddo
762
763   ztemp(:)=zpspsk(:,l)*ztla(:,l)
764   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
765
766   do ig=1,ngrid
767      if (activetmp(ig)) then
768! on ecrit de maniere conservative (sat ou non)
769!          T = Tl +Lv/Cp ql
770           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
771           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
772           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
773!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
774           zha(ig,l) = ztva(ig,l)
775           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
776     &              -zqla(ig,l))-zqla(ig,l))
777           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
778           zdz=zlev(ig,l+1)-zlev(ig,l)
779           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
780
781            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
782            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
783            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
784      endif
785   enddo
786
787        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
788!
789!---------------------------------------------------------------------------
790!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
791!---------------------------------------------------------------------------
792
793   nbpb=0
794   do ig=1,ngrid
795            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
796!               stop'On tombe sur le cas particulier de thermcell_dry'
797!               print*,'On tombe sur le cas particulier de thermcell_plume'
798                nbpb=nbpb+1
799                zw2(ig,l+1)=0.
800                linter(ig)=l+1
801            endif
802
803        if (zw2(ig,l+1).lt.0.) then
804           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
805     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
806           zw2(ig,l+1)=0.
807        endif
808
809           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
810
811        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
812!   lmix est le niveau de la couche ou w (wa_moy) est maximum
813!on rajoute le calcul de lmix_bis
814            if (zqla(ig,l).lt.1.e-10) then
815               lmix_bis(ig)=l+1
816            endif
817            lmix(ig)=l+1
818            wmaxa(ig)=wa_moy(ig,l+1)
819        endif
820   enddo
821
822   if (nbpb>0) then
823   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
824   endif
825
826!=========================================================================
827! FIN DE LA BOUCLE VERTICALE
828      enddo
829!=========================================================================
830
831!on recalcule alim_star_tot
832       do ig=1,ngrid
833          alim_star_tot(ig)=0.
834       enddo
835       do ig=1,ngrid
836          do l=1,lalim(ig)-1
837          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
838          enddo
839       enddo
840       
841
842        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
843
844#undef wrgrads_thermcell
845#ifdef wrgrads_thermcell
846         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
847         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
848         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
849         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
850         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
851         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
852         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
853#endif
854
855
856     return
857     end
858
Note: See TracBrowser for help on using the repository browser.