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

Last change on this file was 5434, checked in by fhourdin, 15 hours ago

Superessing CPP in lmdz_*

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