source: trunk/LMDZ.EARTH/libf/phylmd/thermcell_plume.F90 @ 357

Last change on this file since 357 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 27.3 KB
Line 
1!
2! $Id: thermcell_plume.F90 1403 2010-07-01 09:02:53Z fairhead $
3!
4      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
5     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
6     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
7     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8     &           ,lev_out,lunout1,igout)
9
10!--------------------------------------------------------------------------
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
103! Initialisations des variables reeles
104if (1==1) then
105      ztva(:,:)=ztv(:,:)
106      ztva_est(:,:)=ztva(:,:)
107      ztla(:,:)=zthl(:,:)
108      zqta(:,:)=po(:,:)
109      zha(:,:) = ztva(:,:)
110else
111      ztva(:,:)=0.
112      ztva_est(:,:)=0.
113      ztla(:,:)=0.
114      zqta(:,:)=0.
115      zha(:,:) =0.
116endif
117
118      zqla_est(:,:)=0.
119      zqsatth(:,:)=0.
120      zqla(:,:)=0.
121      detr_star(:,:)=0.
122      entr_star(:,:)=0.
123      alim_star(:,:)=0.
124      alim_star_tot(:)=0.
125      csc(:,:)=0.
126      detr(:,:)=0.
127      entr(:,:)=0.
128      zw2(:,:)=0.
129      w_est(:,:)=0.
130      f_star(:,:)=0.
131      wa_moy(:,:)=0.
132      linter(:)=1.
133      linter(:)=1.
134
135! Initialisation des variables entieres
136      lmix(:)=1
137      lmix_bis(:)=2
138      wmaxa(:)=0.
139      lalim(:)=1
140
141!-------------------------------------------------------------------------
142! On ne considere comme actif que les colonnes dont les deux premieres
143! couches sont instables.
144!-------------------------------------------------------------------------
145      active(:)=ztv(:,1)>ztv(:,2)
146
147!-------------------------------------------------------------------------
148! Definition de l'alimentation a l'origine dans thermcell_init
149!-------------------------------------------------------------------------
150      do l=1,klev-1
151         do ig=1,ngrid
152            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
153               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
154     &                       *sqrt(zlev(ig,l+1))
155               lalim(:)=l+1
156               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
157            endif
158         enddo
159      enddo
160      do l=1,klev
161         do ig=1,ngrid
162            if (alim_star_tot(ig) > 1.e-10 ) then
163               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
164            endif
165         enddo
166      enddo
167      alim_star_tot(:)=1.
168
169
170!------------------------------------------------------------------------------
171! Calcul dans la premiere couche
172! On decide dans cette version que le thermique n'est actif que si la premiere
173! couche est instable.
174! Pourrait etre change si on veut que le thermiques puisse se déclencher
175! dans une couche l>1
176!------------------------------------------------------------------------------
177do ig=1,ngrid
178! Le panache va prendre au debut les caracteristiques de l'air contenu
179! dans cette couche.
180    if (active(ig)) then
181    ztla(ig,1)=zthl(ig,1)
182    zqta(ig,1)=po(ig,1)
183    zqla(ig,1)=zl(ig,1)
184!cr: attention, prise en compte de f*(1)=1
185    f_star(ig,2)=alim_star(ig,1)
186    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
187&                     *(zlev(ig,2)-zlev(ig,1))  &
188&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
189    w_est(ig,2)=zw2(ig,2)
190    endif
191enddo
192!
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        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
346! on ecrit de maniere conservative (sat ou non)
347!          T = Tl +Lv/Cp ql
348           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
349           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
350!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
351           zha(ig,l) = ztva(ig,l)
352           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
353     &              -zqla(ig,l))-zqla(ig,l))
354
355!on ecrit zqsat
356           zqsatth(ig,l)=qsatbef 
357
358!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
359!          zw2(ig,l+1)=&
360!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
361!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
362!     &                 *1./(1.+fact_gamma) &
363!     &                 *(zlev(ig,l+1)-zlev(ig,l))
364!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
365! La meme en plus modulaire :
366           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
367           zdz=zlev(ig,l+1)-zlev(ig,l)
368
369
370           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
371
372if (1==0) then
373           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
374           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
375           zw2(ig,l+1)=zw2modif+zdw2
376else
377           zdrag=fact_epsilon/(zalpha**expa)
378           zw2fact=zbuoy/zdrag*a1
379           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
380      &    +zw2fact
381
382
383endif
384
385      endif
386   enddo
387
388        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
389!
390!---------------------------------------------------------------------------
391!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
392!---------------------------------------------------------------------------
393
394   do ig=1,ngrid
395            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
396!               stop'On tombe sur le cas particulier de thermcell_dry'
397                print*,'On tombe sur le cas particulier de thermcell_plume'
398                zw2(ig,l+1)=0.
399                linter(ig)=l+1
400            endif
401
402        if (zw2(ig,l+1).lt.0.) then
403           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
404     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
405           zw2(ig,l+1)=0.
406        endif
407
408           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
409
410        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
411!   lmix est le niveau de la couche ou w (wa_moy) est maximum
412!on rajoute le calcul de lmix_bis
413            if (zqla(ig,l).lt.1.e-10) then
414               lmix_bis(ig)=l+1
415            endif
416            lmix(ig)=l+1
417            wmaxa(ig)=wa_moy(ig,l+1)
418        endif
419   enddo
420
421!=========================================================================
422! FIN DE LA BOUCLE VERTICALE
423      enddo
424!=========================================================================
425
426!on recalcule alim_star_tot
427       do ig=1,ngrid
428          alim_star_tot(ig)=0.
429       enddo
430       do ig=1,ngrid
431          do l=1,lalim(ig)-1
432          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
433          enddo
434       enddo
435       
436
437        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
438
439
440      return
441      end
442
443!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
444!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
445!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449 SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
450&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
451&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
452&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
453&           ,lev_out,lunout1,igout)
454
455!--------------------------------------------------------------------------
456!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
457! Version conforme a l'article de Rio et al. 2010.
458! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
459!--------------------------------------------------------------------------
460
461      IMPLICIT NONE
462
463#include "YOMCST.h"
464#include "YOETHF.h"
465#include "FCTTRE.h"
466#include "iniprint.h"
467#include "thermcell.h"
468
469      INTEGER itap
470      INTEGER lunout1,igout
471      INTEGER ngrid,klev
472      REAL ptimestep
473      REAL ztv(ngrid,klev)
474      REAL zthl(ngrid,klev)
475      REAL po(ngrid,klev)
476      REAL zl(ngrid,klev)
477      REAL rhobarz(ngrid,klev)
478      REAL zlev(ngrid,klev+1)
479      REAL pplev(ngrid,klev+1)
480      REAL pphi(ngrid,klev)
481      REAL zpspsk(ngrid,klev)
482      REAL alim_star(ngrid,klev)
483      REAL f0(ngrid)
484      INTEGER lalim(ngrid)
485      integer lev_out                           ! niveau pour les print
486   
487      real alim_star_tot(ngrid)
488
489      REAL ztva(ngrid,klev)
490      REAL ztla(ngrid,klev)
491      REAL zqla(ngrid,klev)
492      REAL zqta(ngrid,klev)
493      REAL zha(ngrid,klev)
494
495      REAL detr_star(ngrid,klev)
496      REAL coefc
497      REAL entr_star(ngrid,klev)
498      REAL detr(ngrid,klev)
499      REAL entr(ngrid,klev)
500
501      REAL csc(ngrid,klev)
502
503      REAL zw2(ngrid,klev+1)
504      REAL w_est(ngrid,klev+1)
505      REAL f_star(ngrid,klev+1)
506      REAL wa_moy(ngrid,klev+1)
507
508      REAL ztva_est(ngrid,klev)
509      REAL zqla_est(ngrid,klev)
510      REAL zqsatth(ngrid,klev)
511      REAL zta_est(ngrid,klev)
512      REAL ztemp(ngrid),zqsat(ngrid)
513      REAL zdw2
514      REAL zw2modif
515      REAL zw2fact
516      REAL zeps(ngrid,klev)
517
518      REAL linter(ngrid)
519      INTEGER lmix(ngrid)
520      INTEGER lmix_bis(ngrid)
521      REAL    wmaxa(ngrid)
522
523      INTEGER ig,l,k
524
525      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
526      real zbuoybis
527      real zcor,zdelta,zcvm5,qlbef,zdz2
528      real betalpha,zbetalpha
529      real eps, afact
530      REAL REPS,RLvCp,DDT0
531      PARAMETER (DDT0=.01)
532      logical Zsat
533      LOGICAL active(ngrid),activetmp(ngrid)
534      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
535      REAL c2(ngrid,klev)
536      Zsat=.false.
537! Initialisation
538
539      RLvCp = RLVTT/RCPD
540      fact_epsilon=0.002
541      betalpha=0.9
542      afact=2./3.           
543
544      zbetalpha=betalpha/(1.+betalpha)
545
546
547! Initialisations des variables reeles
548if (1==0) then
549      ztva(:,:)=ztv(:,:)
550      ztva_est(:,:)=ztva(:,:)
551      ztla(:,:)=zthl(:,:)
552      zqta(:,:)=po(:,:)
553      zha(:,:) = ztva(:,:)
554else
555      ztva(:,:)=0.
556      ztva_est(:,:)=0.
557      ztla(:,:)=0.
558      zqta(:,:)=0.
559      zha(:,:) =0.
560endif
561
562      zqla_est(:,:)=0.
563      zqsatth(:,:)=0.
564      zqla(:,:)=0.
565      detr_star(:,:)=0.
566      entr_star(:,:)=0.
567      alim_star(:,:)=0.
568      alim_star_tot(:)=0.
569      csc(:,:)=0.
570      detr(:,:)=0.
571      entr(:,:)=0.
572      zw2(:,:)=0.
573      zbuoy(:,:)=0.
574      gamma(:,:)=0.
575      zeps(:,:)=0.
576      w_est(:,:)=0.
577      f_star(:,:)=0.
578      wa_moy(:,:)=0.
579      linter(:)=1.
580!     linter(:)=1.
581! Initialisation des variables entieres
582      lmix(:)=1
583      lmix_bis(:)=2
584      wmaxa(:)=0.
585      lalim(:)=1
586
587
588!-------------------------------------------------------------------------
589! On ne considere comme actif que les colonnes dont les deux premieres
590! couches sont instables.
591!-------------------------------------------------------------------------
592      active(:)=ztv(:,1)>ztv(:,2)
593
594!-------------------------------------------------------------------------
595! Definition de l'alimentation a l'origine dans thermcell_init
596!-------------------------------------------------------------------------
597      do l=1,klev-1
598         do ig=1,ngrid
599            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
600               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
601     &                       *sqrt(zlev(ig,l+1))
602               lalim(ig)=l+1
603               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
604            endif
605         enddo
606      enddo
607      do l=1,klev
608         do ig=1,ngrid
609            if (alim_star_tot(ig) > 1.e-10 ) then
610               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
611            endif
612         enddo
613      enddo
614      alim_star_tot(:)=1.
615
616
617
618!------------------------------------------------------------------------------
619! Calcul dans la premiere couche
620! On decide dans cette version que le thermique n'est actif que si la premiere
621! couche est instable.
622! Pourrait etre change si on veut que le thermiques puisse se déclencher
623! dans une couche l>1
624!------------------------------------------------------------------------------
625do ig=1,ngrid
626! Le panache va prendre au debut les caracteristiques de l'air contenu
627! dans cette couche.
628    if (active(ig)) then
629    ztla(ig,1)=zthl(ig,1)
630    zqta(ig,1)=po(ig,1)
631    zqla(ig,1)=zl(ig,1)
632!cr: attention, prise en compte de f*(1)=1
633    f_star(ig,2)=alim_star(ig,1)
634    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
635&                     *(zlev(ig,2)-zlev(ig,1))  &
636&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
637    w_est(ig,2)=zw2(ig,2)
638    endif
639enddo
640!
641
642!==============================================================================
643!boucle de calcul de la vitesse verticale dans le thermique
644!==============================================================================
645do l=2,klev-1
646!==============================================================================
647
648
649! On decide si le thermique est encore actif ou non
650! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
651    do ig=1,ngrid
652       active(ig)=active(ig) &
653&                 .and. zw2(ig,l)>1.e-10 &
654&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
655    enddo
656
657
658
659!---------------------------------------------------------------------------
660! calcul des proprietes thermodynamiques et de la vitesse de la couche l
661! sans tenir compte du detrainement et de l'entrainement dans cette
662! couche
663! C'est a dire qu'on suppose
664! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
665! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
666! avant) a l'alimentation pour avoir un calcul plus propre
667!---------------------------------------------------------------------------
668
669   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
670   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
671
672    do ig=1,ngrid
673!       print*,'active',active(ig),ig,l
674        if(active(ig)) then
675        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
676        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
677        zta_est(ig,l)=ztva_est(ig,l)
678        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
679        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
680     &      -zqla_est(ig,l))-zqla_est(ig,l))
681
682!------------------------------------------------
683!AJAM:nouveau calcul de w² 
684!------------------------------------------------
685              zdz=zlev(ig,l+1)-zlev(ig,l)
686              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
687
688              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
689              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
690              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
691 
692
693             if (w_est(ig,l+1).lt.0.) then
694                w_est(ig,l+1)=zw2(ig,l)
695             endif
696       endif
697    enddo
698
699
700!-------------------------------------------------
701!calcul des taux d'entrainement et de detrainement
702!-------------------------------------------------
703
704     do ig=1,ngrid
705        if (active(ig)) then
706
707          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
708          zw2m=w_est(ig,l+1)
709          zdz=zlev(ig,l+1)-zlev(ig,l)
710          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
711!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
712          zbuoybis=zbuoy(ig,l)
713          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
714          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
715
716         
717          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
718    &     afact*zbuoybis/zw2m - fact_epsilon )
719
720
721          detr_star(ig,l)=f_star(ig,l)*zdz                        &
722    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
723    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
724         
725! En dessous de lalim, on prend le max de alim_star et entr_star pour
726! alim_star et 0 sinon
727        if (l.lt.lalim(ig)) then
728          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
729          entr_star(ig,l)=0.
730        endif
731
732! Calcul du flux montant normalise
733      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
734     &              -detr_star(ig,l)
735
736      endif
737   enddo
738
739
740!----------------------------------------------------------------------------
741!calcul de la vitesse verticale en melangeant Tl et qt du thermique
742!---------------------------------------------------------------------------
743   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
744   do ig=1,ngrid
745       if (activetmp(ig)) then
746           Zsat=.false.
747           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
748     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
749     &            /(f_star(ig,l+1)+detr_star(ig,l))
750           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
751     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
752     &            /(f_star(ig,l+1)+detr_star(ig,l))
753
754        endif
755    enddo
756
757   ztemp(:)=zpspsk(:,l)*ztla(:,l)
758   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
759
760   do ig=1,ngrid
761      if (activetmp(ig)) then
762        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
763! on ecrit de maniere conservative (sat ou non)
764!          T = Tl +Lv/Cp ql
765           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
766           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
767           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
768!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
769           zha(ig,l) = ztva(ig,l)
770           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
771     &              -zqla(ig,l))-zqla(ig,l))
772           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
773           zdz=zlev(ig,l+1)-zlev(ig,l)
774           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
775
776            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
777            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
778            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
779      endif
780   enddo
781
782        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
783!
784!---------------------------------------------------------------------------
785!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
786!---------------------------------------------------------------------------
787
788   do ig=1,ngrid
789            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
790!               stop'On tombe sur le cas particulier de thermcell_dry'
791                print*,'On tombe sur le cas particulier de thermcell_plume'
792                zw2(ig,l+1)=0.
793                linter(ig)=l+1
794            endif
795
796        if (zw2(ig,l+1).lt.0.) then
797           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
798     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
799           zw2(ig,l+1)=0.
800        endif
801
802           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
803
804        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
805!   lmix est le niveau de la couche ou w (wa_moy) est maximum
806!on rajoute le calcul de lmix_bis
807            if (zqla(ig,l).lt.1.e-10) then
808               lmix_bis(ig)=l+1
809            endif
810            lmix(ig)=l+1
811            wmaxa(ig)=wa_moy(ig,l+1)
812        endif
813   enddo
814
815!=========================================================================
816! FIN DE LA BOUCLE VERTICALE
817      enddo
818!=========================================================================
819
820!on recalcule alim_star_tot
821       do ig=1,ngrid
822          alim_star_tot(ig)=0.
823       enddo
824       do ig=1,ngrid
825          do l=1,lalim(ig)-1
826          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
827          enddo
828       enddo
829       
830
831        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
832
833#undef wrgrads_thermcell
834#ifdef wrgrads_thermcell
835         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
836         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
837         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
838         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
839         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
840         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
841         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
842#endif
843
844
845     return
846     end
847
Note: See TracBrowser for help on using the repository browser.