source: LMDZ6/trunk/libf/phylmd/thermcell_plume.F90 @ 4305

Last change on this file since 4305 was 4162, checked in by fhourdin, 3 years ago

Bug fix
File size: 16.3 KB
Line 
1!
2! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
3!
4      SUBROUTINE thermcell_plume(itap,ngrid,nlay,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!     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
10!--------------------------------------------------------------------------
11! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
12!
13!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
14!   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
15!   thermcell_plume_6A is activate for flag_thermas_ed < 10
16!   thermcell_plume_5B for flag_thermas_ed < 20
17!   thermcell_plume for flag_thermals_ed>= 20
18!   Various options are controled by the flag_thermals_ed parameter
19!   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
20!   = 21 : the Jam strato-cumulus modif is not activated in detrainment
21!   = 29 : an other way to compute the modified buoyancy (to be tested)
22!--------------------------------------------------------------------------
23
24       USE thermcell_ini_mod, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
25       USE thermcell_ini_mod, ONLY: fact_epsilon, betalpha, afact, fact_shell
26       USE thermcell_ini_mod, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
27       USE thermcell_ini_mod, ONLY: mix0, thermals_flag_alim
28
29       IMPLICIT NONE
30
31      integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay
32      real,intent(in) :: ptimestep
33      real,intent(in),dimension(ngrid,nlay) :: ztv
34      real,intent(in),dimension(ngrid,nlay) :: zthl
35      real,intent(in),dimension(ngrid,nlay) :: po
36      real,intent(in),dimension(ngrid,nlay) :: zl
37      real,intent(in),dimension(ngrid,nlay) :: rhobarz
38      real,intent(in),dimension(ngrid,nlay+1) :: zlev
39      real,intent(in),dimension(ngrid,nlay+1) :: pplev
40      real,intent(in),dimension(ngrid,nlay) :: pphi
41      real,intent(in),dimension(ngrid,nlay) :: zpspsk
42      real,intent(in),dimension(ngrid) :: f0
43
44      integer,intent(out) :: lalim(ngrid)
45      real,intent(out),dimension(ngrid,nlay) :: alim_star
46      real,intent(out),dimension(ngrid) :: alim_star_tot
47      real,intent(out),dimension(ngrid,nlay) :: detr_star
48      real,intent(out),dimension(ngrid,nlay) :: entr_star
49      real,intent(out),dimension(ngrid,nlay+1) :: f_star
50      real,intent(out),dimension(ngrid,nlay) :: csc
51      real,intent(out),dimension(ngrid,nlay) :: ztva
52      real,intent(out),dimension(ngrid,nlay) :: ztla
53      real,intent(out),dimension(ngrid,nlay) :: zqla
54      real,intent(out),dimension(ngrid,nlay) :: zqta
55      real,intent(out),dimension(ngrid,nlay) :: zha
56      real,intent(out),dimension(ngrid,nlay+1) :: zw2
57      real,intent(out),dimension(ngrid,nlay+1) :: w_est
58      real,intent(out),dimension(ngrid,nlay) :: ztva_est
59      real,intent(out),dimension(ngrid,nlay) :: zqsatth
60      integer,intent(out),dimension(ngrid) :: lmix(ngrid)
61      integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid)
62      real,intent(out),dimension(ngrid) :: linter(ngrid)
63
64
65      REAL,dimension(ngrid,nlay+1) :: wa_moy
66      REAL,dimension(ngrid,nlay) :: entr,detr
67      REAL,dimension(ngrid,nlay) :: ztv_est
68      REAL,dimension(ngrid,nlay) :: zqla_est
69      REAL,dimension(ngrid,nlay) :: zta_est
70      REAL,dimension(ngrid) :: ztemp,zqsat
71      REAL zdw2,zdw2bis
72      REAL zw2modif
73      REAL zw2fact,zw2factbis
74      REAL,dimension(ngrid,nlay) :: zeps
75
76      REAL,dimension(ngrid) :: wmaxa
77
78      INTEGER ig,l,k,lt,it,lm,nbpb
79
80      real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
81      real zdz,zalpha,zw2m
82      real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
83      real zdz2,zdz3,lmel,entrbis,zdzbis
84      real,dimension(ngrid) :: d_temp
85      real ztv1,ztv2,factinv,zinv,zlmel
86      real zlmelup,zlmeldwn,zlt,zltdwn,zltup
87      real atv1,atv2,btv1,btv2
88      real ztv_est1,ztv_est2
89      real zcor,zdelta,zcvm5,qlbef
90      real zbetalpha, coefzlmel
91      real eps
92      logical Zsat
93      LOGICAL,dimension(ngrid) :: active,activetmp
94      REAL fact_gamma,fact_gamma2,fact_epsilon2
95
96
97      REAL,dimension(ngrid,nlay) :: c2
98
99      if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
100      Zsat=.false.
101! Initialisation
102
103
104      zbetalpha=betalpha/(1.+betalpha)
105
106
107! Initialisations des variables r?elles
108if (1==1) then
109      ztva(:,:)=ztv(:,:)
110      ztva_est(:,:)=ztva(:,:)
111      ztv_est(:,:)=ztv(:,:)
112      ztla(:,:)=zthl(:,:)
113      zqta(:,:)=po(:,:)
114      zqla(:,:)=0.
115      zha(:,:) = ztva(:,:)
116else
117      ztva(:,:)=0.
118      ztv_est(:,:)=0.
119      ztva_est(:,:)=0.
120      ztla(:,:)=0.
121      zqta(:,:)=0.
122      zha(:,:) =0.
123endif
124
125      zqla_est(:,:)=0.
126      zqsatth(:,:)=0.
127      zqla(:,:)=0.
128      detr_star(:,:)=0.
129      entr_star(:,:)=0.
130      alim_star(:,:)=0.
131      alim_star_tot(:)=0.
132      csc(:,:)=0.
133      detr(:,:)=0.
134      entr(:,:)=0.
135      zw2(:,:)=0.
136      zbuoy(:,:)=0.
137      zbuoyjam(:,:)=0.
138      gamma(:,:)=0.
139      zeps(:,:)=0.
140      w_est(:,:)=0.
141      f_star(:,:)=0.
142      wa_moy(:,:)=0.
143      linter(:)=1.
144!     linter(:)=1.
145! Initialisation des variables entieres
146      lmix(:)=1
147      lmix_bis(:)=2
148      wmaxa(:)=0.
149
150
151!-------------------------------------------------------------------------
152! On ne considere comme actif que les colonnes dont les deux premieres
153! couches sont instables.
154!-------------------------------------------------------------------------
155
156      active(:)=ztv(:,1)>ztv(:,2)
157      d_temp(:)=0. ! Pour activer un contraste de temperature a la base
158                   ! du panache
159!  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
160      CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
161
162!------------------------------------------------------------------------------
163! Calcul dans la premiere couche
164! On decide dans cette version que le thermique n'est actif que si la premiere
165! couche est instable.
166! Pourrait etre change si on veut que le thermiques puisse se d??clencher
167! dans une couche l>1
168!------------------------------------------------------------------------------
169do ig=1,ngrid
170! Le panache va prendre au debut les caracteristiques de l'air contenu
171! dans cette couche.
172    if (active(ig)) then
173    ztla(ig,1)=zthl(ig,1)
174    zqta(ig,1)=po(ig,1)
175    zqla(ig,1)=zl(ig,1)
176!cr: attention, prise en compte de f*(1)=1
177    f_star(ig,2)=alim_star(ig,1)
178    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
179&                     *(zlev(ig,2)-zlev(ig,1))  &
180&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
181    w_est(ig,2)=zw2(ig,2)
182    endif
183enddo
184!
185
186!==============================================================================
187!boucle de calcul de la vitesse verticale dans le thermique
188!==============================================================================
189do l=2,nlay-1
190!==============================================================================
191
192
193! On decide si le thermique est encore actif ou non
194! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
195    do ig=1,ngrid
196       active(ig)=active(ig) &
197&                 .and. zw2(ig,l)>1.e-10 &
198&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
199    enddo
200
201
202
203!---------------------------------------------------------------------------
204! calcul des proprietes thermodynamiques et de la vitesse de la couche l
205! sans tenir compte du detrainement et de l'entrainement dans cette
206! couche
207! C'est a dire qu'on suppose
208! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
209! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
210! avant) a l'alimentation pour avoir un calcul plus propre
211!---------------------------------------------------------------------------
212
213   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
214   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
215    do ig=1,ngrid
216!       print*,'active',active(ig),ig,l
217        if(active(ig)) then
218        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
219        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
220        zta_est(ig,l)=ztva_est(ig,l)
221        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
222        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
223     &      -zqla_est(ig,l))-zqla_est(ig,l))
224 
225
226!Modif AJAM
227
228        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
229        zdz=zlev(ig,l+1)-zlev(ig,l)         
230        lmel=fact_thermals_ed_dz*zlev(ig,l)
231!        lmel=0.09*zlev(ig,l)
232        zlmel=zlev(ig,l)+lmel
233        zlmelup=zlmel+(zdz/2)
234        zlmeldwn=zlmel-(zdz/2)
235
236        lt=l+1
237        zlt=zlev(ig,lt)
238        zdz3=zlev(ig,lt+1)-zlt
239        zltdwn=zlt-zdz3/2
240        zltup=zlt+zdz3/2
241         
242!=========================================================================
243! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
244!=========================================================================
245
246!--------------------------------------------------
247        lt=l+1
248        zlt=zlev(ig,lt)
249        zdz2=zlev(ig,lt)-zlev(ig,l)
250
251        do while (lmel.gt.zdz2)
252           lt=lt+1
253           zlt=zlev(ig,lt)
254           zdz2=zlt-zlev(ig,l)
255        enddo
256        zdz3=zlev(ig,lt+1)-zlt
257        zltdwn=zlev(ig,lt)-zdz3/2
258        zlmelup=zlmel+(zdz/2)
259        coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
260        zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- &
261    &   ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
262    &   ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
263
264!------------------------------------------------
265!AJAM:nouveau calcul de w? 
266!------------------------------------------------
267        zdz=zlev(ig,l+1)-zlev(ig,l)
268        zdzbis=zlev(ig,l)-zlev(ig,l-1)
269        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
270        zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
271        zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
272        zdw2=afact*zbuoy(ig,l)/fact_epsilon
273        zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
274        lm=Max(1,l-2)
275        w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
276       endif
277    enddo
278
279
280!-------------------------------------------------
281!calcul des taux d'entrainement et de detrainement
282!-------------------------------------------------
283
284     do ig=1,ngrid
285        if (active(ig)) then
286
287!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
288          zw2m=w_est(ig,l+1)
289          zdz=zlev(ig,l+1)-zlev(ig,l)
290          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
291          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
292          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
293
294!=========================================================================
295! 4. Calcul de l'entrainement et du detrainement
296!=========================================================================
297
298          detr_star(ig,l)=f_star(ig,l)*zdz             &
299    &     *( mix0 * 0.1 / (zalpha+0.001)               &
300    &     + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
301    &     + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
302
303          if ( iflag_thermals_ed == 20 ) then
304             entr_star(ig,l)=f_star(ig,l)*zdz* (         &
305    &          mix0 * 0.1 / (zalpha+0.001)               &
306    &        + zbetalpha*MAX(entr_min,                   &
307    &        afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
308          else
309             entr_star(ig,l)=f_star(ig,l)*zdz* (         &
310    &          mix0 * 0.1 / (zalpha+0.001)               &
311    &        + zbetalpha*MAX(entr_min,                   &
312    &        afact*zbuoy(ig,l)/zw2m - fact_epsilon))
313          endif
314         
315! En dessous de lalim, on prend le max de alim_star et entr_star pour
316! alim_star et 0 sinon
317        if (l.lt.lalim(ig)) then
318          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
319          entr_star(ig,l)=0.
320        endif
321        f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
322     &              -detr_star(ig,l)
323
324      endif
325   enddo
326
327
328!============================================================================
329! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
330!===========================================================================
331
332   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
333   do ig=1,ngrid
334       if (activetmp(ig)) then
335           Zsat=.false.
336           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
337     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
338     &            /(f_star(ig,l+1)+detr_star(ig,l))
339           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
340     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
341     &            /(f_star(ig,l+1)+detr_star(ig,l))
342
343        endif
344    enddo
345
346   ztemp(:)=zpspsk(:,l)*ztla(:,l)
347   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
348   do ig=1,ngrid
349      if (activetmp(ig)) then
350! on ecrit de maniere conservative (sat ou non)
351!          T = Tl +Lv/Cp ql
352           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
353           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
354           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
355!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
356           zha(ig,l) = ztva(ig,l)
357           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
358     &              -zqla(ig,l))-zqla(ig,l))
359           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
360           zdz=zlev(ig,l+1)-zlev(ig,l)
361           zdzbis=zlev(ig,l)-zlev(ig,l-1)
362           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
363           zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
364           zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
365           zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
366           zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
367           zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
368      endif
369   enddo
370
371   if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
372!
373!===========================================================================
374! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
375!===========================================================================
376
377   nbpb=0
378   do ig=1,ngrid
379            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
380!               stop'On tombe sur le cas particulier de thermcell_dry'
381!               print*,'On tombe sur le cas particulier de thermcell_plume'
382                nbpb=nbpb+1
383                zw2(ig,l+1)=0.
384                linter(ig)=l+1
385            endif
386
387        if (zw2(ig,l+1).lt.0.) then
388           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
389     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
390           zw2(ig,l+1)=0.
391!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
392        elseif (f_star(ig,l+1).lt.0.) then
393           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
394     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
395           zw2(ig,l+1)=0.
396!fin CR:04/05/12
397        endif
398
399           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
400
401        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
402!   lmix est le niveau de la couche ou w (wa_moy) est maximum
403!on rajoute le calcul de lmix_bis
404            if (zqla(ig,l).lt.1.e-10) then
405               lmix_bis(ig)=l+1
406            endif
407            lmix(ig)=l+1
408            wmaxa(ig)=wa_moy(ig,l+1)
409        endif
410   enddo
411
412   if (nbpb>0) then
413   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
414   endif
415
416!=========================================================================
417! FIN DE LA BOUCLE VERTICALE
418      enddo
419!=========================================================================
420
421!on recalcule alim_star_tot
422       do ig=1,ngrid
423          alim_star_tot(ig)=0.
424       enddo
425       do ig=1,ngrid
426          do l=1,lalim(ig)-1
427          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
428          enddo
429       enddo
430       
431
432        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
433
434#undef wrgrads_thermcell
435#ifdef wrgrads_thermcell
436         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
437         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
438         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
439         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
440         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
441         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
442         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
443#endif
444
445
446 RETURN
447     end
Note: See TracBrowser for help on using the repository browser.