source: LMDZ5/branches/LMDZ5_AR5/libf/phylmd/thermcell_plume.F90 @ 5434

Last change on this file since 5434 was 1503, checked in by idelkadi, 14 years ago

Correction de bug et vectorisation de thermcell_plume.F90

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.4 KB
Line 
1!
2! $Id: thermcell_plume.F90 1503 2011-03-23 11:57:52Z fhourdin $
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! Initialisation
93      RLvCp = RLVTT/RCPD
94     
95     
96         fact_epsilon=0.002
97         a1=2./3.
98         fact_gamma=0.9
99         zfact=fact_gamma/(1+fact_gamma)
100         fact_gamma2=zfact
101         expa=0.
102     
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(ig)=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! 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
372!if (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
376!else
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!endif
382
383      endif
384   enddo
385
386        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
387!
388!---------------------------------------------------------------------------
389!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
390!---------------------------------------------------------------------------
391
392   nbpb=0
393   do ig=1,ngrid
394            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
395!               stop'On tombe sur le cas particulier de thermcell_dry'
396!               print*,'On tombe sur le cas particulier de thermcell_plume'
397                nbpb=nbpb+1
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   if (nbpb>0) then
422   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
423   endif
424
425!=========================================================================
426! FIN DE LA BOUCLE VERTICALE
427      enddo
428!=========================================================================
429
430!on recalcule alim_star_tot
431       do ig=1,ngrid
432          alim_star_tot(ig)=0.
433       enddo
434       do ig=1,ngrid
435          do l=1,lalim(ig)-1
436          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
437          enddo
438       enddo
439       
440
441        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
442
443
444      return
445      end
446
447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
453 SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
454&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
455&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
456&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
457&           ,lev_out,lunout1,igout)
458
459!--------------------------------------------------------------------------
460!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
461! Version conforme a l'article de Rio et al. 2010.
462! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
463!--------------------------------------------------------------------------
464
465      IMPLICIT NONE
466
467#include "YOMCST.h"
468#include "YOETHF.h"
469#include "FCTTRE.h"
470#include "iniprint.h"
471#include "thermcell.h"
472
473      INTEGER itap
474      INTEGER lunout1,igout
475      INTEGER ngrid,klev
476      REAL ptimestep
477      REAL ztv(ngrid,klev)
478      REAL zthl(ngrid,klev)
479      REAL po(ngrid,klev)
480      REAL zl(ngrid,klev)
481      REAL rhobarz(ngrid,klev)
482      REAL zlev(ngrid,klev+1)
483      REAL pplev(ngrid,klev+1)
484      REAL pphi(ngrid,klev)
485      REAL zpspsk(ngrid,klev)
486      REAL alim_star(ngrid,klev)
487      REAL f0(ngrid)
488      INTEGER lalim(ngrid)
489      integer lev_out                           ! niveau pour les print
490      integer nbpb
491   
492      real alim_star_tot(ngrid)
493
494      REAL ztva(ngrid,klev)
495      REAL ztla(ngrid,klev)
496      REAL zqla(ngrid,klev)
497      REAL zqta(ngrid,klev)
498      REAL zha(ngrid,klev)
499
500      REAL detr_star(ngrid,klev)
501      REAL coefc
502      REAL entr_star(ngrid,klev)
503      REAL detr(ngrid,klev)
504      REAL entr(ngrid,klev)
505
506      REAL csc(ngrid,klev)
507
508      REAL zw2(ngrid,klev+1)
509      REAL w_est(ngrid,klev+1)
510      REAL f_star(ngrid,klev+1)
511      REAL wa_moy(ngrid,klev+1)
512
513      REAL ztva_est(ngrid,klev)
514      REAL zqla_est(ngrid,klev)
515      REAL zqsatth(ngrid,klev)
516      REAL zta_est(ngrid,klev)
517      REAL ztemp(ngrid),zqsat(ngrid)
518      REAL zdw2
519      REAL zw2modif
520      REAL zw2fact
521      REAL zeps(ngrid,klev)
522
523      REAL linter(ngrid)
524      INTEGER lmix(ngrid)
525      INTEGER lmix_bis(ngrid)
526      REAL    wmaxa(ngrid)
527
528      INTEGER ig,l,k
529
530      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
531      real zbuoybis
532      real zcor,zdelta,zcvm5,qlbef,zdz2
533      real betalpha,zbetalpha
534      real eps, afact
535      REAL REPS,RLvCp,DDT0
536      PARAMETER (DDT0=.01)
537      logical Zsat
538      LOGICAL active(ngrid),activetmp(ngrid)
539      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
540      REAL c2(ngrid,klev)
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               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     do ig=1,ngrid
710        if (active(ig)) then
711
712          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
713          zw2m=w_est(ig,l+1)
714          zdz=zlev(ig,l+1)-zlev(ig,l)
715          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
716!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
717          zbuoybis=zbuoy(ig,l)
718          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
719          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
720
721         
722          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
723    &     afact*zbuoybis/zw2m - fact_epsilon )
724
725
726          detr_star(ig,l)=f_star(ig,l)*zdz                        &
727    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
728    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
729         
730! En dessous de lalim, on prend le max de alim_star et entr_star pour
731! alim_star et 0 sinon
732        if (l.lt.lalim(ig)) then
733          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
734          entr_star(ig,l)=0.
735        endif
736
737! Calcul du flux montant normalise
738      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
739     &              -detr_star(ig,l)
740
741      endif
742   enddo
743
744
745!----------------------------------------------------------------------------
746!calcul de la vitesse verticale en melangeant Tl et qt du thermique
747!---------------------------------------------------------------------------
748   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
749   do ig=1,ngrid
750       if (activetmp(ig)) then
751           Zsat=.false.
752           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
753     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
754     &            /(f_star(ig,l+1)+detr_star(ig,l))
755           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
756     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
757     &            /(f_star(ig,l+1)+detr_star(ig,l))
758
759        endif
760    enddo
761
762   ztemp(:)=zpspsk(:,l)*ztla(:,l)
763   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
764
765   do ig=1,ngrid
766      if (activetmp(ig)) then
767! on ecrit de maniere conservative (sat ou non)
768!          T = Tl +Lv/Cp ql
769           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
770           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
771           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
772!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
773           zha(ig,l) = ztva(ig,l)
774           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
775     &              -zqla(ig,l))-zqla(ig,l))
776           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
777           zdz=zlev(ig,l+1)-zlev(ig,l)
778           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
779
780            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
781            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
782            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
783      endif
784   enddo
785
786        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
787!
788!---------------------------------------------------------------------------
789!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
790!---------------------------------------------------------------------------
791
792   nbpb=0
793   do ig=1,ngrid
794            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
795!               stop'On tombe sur le cas particulier de thermcell_dry'
796!               print*,'On tombe sur le cas particulier de thermcell_plume'
797                nbpb=nbpb+1
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   if (nbpb>0) then
822   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
823   endif
824
825!=========================================================================
826! FIN DE LA BOUCLE VERTICALE
827      enddo
828!=========================================================================
829
830!on recalcule alim_star_tot
831       do ig=1,ngrid
832          alim_star_tot(ig)=0.
833       enddo
834       do ig=1,ngrid
835          do l=1,lalim(ig)-1
836          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
837          enddo
838       enddo
839       
840
841        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
842
843#undef wrgrads_thermcell
844#ifdef wrgrads_thermcell
845         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
846         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
847         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
848         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
849         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
850         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
851         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
852#endif
853
854
855     return
856     end
857
Note: See TracBrowser for help on using the repository browser.