source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume_5B.f90

Last change on this file was 5852, checked in by fhourdin, 7 weeks ago

Separation thermcell_plume 5B/6A
File size: 12.7 KB
Line 
1MODULE lmdz_thermcell_plume_5B
2CONTAINS
3 SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
4&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
5&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
6&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
7&           ,lev_out,lunout1,igout)
8!&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
9
10!--------------------------------------------------------------------------
11!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12! Version conforme a l'article de Rio et al. 2010.
13! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
14!--------------------------------------------------------------------------
15
16      USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
17       USE lmdz_thermcell_qsat, ONLY : thermcell_qsat
18      IMPLICIT NONE
19
20      INTEGER itap
21      INTEGER lunout1,igout
22      INTEGER ngrid,nlay
23      REAL ptimestep
24      REAL ztv(ngrid,nlay)
25      REAL zthl(ngrid,nlay)
26      REAL, INTENT(IN) :: po(ngrid,nlay)
27      REAL zl(ngrid,nlay)
28      REAL rhobarz(ngrid,nlay)
29      REAL zlev(ngrid,nlay+1)
30      REAL pplev(ngrid,nlay+1)
31      REAL pphi(ngrid,nlay)
32      REAL zpspsk(ngrid,nlay)
33      REAL alim_star(ngrid,nlay)
34      REAL f0(ngrid)
35      INTEGER lalim(ngrid)
36      integer lev_out                           ! niveau pour les print
37      integer nbpb
38   
39      real alim_star_tot(ngrid)
40
41      REAL ztva(ngrid,nlay)
42      REAL ztla(ngrid,nlay)
43      REAL zqla(ngrid,nlay)
44      REAL zqta(ngrid,nlay)
45      REAL zha(ngrid,nlay)
46
47      REAL detr_star(ngrid,nlay)
48      REAL coefc
49      REAL entr_star(ngrid,nlay)
50      REAL detr(ngrid,nlay)
51      REAL entr(ngrid,nlay)
52
53      REAL csc(ngrid,nlay)
54
55      REAL zw2(ngrid,nlay+1)
56      REAL w_est(ngrid,nlay+1)
57      REAL f_star(ngrid,nlay+1)
58      REAL wa_moy(ngrid,nlay+1)
59
60      REAL ztva_est(ngrid,nlay)
61      REAL zqla_est(ngrid,nlay)
62      REAL zqsatth(ngrid,nlay)
63      REAL zta_est(ngrid,nlay)
64      REAL zbuoyjam(ngrid,nlay)
65      REAL ztemp(ngrid),zqsat(ngrid)
66      REAL zdw2
67      REAL zw2modif
68      REAL zw2fact
69      REAL zeps(ngrid,nlay)
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,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
79      real zbuoybis
80      real zcor,zdelta,zcvm5,qlbef,zdz2
81      real betalpha,zbetalpha
82      real eps, afact
83      logical Zsat
84      LOGICAL active(ngrid),activetmp(ngrid)
85      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
86      REAL c2(ngrid,nlay)
87      Zsat=.false.
88! Initialisation
89
90      fact_epsilon=0.002
91      betalpha=0.9
92      afact=2./3.           
93
94      zbetalpha=betalpha/(1.+betalpha)
95
96
97! Initialisations des variables reeles
98if (1==1) then
99      ztva(:,:)=ztv(:,:)
100      ztva_est(:,:)=ztva(:,:)
101      ztla(:,:)=zthl(:,:)
102      zqta(:,:)=po(:,:)
103      zha(:,:) = ztva(:,:)
104else
105      ztva(:,:)=0.
106      ztva_est(:,:)=0.
107      ztla(:,:)=0.
108      zqta(:,:)=0.
109      zha(:,:) =0.
110endif
111
112      zqla_est(:,:)=0.
113      zqsatth(:,:)=0.
114      zqla(:,:)=0.
115      detr_star(:,:)=0.
116      entr_star(:,:)=0.
117      alim_star(:,:)=0.
118      alim_star_tot(:)=0.
119      csc(:,:)=0.
120      detr(:,:)=0.
121      entr(:,:)=0.
122      zw2(:,:)=0.
123      zbuoy(:,:)=0.
124      zbuoyjam(:,:)=0.
125      gamma(:,:)=0.
126      zeps(:,:)=0.
127      w_est(:,:)=0.
128      f_star(:,:)=0.
129      wa_moy(:,:)=0.
130      linter(:)=1.
131!     linter(:)=1.
132! Initialisation des variables entieres
133      lmix(:)=1
134      lmix_bis(:)=2
135      wmaxa(:)=0.
136      lalim(:)=1
137
138
139!-------------------------------------------------------------------------
140! On ne considere comme actif que les colonnes dont les deux premieres
141! couches sont instables.
142!-------------------------------------------------------------------------
143      active(:)=ztv(:,1)>ztv(:,2)
144
145!-------------------------------------------------------------------------
146! Definition de l'alimentation
147!-------------------------------------------------------------------------
148      do l=1,nlay-1
149         do ig=1,ngrid
150            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
151               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
152     &                       *sqrt(zlev(ig,l+1))
153               lalim(ig)=l+1
154               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
155            endif
156         enddo
157      enddo
158      do l=1,nlay
159         do ig=1,ngrid
160            if (alim_star_tot(ig) > 1.e-10 ) then
161               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
162            endif
163         enddo
164      enddo
165      alim_star_tot(:)=1.
166
167
168
169!------------------------------------------------------------------------------
170! Calcul dans la premiere couche
171! On decide dans cette version que le thermique n'est actif que si la premiere
172! couche est instable.
173! Pourrait etre change si on veut que le thermiques puisse se d??clencher
174! dans une couche l>1
175!------------------------------------------------------------------------------
176do ig=1,ngrid
177! Le panache va prendre au debut les caracteristiques de l'air contenu
178! dans cette couche.
179    if (active(ig)) then
180    ztla(ig,1)=zthl(ig,1)
181    zqta(ig,1)=po(ig,1)
182    zqla(ig,1)=zl(ig,1)
183!cr: attention, prise en compte de f*(1)=1
184    f_star(ig,2)=alim_star(ig,1)
185    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
186&                     *(zlev(ig,2)-zlev(ig,1))  &
187&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
188    w_est(ig,2)=zw2(ig,2)
189    endif
190enddo
191!
192
193!==============================================================================
194!boucle de calcul de la vitesse verticale dans le thermique
195!==============================================================================
196do l=2,nlay-1
197!==============================================================================
198
199
200! On decide si le thermique est encore actif ou non
201! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
202    do ig=1,ngrid
203       active(ig)=active(ig) &
204&                 .and. zw2(ig,l)>1.e-10 &
205&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
206    enddo
207
208
209
210!---------------------------------------------------------------------------
211! calcul des proprietes thermodynamiques et de la vitesse de la couche l
212! sans tenir compte du detrainement et de l'entrainement dans cette
213! couche
214! C'est a dire qu'on suppose
215! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
216! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
217! avant) a l'alimentation pour avoir un calcul plus propre
218!---------------------------------------------------------------------------
219
220   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
221   call thermcell_qsat(ngrid, 1, active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
222
223    do ig=1,ngrid
224!       print*,'active',active(ig),ig,l
225        if(active(ig)) then
226        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
227        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
228        zta_est(ig,l)=ztva_est(ig,l)
229        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
230        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
231     &      -zqla_est(ig,l))-zqla_est(ig,l))
232
233!------------------------------------------------
234!AJAM:nouveau calcul de w? 
235!------------------------------------------------
236              zdz=zlev(ig,l+1)-zlev(ig,l)
237              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
238
239              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
240              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
241              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
242 
243
244             if (w_est(ig,l+1).lt.0.) then
245                w_est(ig,l+1)=zw2(ig,l)
246             endif
247       endif
248    enddo
249
250
251!-------------------------------------------------
252!calcul des taux d'entrainement et de detrainement
253!-------------------------------------------------
254
255     do ig=1,ngrid
256        if (active(ig)) then
257
258          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
259          zw2m=w_est(ig,l+1)
260          zdz=zlev(ig,l+1)-zlev(ig,l)
261          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
262!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
263          zbuoybis=zbuoy(ig,l)
264          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
265          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
266
267         
268          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
269    &     afact*zbuoybis/zw2m - fact_epsilon )
270
271
272          detr_star(ig,l)=f_star(ig,l)*zdz                        &
273    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
274    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
275         
276! En dessous de lalim, on prend le max de alim_star et entr_star pour
277! alim_star et 0 sinon
278        if (l.lt.lalim(ig)) then
279          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
280          entr_star(ig,l)=0.
281        endif
282
283! Calcul du flux montant normalise
284      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
285     &              -detr_star(ig,l)
286
287      endif
288   enddo
289
290
291!----------------------------------------------------------------------------
292!calcul de la vitesse verticale en melangeant Tl et qt du thermique
293!---------------------------------------------------------------------------
294   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
295   do ig=1,ngrid
296       if (activetmp(ig)) then
297           Zsat=.false.
298           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
299     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
300     &            /(f_star(ig,l+1)+detr_star(ig,l))
301           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
302     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
303     &            /(f_star(ig,l+1)+detr_star(ig,l))
304
305        endif
306    enddo
307
308   ztemp(:)=zpspsk(:,l)*ztla(:,l)
309   call thermcell_qsat(ngrid, 1, activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
310
311   do ig=1,ngrid
312      if (activetmp(ig)) then
313! on ecrit de maniere conservative (sat ou non)
314!          T = Tl +Lv/Cp ql
315           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
316           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
317           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
318!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
319           zha(ig,l) = ztva(ig,l)
320           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
321     &              -zqla(ig,l))-zqla(ig,l))
322           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
323           zdz=zlev(ig,l+1)-zlev(ig,l)
324           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
325
326            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
327            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
328            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
329      endif
330   enddo
331
332        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
333!
334!---------------------------------------------------------------------------
335!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
336!---------------------------------------------------------------------------
337
338   nbpb=0
339   do ig=1,ngrid
340            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
341!               stop'On tombe sur le cas particulier de thermcell_dry'
342!               print*,'On tombe sur le cas particulier de thermcell_plume'
343                nbpb=nbpb+1
344                zw2(ig,l+1)=0.
345                linter(ig)=l+1
346            endif
347
348        if (zw2(ig,l+1).lt.0.) then
349           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
350     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
351           zw2(ig,l+1)=0.
352        elseif (f_star(ig,l+1).lt.0.) then
353           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
354     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
355!           print*,"linter plume", linter(ig)
356           zw2(ig,l+1)=0.
357        endif
358
359           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
360
361        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
362!   lmix est le niveau de la couche ou w (wa_moy) est maximum
363!on rajoute le calcul de lmix_bis
364            if (zqla(ig,l).lt.1.e-10) then
365               lmix_bis(ig)=l+1
366            endif
367            lmix(ig)=l+1
368            wmaxa(ig)=wa_moy(ig,l+1)
369        endif
370   enddo
371
372   if (nbpb>0) then
373   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
374   endif
375
376!=========================================================================
377! FIN DE LA BOUCLE VERTICALE
378      enddo
379!=========================================================================
380
381!on recalcule alim_star_tot
382       do ig=1,ngrid
383          alim_star_tot(ig)=0.
384       enddo
385       do ig=1,ngrid
386          do l=1,lalim(ig)-1
387          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
388          enddo
389       enddo
390       
391
392        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
393
394     return
395     END SUBROUTINE thermcell_plume_5B
396END MODULE lmdz_thermcell_plume_5B
Note: See TracBrowser for help on using the repository browser.