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

Last change on this file since 4093 was 4093, checked in by fhourdin, 2 years ago

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