source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90 @ 5225

Last change on this file since 5225 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

  • calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

  • fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

  • physiq.F : appel modifie a calltherm et fisrtilp
  • thermcell_main.F90 : changement dans les flag_thermals_ed
  • thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

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