source: LMDZ4/trunk/libf/phylmd/thermcell_plume.F90 @ 950

Last change on this file since 950 was 938, checked in by lmdzadmin, 16 years ago

Enleve prints par defaut
IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.4 KB
Line 
1      SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz,  &
2     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
3     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
4     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
5
6!--------------------------------------------------------------------------
7!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
8!--------------------------------------------------------------------------
9
10      IMPLICIT NONE
11
12#include "YOMCST.h"
13#include "YOETHF.h"
14#include "FCTTRE.h"
15#include "iniprint.h"
16
17      INTEGER ngrid,klev
18      REAL ztv(ngrid,klev)
19      REAL zthl(ngrid,klev)
20      REAL po(ngrid,klev)
21      REAL zl(ngrid,klev)
22      REAL rhobarz(ngrid,klev)
23      REAL zlev(ngrid,klev+1)
24      REAL pplev(ngrid,klev+1)
25      REAL pphi(ngrid,klev)
26      REAL zpspsk(ngrid,klev)
27      REAL alim_star(ngrid,klev)
28      REAL zmax_sec(ngrid)
29      REAL f0(ngrid)
30      REAL l_mix
31      REAL r_aspect
32      INTEGER lalim(ngrid)
33      integer lev_out                           ! niveau pour les print
34
35      REAL ztva(ngrid,klev)
36      REAL ztla(ngrid,klev)
37      REAL zqla(ngrid,klev)
38      REAL zqta(ngrid,klev)
39      REAL zha(ngrid,klev)
40
41      REAL detr_star(ngrid,klev)
42      REAL entr_star(ngrid,klev)
43      REAL detr(ngrid,klev)
44      REAL entr(ngrid,klev)
45
46      REAL zw2(ngrid,klev+1)
47      REAL w_est(ngrid,klev+1)
48      REAL f_star(ngrid,klev+1)
49      REAL wa_moy(ngrid,klev+1)
50
51      REAL ztva_est(ngrid,klev)
52      REAL zqla_est(ngrid,klev)
53      REAL zqsatth(ngrid,klev)
54
55      REAL linter(ngrid)
56      INTEGER lmix(ngrid)
57      REAL    wmaxa(ngrid)
58
59      INTEGER ig,l,k
60
61      real zcor,zdelta,zcvm5,qlbef
62      real Tbef,qsatbef
63      real dqsat_dT,DT,num,denom
64      REAL REPS,RLvCp,DDT0
65      PARAMETER (DDT0=.01)
66      logical Zsat
67
68      Zsat=.false.
69! Initialisation
70      RLvCp = RLVTT/RCPD
71     
72      do l=1,klev
73         do ig=1,ngrid
74            zqla_est(ig,l)=0.
75            ztva_est(ig,l)=ztva(ig,l)
76            zqsatth(ig,l)=0.
77         enddo
78      enddo
79
80!AM:initialisations du thermique
81      do k=1,klev
82         do ig=1,ngrid
83            ztva(ig,k)=ztv(ig,k)
84            ztla(ig,k)=zthl(ig,k)
85            zqla(ig,k)=0.
86            zqta(ig,k)=po(ig,k)
87         enddo
88      enddo
89      do k=1,klev
90        do ig=1,ngrid
91           detr_star(ig,k)=0.
92           entr_star(ig,k)=0.
93           detr(ig,k)=0.
94           entr(ig,k)=0.
95        enddo
96      enddo
97      if (prt_level.ge.1) print*,'7 OK convect8'
98      do k=1,klev+1
99         do ig=1,ngrid
100            zw2(ig,k)=0.
101            w_est(ig,k)=0.
102            f_star(ig,k)=0.
103            wa_moy(ig,k)=0.
104         enddo
105      enddo
106
107      if (prt_level.ge.1) print*,'8 OK convect8'
108      do ig=1,ngrid
109         linter(ig)=1.
110         lmix(ig)=1
111         wmaxa(ig)=0.
112      enddo
113
114!-----------------------------------------------------------------------------------
115!boucle de calcul de la vitesse verticale dans le thermique
116!-----------------------------------------------------------------------------------
117      do l=1,klev-1
118         do ig=1,ngrid
119            if (ztv(ig,l).gt.ztv(ig,l+1)  &
120     &         .and.alim_star(ig,l).gt.1.e-10  &
121     &         .and.zw2(ig,l).lt.1e-10) then
122               ztla(ig,l)=zthl(ig,l)
123               zqta(ig,l)=po(ig,l)
124               zqla(ig,l)=zl(ig,l)
125               f_star(ig,l+1)=alim_star(ig,l)
126               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
127     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
128     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
129               w_est(ig,l+1)=zw2(ig,l+1)
130!
131            else if ((zw2(ig,l).ge.1e-10).and.  &
132     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
133!estimation du detrainement a partir de la geometrie du pas precedent
134!tests sur la definition du detr
135!calcul de detr_star et entr_star
136
137
138
139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140! FH le test miraculeux de Catherine ? Le bout du tunel ?
141!               w_est(ig,3)=zw2(ig,2)*  &
142!    &                   ((f_star(ig,2))**2)  &
143!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
144!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
145!    &                   *(zlev(ig,3)-zlev(ig,2))
146!     if (l.gt.2) then
147!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
148
149      if (l.gt.2) then
150          w_est(ig,3)=zw2(ig,2)* &
151     &      ((f_star(ig,2))**2) &
152     &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
153     &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
154     &      *(zlev(ig,3)-zlev(ig,2))
155
156
157!---------------------------------------------------------------------------
158!calcul de l entrainement et du detrainement lateral
159!---------------------------------------------------------------------------
160!
161!test:estimation de ztva_new_est sans entrainement
162               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
163               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
164               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
165               qsatbef=MIN(0.5,qsatbef)
166               zcor=1./(1.-retv*qsatbef)
167               qsatbef=qsatbef*zcor
168               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
169               if (Zsat) then
170               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
171               DT = 0.5*RLvCp*qlbef
172               do while (abs(DT).gt.DDT0)
173                 Tbef=Tbef+DT
174                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
175                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
176                 qsatbef=MIN(0.5,qsatbef)
177                 zcor=1./(1.-retv*qsatbef)
178                 qsatbef=qsatbef*zcor
179                 qlbef=zqta(ig,l-1)-qsatbef
180
181                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
182                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
183                 zcor=1./(1.-retv*qsatbef)
184                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
185                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
186                 denom=1.+RLvCp*dqsat_dT
187                 DT=num/denom
188               enddo
189                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
190               endif
191        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
192        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
193        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
194     &      -zqla_est(ig,l))-zqla_est(ig,l))
195
196             w_est(ig,l+1)=zw2(ig,l)*  &
197     &                   ((f_star(ig,l))**2)  &
198     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
199     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
200     &                   *(zlev(ig,l+1)-zlev(ig,l))
201             if (w_est(ig,l+1).lt.0.) then
202                w_est(ig,l+1)=zw2(ig,l)
203             endif
204!
205!calcul du detrainement
206          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
207     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
208     &       (zqla(ig,l-1).lt.1.e-10)) then
209             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
210     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
211     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
212     &       /(r_aspect*zmax_sec(ig)))
213       if (prt_level.ge.20) print*,'coucou calcul detr 1'
214          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
215     &            (zqla(ig,l-1).lt.1.e-10)) then
216             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
217     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
218     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
219     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
220     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
221     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
222     &       *((zmax_sec(ig)-zlev(ig,l))/  &
223     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
224        if (prt_level.ge.20) print*,'coucou calcul detr 2'
225          else
226             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
227     &                      *(zlev(ig,l+1)-zlev(ig,l))
228        if (prt_level.ge.20) print*,'coucou calcul detr 3'
229             
230          endif
231          detr_star(ig,l)=detr_star(ig,l)/f0(ig)
232!
233!calcul de entr_star
234!
235! Deplacement du calcul de entr_star pour eviter d'avoir aussi
236! entr_star > fstar.
237! Redeplacer suite a la transformation du cas detr>f
238! FH
239        if(l.gt.lalim(ig)) then
240         entr_star(ig,l)=0.4*detr_star(ig,l)
241        else
242
243! FH :
244! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
245! en haut de la couche d'alimentation.
246! A remettre en questoin a la premiere occasion mais ca peut aider a
247! ecrire un code robuste.
248! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
249! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
250! d'alimentation, ce qui n'est pas forcement heureux.
251         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
252         detr_star(ig,l)=entr_star(ig,l)
253        endif
254
255!
256        if (detr_star(ig,l).gt.f_star(ig,l)) then
257
258!  Ce test est là entre autres parce qu'on passe par des valeurs
259!  delirantes de detr_star.
260!  ca vaut sans doute le coup de verifier pourquoi.
261
262           detr_star(ig,l)=f_star(ig,l)
263           if (l.gt.lalim(ig)+1) then
264               entr_star(ig,l)=0.
265               alim_star(ig,l)=0.
266! FH ajout pour forcer a stoper le thermique juste sous le sommet
267! de la couche (voir calcul de finter)
268               zw2(ig,l+1)=-1.e-10
269               linter(ig)=l+1
270            else
271               entr_star(ig,l)=detr_star(ig,l)
272            endif
273        endif
274
275      else
276         detr_star(ig,l)=0.
277         entr_star(ig,l)=0.
278      endif
279
280!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281! FH inutile si on conserve comme on l'a fait plus haut entr=detr
282! dans la couche d'alimentation
283!pas d entrainement dans la couche alim
284!      if ((l.le.lalim(ig))) then
285!           entr_star(ig,l)=0.
286!      endif
287!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288!
289!prise en compte du detrainement et de l entrainement dans le calcul du flux
290      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
291     &              -detr_star(ig,l)
292
293!test sur le signe de f_star
294       if (f_star(ig,l+1).gt.1.e-10) then
295!----------------------------------------------------------------------------
296!calcul de la vitesse verticale en melangeant Tl et qt du thermique
297!---------------------------------------------------------------------------
298!
299       Zsat=.false.
300       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
301     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
302     &            /(f_star(ig,l+1)+detr_star(ig,l))
303!
304       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
305     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
306     &            /(f_star(ig,l+1)+detr_star(ig,l))
307
308               Tbef=ztla(ig,l)*zpspsk(ig,l)
309               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
310               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
311               qsatbef=MIN(0.5,qsatbef)
312               zcor=1./(1.-retv*qsatbef)
313               qsatbef=qsatbef*zcor
314               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
315               if (Zsat) then
316               qlbef=max(0.,zqta(ig,l)-qsatbef)
317               DT = 0.5*RLvCp*qlbef
318               do while (abs(DT).gt.DDT0)
319                 Tbef=Tbef+DT
320                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
321                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
322                 qsatbef=MIN(0.5,qsatbef)
323                 zcor=1./(1.-retv*qsatbef)
324                 qsatbef=qsatbef*zcor
325                 qlbef=zqta(ig,l)-qsatbef
326
327                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
328                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
329                 zcor=1./(1.-retv*qsatbef)
330                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
331                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
332                 denom=1.+RLvCp*dqsat_dT
333                 DT=num/denom
334              enddo
335                 zqla(ig,l) = max(0.,qlbef)
336              endif
337!   
338! on ecrit de maniere conservative (sat ou non)
339!          T = Tl +Lv/Cp ql
340           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
341           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
342!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
343           zha(ig,l) = ztva(ig,l)
344           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
345     &              -zqla(ig,l))-zqla(ig,l))
346
347!on ecrit zqsat
348           zqsatth(ig,l)=qsatbef 
349!calcul de vitesse
350           zw2(ig,l+1)=zw2(ig,l)*  &
351     &                 ((f_star(ig,l))**2)  &
352     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
353     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
354     &                 *(zlev(ig,l+1)-zlev(ig,l))
355               
356            endif
357        endif
358!
359!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
360
361            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
362!               stop'On tombe sur le cas particulier de thermcell_dry'
363                print*,'On tombe sur le cas particulier de thermcell_plume'
364                zw2(ig,l+1)=0.
365                linter(ig)=l+1
366            endif
367
368
369        if (zw2(ig,l+1).lt.0.) then
370           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
371     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
372           zw2(ig,l+1)=0.
373        endif
374
375           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
376
377        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
378!   lmix est le niveau de la couche ou w (wa_moy) est maximum
379            lmix(ig)=l+1
380            wmaxa(ig)=wa_moy(ig,l+1)
381        endif
382        enddo
383      enddo
384
385
386      return
387      end
Note: See TracBrowser for help on using the repository browser.