source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/thermcell_plume.F90 @ 3655

Last change on this file since 3655 was 1057, checked in by lmdzadmin, 16 years ago

Un peu de nettoyage
ACA/FH/IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.8 KB
Line 
1      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
2     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
3     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
4     &           ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
5     &           ,lev_out,lunout1,igout)
6
7!--------------------------------------------------------------------------
8!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
9!--------------------------------------------------------------------------
10
11      IMPLICIT NONE
12
13#include "YOMCST.h"
14#include "YOETHF.h"
15#include "FCTTRE.h"
16#include "iniprint.h"
17#include "thermcell.h"
18
19      INTEGER itap
20      INTEGER lunout1,igout
21      INTEGER ngrid,klev
22      REAL ptimestep
23      REAL ztv(ngrid,klev)
24      REAL zthl(ngrid,klev)
25      REAL po(ngrid,klev)
26      REAL zl(ngrid,klev)
27      REAL rhobarz(ngrid,klev)
28      REAL zlev(ngrid,klev+1)
29      REAL pplev(ngrid,klev+1)
30      REAL pphi(ngrid,klev)
31      REAL zpspsk(ngrid,klev)
32      REAL alim_star(ngrid,klev)
33      REAL zmax_sec(ngrid)
34      REAL f0(ngrid)
35      REAL l_mix
36      REAL r_aspect
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 zqla0(ngrid,klev)
47      REAL zqta(ngrid,klev)
48      REAL zha(ngrid,klev)
49
50      REAL detr_star(ngrid,klev)
51      REAL coefc
52      REAL detr_stara(ngrid,klev)
53      REAL detr_starb(ngrid,klev)
54      REAL detr_starc(ngrid,klev)
55      REAL detr_star0(ngrid,klev)
56      REAL detr_star1(ngrid,klev)
57      REAL detr_star2(ngrid,klev)
58
59      REAL entr_star(ngrid,klev)
60      REAL entr_star1(ngrid,klev)
61      REAL entr_star2(ngrid,klev)
62      REAL detr(ngrid,klev)
63      REAL entr(ngrid,klev)
64
65      REAL zw2(ngrid,klev+1)
66      REAL w_est(ngrid,klev+1)
67      REAL f_star(ngrid,klev+1)
68      REAL wa_moy(ngrid,klev+1)
69
70      REAL ztva_est(ngrid,klev)
71      REAL zqla_est(ngrid,klev)
72      REAL zqsatth(ngrid,klev)
73      REAL zta_est(ngrid,klev)
74
75      REAL linter(ngrid)
76      INTEGER lmix(ngrid)
77      INTEGER lmix_bis(ngrid)
78      REAL    wmaxa(ngrid)
79
80      INTEGER ig,l,k
81
82      real zcor,zdelta,zcvm5,qlbef
83      real Tbef,qsatbef
84      real dqsat_dT,DT,num,denom
85      REAL REPS,RLvCp,DDT0
86      PARAMETER (DDT0=.01)
87      logical Zsat
88      REAL fact_gamma,fact_epsilon
89      REAL c2(ngrid,klev)
90
91      Zsat=.false.
92! Initialisation
93      RLvCp = RLVTT/RCPD
94     
95      if (iflag_thermals_ed==0) then
96         fact_gamma=1.
97         fact_epsilon=1.
98      else if (iflag_thermals_ed==1)  then
99         fact_gamma=1.
100         fact_epsilon=1.
101      else if (iflag_thermals_ed==2)  then
102         fact_gamma=1.
103         fact_epsilon=2.
104      endif
105
106      do l=1,klev
107         do ig=1,ngrid
108            zqla_est(ig,l)=0.
109            ztva_est(ig,l)=ztva(ig,l)
110            zqsatth(ig,l)=0.
111         enddo
112      enddo
113
114!CR: attention test couche alim
115!     do l=2,klev
116!     do ig=1,ngrid
117!        alim_star(ig,l)=0.
118!     enddo
119!     enddo
120!AM:initialisations du thermique
121      do k=1,klev
122         do ig=1,ngrid
123            ztva(ig,k)=ztv(ig,k)
124            ztla(ig,k)=zthl(ig,k)
125            zqla(ig,k)=0.
126            zqta(ig,k)=po(ig,k)
127!
128            ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
129            ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
130            zha(ig,k) = ztva(ig,k)
131!
132         enddo
133      enddo
134      do k=1,klev
135        do ig=1,ngrid
136           detr_star(ig,k)=0.
137           entr_star(ig,k)=0.
138
139           detr_stara(ig,k)=0.
140           detr_starb(ig,k)=0.
141           detr_starc(ig,k)=0.
142           detr_star0(ig,k)=0.
143           zqla0(ig,k)=0.
144           detr_star1(ig,k)=0.
145           detr_star2(ig,k)=0.
146           entr_star1(ig,k)=0.
147           entr_star2(ig,k)=0.
148
149           detr(ig,k)=0.
150           entr(ig,k)=0.
151        enddo
152      enddo
153      if (prt_level.ge.1) print*,'7 OK convect8'
154      do k=1,klev+1
155         do ig=1,ngrid
156            zw2(ig,k)=0.
157            w_est(ig,k)=0.
158            f_star(ig,k)=0.
159            wa_moy(ig,k)=0.
160         enddo
161      enddo
162
163      if (prt_level.ge.1) print*,'8 OK convect8'
164      do ig=1,ngrid
165         linter(ig)=1.
166         lmix(ig)=1
167         lmix_bis(ig)=2
168         wmaxa(ig)=0.
169      enddo
170
171!-----------------------------------------------------------------------------------
172!boucle de calcul de la vitesse verticale dans le thermique
173!-----------------------------------------------------------------------------------
174      do l=1,klev-1
175         do ig=1,ngrid
176
177
178
179! Calcul dans la premiere couche active du thermique (ce qu'on teste
180! en disant que la couche est instable et que w2 en bas de la couche
181! est nulle.
182
183            if (ztv(ig,l).gt.ztv(ig,l+1)  &
184     &         .and.alim_star(ig,l).gt.1.e-10  &
185     &         .and.zw2(ig,l).lt.1e-10) then
186
187
188! Le panache va prendre au debut les caracteristiques de l'air contenu
189! dans cette couche.
190               ztla(ig,l)=zthl(ig,l)
191               zqta(ig,l)=po(ig,l)
192               zqla(ig,l)=zl(ig,l)
193               f_star(ig,l+1)=alim_star(ig,l)
194
195               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
196     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
197     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
198               w_est(ig,l+1)=zw2(ig,l+1)
199!
200
201
202            else if ((zw2(ig,l).ge.1e-10).and.  &
203     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
204!estimation du detrainement a partir de la geometrie du pas precedent
205!tests sur la definition du detr
206!calcul de detr_star et entr_star
207
208
209
210!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
211! FH le test miraculeux de Catherine ? Le bout du tunel ?
212!               w_est(ig,3)=zw2(ig,2)*  &
213!    &                   ((f_star(ig,2))**2)  &
214!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
215!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
216!    &                   *(zlev(ig,3)-zlev(ig,2))
217!     if (l.gt.2) then
218!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
219
220
221
222! Premier calcul de la vitesse verticale a partir de la temperature
223! potentielle virtuelle
224
225! FH CESTQUOI CA ????
226#define int1d2
227!#undef int1d2
228#ifdef int1d2
229      if (l.ge.2) then
230#else
231      if (l.gt.2) then
232#endif
233
234      if (1.eq.1) then
235          w_est(ig,3)=zw2(ig,2)* &
236     &      ((f_star(ig,2))**2) &
237     &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
238     &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
239!     &      *1./3. &
240     &      *(zlev(ig,3)-zlev(ig,2))
241       endif
242
243
244!---------------------------------------------------------------------------
245!calcul de l entrainement et du detrainement lateral
246!---------------------------------------------------------------------------
247!
248!test:estimation de ztva_new_est sans entrainement
249
250               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
251               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
252               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
253               qsatbef=MIN(0.5,qsatbef)
254               zcor=1./(1.-retv*qsatbef)
255               qsatbef=qsatbef*zcor
256               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
257               if (Zsat) then
258               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
259               DT = 0.5*RLvCp*qlbef
260               do while (abs(DT).gt.DDT0)
261                 Tbef=Tbef+DT
262                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
263                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
264                 qsatbef=MIN(0.5,qsatbef)
265                 zcor=1./(1.-retv*qsatbef)
266                 qsatbef=qsatbef*zcor
267                 qlbef=zqta(ig,l-1)-qsatbef
268
269                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
270                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
271                 zcor=1./(1.-retv*qsatbef)
272                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
273                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
274                 denom=1.+RLvCp*dqsat_dT
275                 DT=num/denom
276               enddo
277                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
278               endif
279        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
280        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
281        zta_est(ig,l)=ztva_est(ig,l)
282        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
283     &      -zqla_est(ig,l))-zqla_est(ig,l))
284
285             w_est(ig,l+1)=zw2(ig,l)*  &
286     &                   ((f_star(ig,l))**2)  &
287     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
288     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
289!     &                   *1./3. &
290     &                   *(zlev(ig,l+1)-zlev(ig,l))
291             if (w_est(ig,l+1).lt.0.) then
292                w_est(ig,l+1)=zw2(ig,l)
293             endif
294!
295!calcul du detrainement
296!=======================
297
298!CR:on vire les modifs
299         if (iflag_thermals_ed==0) then
300
301! Modifications du calcul du detrainement.
302! Dans la version de la these de Catherine, on passe brusquement
303! de la version seche a la version nuageuse pour le detrainement
304! ce qui peut occasioner des oscillations.
305! dans la nouvelle version, on commence par calculer un detrainement sec.
306! Puis un autre en cas de nuages.
307! Puis on combine les deux lineairement en fonction de la quantite d'eau.
308
309#define int1d3
310!#undef int1d3
311#define RIO_TH
312#ifdef RIO_TH
313!1. Cas non nuageux
314! 1.1 on est sous le zmax_sec et w croit
315          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
316     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
317#ifdef int1d3
318     &       (zqla_est(ig,l).lt.1.e-10)) then
319#else
320     &       (zqla(ig,l-1).lt.1.e-10)) then
321#endif
322             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
323     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
324     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
325     &       /(r_aspect*zmax_sec(ig)))
326             detr_stara(ig,l)=detr_star(ig,l)
327
328       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
329
330! 1.2 on est sous le zmax_sec et w decroit
331          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
332#ifdef int1d3
333     &            (zqla_est(ig,l).lt.1.e-10)) then
334#else
335     &            (zqla(ig,l-1).lt.1.e-10)) then
336#endif
337             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
338     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
339     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
340     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
341     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
342     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
343     &       *((zmax_sec(ig)-zlev(ig,l))/  &
344     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
345             detr_starb(ig,l)=detr_star(ig,l)
346
347        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
348
349          else
350
351! 1.3 dans les autres cas
352             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
353     &                      *(zlev(ig,l+1)-zlev(ig,l))
354             detr_starc(ig,l)=detr_star(ig,l)
355
356        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
357             
358          endif
359
360#else
361
362! 1.1 on est sous le zmax_sec et w croit
363          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
364     &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
365             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
366     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
367     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
368     &       /(r_aspect*zmax_sec(ig)))
369
370       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
371
372! 1.2 on est sous le zmax_sec et w decroit
373          else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
374             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
375     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
376     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
377     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
378     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
379     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
380     &       *((zmax_sec(ig)-zlev(ig,l))/  &
381     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
382       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
383
384          else
385             detr_star=0.
386          endif
387
388! 1.3 dans les autres cas
389          detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
390     &                      *(zlev(ig,l+1)-zlev(ig,l))
391
392          coefc=min(zqla(ig,l-1)/1.e-3,1.)
393          if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
394          coefc=1.
395! il semble qu'il soit important de baser le calcul sur
396! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
397          detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
398
399        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
400
401#endif
402
403
404        if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
405!IM 730508 beg
406!        if(itap.GE.7200) THEN
407!         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
408!        endif
409!IM 730508 end
410         
411         zqla0(ig,l)=zqla_est(ig,l)
412         detr_star0(ig,l)=detr_star(ig,l)
413!IM 060508 beg
414!         if(detr_star(ig,l).GT.1.) THEN
415!          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
416!   &      ,detr_starc(ig,l),coefc
417!         endif
418!IM 060508 end
419!IM 160508 beg
420!IM 160508       IF (f0(ig).NE.0.) THEN
421           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
422!IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
423!IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
424!IM 160508       ELSE
425!IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
426!IM 160508       ENDIF
427!IM 160508 end
428!IM 060508 beg
429!        if(detr_star(ig,l).GT.1.) THEN
430!         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
431!   &     float(1)/f0(ig)
432!        endif
433!IM 060508 end
434        if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
435!
436!calcul de entr_star
437
438! #undef test2
439! #ifdef test2
440! La version test2 destabilise beaucoup le modele.
441! Il semble donc que ca aide d'avoir un entrainement important sous
442! le nuage.
443!         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
444!          entr_star(ig,l)=0.4*detr_star(ig,l)
445!         else
446!          entr_star(ig,l)=0.
447!         endif
448! #else
449!
450! Deplacement du calcul de entr_star pour eviter d'avoir aussi
451! entr_star > fstar.
452! Redeplacer suite a la transformation du cas detr>f
453! FH
454
455        if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
456#define int1d
457!FH 070508 #define int1d4
458!#undef int1d4
459! L'option int1d4 correspond au choix dans le cas ou le detrainement
460! devient trop grand.
461
462#ifdef int1d
463
464#ifdef int1d4
465#else
466       detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
467!FH 070508 plus
468       detr_star(ig,l)=min(detr_star(ig,l),1.)
469#endif
470
471       entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
472
473        if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
474#ifdef int1d4
475! Si le detrainement excede le flux en bas + l'entrainement, le thermique
476! doit disparaitre.
477       if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
478          detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
479          f_star(ig,l+1)=0.
480          linter(ig)=l+1
481          zw2(ig,l+1)=-1.e-10
482       endif
483#endif
484
485
486#else
487
488        if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
489        if(l.gt.lalim(ig)) then
490         entr_star(ig,l)=0.4*detr_star(ig,l)
491        else
492
493! FH :
494! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
495! en haut de la couche d'alimentation.
496! A remettre en questoin a la premiere occasion mais ca peut aider a
497! ecrire un code robuste.
498! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
499! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
500! d'alimentation, ce qui n'est pas forcement heureux.
501
502        if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
503#undef pre_int1c
504#ifdef pre_int1c
505         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
506         detr_star(ig,l)=entr_star(ig,l)
507#else
508         entr_star(ig,l)=0.
509#endif
510
511        endif
512
513#endif
514
515        if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
516        entr_star1(ig,l)=entr_star(ig,l)
517        detr_star1(ig,l)=detr_star(ig,l)
518!
519
520#ifdef int1d
521#else
522        if (detr_star(ig,l).gt.f_star(ig,l)) then
523
524!  Ce test est là entre autres parce qu'on passe par des valeurs
525!  delirantes de detr_star.
526!  ca vaut sans doute le coup de verifier pourquoi.
527
528           detr_star(ig,l)=f_star(ig,l)
529#ifdef pre_int1c
530           if (l.gt.lalim(ig)+1) then
531               entr_star(ig,l)=0.
532               alim_star(ig,l)=0.
533! FH ajout pour forcer a stoper le thermique juste sous le sommet
534! de la couche (voir calcul de finter)
535               zw2(ig,l+1)=-1.e-10
536               linter(ig)=l+1
537            else
538               entr_star(ig,l)=0.4*detr_star(ig,l)
539            endif
540#else
541           entr_star(ig,l)=0.4*detr_star(ig,l)
542#endif
543        endif
544#endif
545
546      else !l > 2
547         detr_star(ig,l)=0.
548         entr_star(ig,l)=0.
549      endif
550
551        entr_star2(ig,l)=entr_star(ig,l)
552        detr_star2(ig,l)=detr_star(ig,l)
553        if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
554
555       endif  ! iflag_thermals_ed==0
556
557!CR:nvlle def de entr_star et detr_star
558      if (iflag_thermals_ed>=1) then
559!      if (l.lt.lalim(ig)) then
560!      if (l.lt.2) then
561!        entr_star(ig,l)=0.
562!        detr_star(ig,l)=0.
563!      else
564!      if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then
565!         entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
566!      else
567!         entr_star(ig,l)=  &
568!     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
569!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
570!     &                *(zlev(ig,l+1)-zlev(ig,l))
571
572 
573         entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &         
574     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
575     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)   &
576     &                *(zlev(ig,l+1)-zlev(ig,l))) &
577     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
578
579        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
580            alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
581            lalim(ig)=lmix_bis(ig)
582            if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
583        endif
584
585        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
586!        c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
587         c2(ig,l)=0.001
588         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
589     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
590     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
591     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
592     &                *(zlev(ig,l+1)-zlev(ig,l)))                    &
593     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
594
595       else
596!         c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
597          c2(ig,l)=0.003
598
599         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
600     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
601     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
602     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
603     &                *(zlev(ig,l+1)-zlev(ig,l))) &
604     &                +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
605       endif
606         
607           
608!        detr_star(ig,l)=detr_star(ig,l)*3.
609!        if (l.lt.lalim(ig)) then
610!          entr_star(ig,l)=0.
611!        endif
612!        if (l.lt.2) then
613!          entr_star(ig,l)=0.
614!          detr_star(ig,l)=0.
615!        endif
616
617
618!      endif
619!      else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
620!      entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1))        &
621!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
622!     &                *(zlev(ig,l+1)-zlev(ig,l))
623!      detr_star(ig,l)=0.002*f_star(ig,l)                         &
624!     &                *(zlev(ig,l+1)-zlev(ig,l))
625!      else
626!      entr_star(ig,l)=0.001*f_star(ig,l)                         &
627!     &                *(zlev(ig,l+1)-zlev(ig,l))
628!      detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1))       &
629!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))       &
630!     &                *(zlev(ig,l+1)-zlev(ig,l))                      &
631!     &                +0.002*f_star(ig,l)                             &
632!     &                *(zlev(ig,l+1)-zlev(ig,l))
633!      endif
634
635      endif   ! iflag_thermals_ed==1
636
637!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
638! FH inutile si on conserve comme on l'a fait plus haut entr=detr
639! dans la couche d'alimentation
640!pas d entrainement dans la couche alim
641!      if ((l.le.lalim(ig))) then
642!           entr_star(ig,l)=0.
643!      endif
644!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
645!
646!prise en compte du detrainement et de l entrainement dans le calcul du flux
647
648      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
649     &              -detr_star(ig,l)
650
651!test sur le signe de f_star
652        if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
653       if (f_star(ig,l+1).gt.1.e-10) then
654!----------------------------------------------------------------------------
655!calcul de la vitesse verticale en melangeant Tl et qt du thermique
656!---------------------------------------------------------------------------
657!
658       Zsat=.false.
659       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
660     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
661     &            /(f_star(ig,l+1)+detr_star(ig,l))
662!
663       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
664     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
665     &            /(f_star(ig,l+1)+detr_star(ig,l))
666
667               Tbef=ztla(ig,l)*zpspsk(ig,l)
668               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
669               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
670               qsatbef=MIN(0.5,qsatbef)
671               zcor=1./(1.-retv*qsatbef)
672               qsatbef=qsatbef*zcor
673               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
674               if (Zsat) then
675               qlbef=max(0.,zqta(ig,l)-qsatbef)
676               DT = 0.5*RLvCp*qlbef
677               do while (abs(DT).gt.DDT0)
678                 Tbef=Tbef+DT
679                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
680                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
681                 qsatbef=MIN(0.5,qsatbef)
682                 zcor=1./(1.-retv*qsatbef)
683                 qsatbef=qsatbef*zcor
684                 qlbef=zqta(ig,l)-qsatbef
685
686                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
687                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
688                 zcor=1./(1.-retv*qsatbef)
689                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
690                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
691                 denom=1.+RLvCp*dqsat_dT
692                 DT=num/denom
693              enddo
694                 zqla(ig,l) = max(0.,qlbef)
695              endif
696!   
697        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
698! on ecrit de maniere conservative (sat ou non)
699!          T = Tl +Lv/Cp ql
700           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
701           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
702!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
703           zha(ig,l) = ztva(ig,l)
704           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
705     &              -zqla(ig,l))-zqla(ig,l))
706
707!on ecrit zqsat
708           zqsatth(ig,l)=qsatbef 
709!calcul de vitesse
710           zw2(ig,l+1)=zw2(ig,l)*  &
711     &                 ((f_star(ig,l))**2)  &
712!  Tests de Catherine
713!     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
714     &      /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
715     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
716     &                 *fact_gamma &
717     &                 *(zlev(ig,l+1)-zlev(ig,l))
718!prise en compte des forces de pression que qd flottabilité<0
719!              zw2(ig,l+1)=zw2(ig,l)*  &
720!     &            1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &       
721!     &                 (f_star(ig,l))**2 &
722!     &                 /(f_star(ig,l)+entr_star(ig,l))**2+ &
723!     &                 (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+  &       
724!     &                 /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
725!     &                 /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
726!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
727!     &                 *1./3. &
728!     &                 *(zlev(ig,l+1)-zlev(ig,l))
729         
730!        write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
731!     &              -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
732!     &               2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
733
734 
735!             zw2(ig,l+1)=zw2(ig,l)*  &
736!     &                 (2.-2.*entr_star(ig,l)/f_star(ig,l)) & 
737!     &                 -zw2(ig,l-1)+  &       
738!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
739!     &                 *1./3. &
740!     &                 *(zlev(ig,l+1)-zlev(ig,l))             
741
742            endif
743        endif
744        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
745!
746!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
747
748            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
749!               stop'On tombe sur le cas particulier de thermcell_dry'
750                print*,'On tombe sur le cas particulier de thermcell_plume'
751                zw2(ig,l+1)=0.
752                linter(ig)=l+1
753            endif
754
755!        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
756        if (zw2(ig,l+1).lt.0.) then
757           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
758     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
759           zw2(ig,l+1)=0.
760        endif
761
762           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
763
764        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
765!   lmix est le niveau de la couche ou w (wa_moy) est maximum
766!on rajoute le calcul de lmix_bis
767            if (zqla(ig,l).lt.1.e-10) then
768               lmix_bis(ig)=l+1
769            endif
770            lmix(ig)=l+1
771            wmaxa(ig)=wa_moy(ig,l+1)
772        endif
773        enddo
774      enddo
775
776!on remplace a* par e* ds premiere couche
777!      if (iflag_thermals_ed.ge.1) then
778!       do ig=1,ngrid
779!       do l=2,klev
780!          if (l.lt.lalim(ig)) then
781!             alim_star(ig,l)=entr_star(ig,l)
782!          endif
783!       enddo
784!       enddo
785!       do ig=1,ngrid
786!          lalim(ig)=lmix_bis(ig)
787!       enddo
788!      endif
789       if (iflag_thermals_ed.ge.1) then
790          do ig=1,ngrid
791             do l=2,lalim(ig)
792                alim_star(ig,l)=entr_star(ig,l)
793                entr_star(ig,l)=0.
794             enddo
795           enddo
796       endif
797        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
798
799!     print*,'thermcell_plume OK'
800
801      return
802      end
Note: See TracBrowser for help on using the repository browser.