source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume.F90 @ 5099

Last change on this file since 5099 was 5099, checked in by abarral, 2 months ago

Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90
Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad
Create lmdz_yoerad.f90
(lint) Remove "!" on otherwise empty line
File size: 16.4 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!boucle de calcul de la vitesse verticale dans le thermique
193!==============================================================================
194do l=2,nlay-1
195!==============================================================================
196
197
198! On decide si le thermique est encore actif ou non
199! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
200    do ig=1,ngrid
201       active(ig)=active(ig) &
202&                 .and. zw2(ig,l)>1.e-10 &
203&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
204    enddo
205
206
207
208!---------------------------------------------------------------------------
209! calcul des proprietes thermodynamiques et de la vitesse de la couche l
210! sans tenir compte du detrainement et de l'entrainement dans cette
211! couche
212! C'est a dire qu'on suppose
213! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
214! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
215! avant) a l'alimentation pour avoir un calcul plus propre
216!---------------------------------------------------------------------------
217
218   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
219   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
220    do ig=1,ngrid
221!       print*,'active',active(ig),ig,l
222        if(active(ig)) then
223        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
224        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
225        zta_est(ig,l)=ztva_est(ig,l)
226        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
227        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
228        -zqla_est(ig,l))-zqla_est(ig,l))
229 
230
231!Modif AJAM
232
233        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
234        zdz=zlev(ig,l+1)-zlev(ig,l)         
235        lmel=fact_thermals_ed_dz*zlev(ig,l)
236!        lmel=0.09*zlev(ig,l)
237        zlmel=zlev(ig,l)+lmel
238        zlmelup=zlmel+(zdz/2)
239        zlmeldwn=zlmel-(zdz/2)
240
241        lt=l+1
242        zlt=zlev(ig,lt)
243        zdz3=zlev(ig,lt+1)-zlt
244        zltdwn=zlt-zdz3/2
245        zltup=zlt+zdz3/2
246         
247!=========================================================================
248! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
249!=========================================================================
250
251!--------------------------------------------------
252        lt=l+1
253        zlt=zlev(ig,lt)
254        zdz2=zlev(ig,lt)-zlev(ig,l)
255
256        do while (lmel>zdz2)
257           lt=lt+1
258           zlt=zlev(ig,lt)
259           zdz2=zlt-zlev(ig,l)
260        enddo
261        zdz3=zlev(ig,lt+1)-zlt
262        zltdwn=zlev(ig,lt)-zdz3/2
263        zlmelup=zlmel+(zdz/2)
264        coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
265        zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- &
266     ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
267     ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
268
269!------------------------------------------------
270!AJAM:nouveau calcul de w? 
271!------------------------------------------------
272        zdz=zlev(ig,l+1)-zlev(ig,l)
273        zdzbis=zlev(ig,l)-zlev(ig,l-1)
274        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
275        zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
276        zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
277        zdw2=afact*zbuoy(ig,l)/fact_epsilon
278        zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
279        lm=Max(1,l-2)
280        w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
281       endif
282    enddo
283
284
285!-------------------------------------------------
286!calcul des taux d'entrainement et de detrainement
287!-------------------------------------------------
288
289     do ig=1,ngrid
290        if (active(ig)) then
291
292!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
293          zw2m=w_est(ig,l+1)
294          zdz=zlev(ig,l+1)-zlev(ig,l)
295          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
296          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
297          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
298
299!=========================================================================
300! 4. Calcul de l'entrainement et du detrainement
301!=========================================================================
302
303          detr_star(ig,l)=f_star(ig,l)*zdz             &
304       *( mix0 * 0.1 / (zalpha+0.001)               &
305       + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
306       + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
307
308          if ( iflag_thermals_ed == 20 ) then
309             entr_star(ig,l)=f_star(ig,l)*zdz* (         &
310            mix0 * 0.1 / (zalpha+0.001)               &
311          + zbetalpha*MAX(entr_min,                   &
312          afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
313          else
314             entr_star(ig,l)=f_star(ig,l)*zdz* (         &
315            mix0 * 0.1 / (zalpha+0.001)               &
316          + zbetalpha*MAX(entr_min,                   &
317          afact*zbuoy(ig,l)/zw2m - fact_epsilon))
318          endif
319         
320! En dessous de lalim, on prend le max de alim_star et entr_star pour
321! alim_star et 0 sinon
322        if (l<lalim(ig)) then
323          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
324          entr_star(ig,l)=0.
325        endif
326        f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
327                -detr_star(ig,l)
328
329      endif
330   enddo
331
332
333!============================================================================
334! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
335!===========================================================================
336
337   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
338   do ig=1,ngrid
339       if (activetmp(ig)) then
340           Zsat=.false.
341           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
342              (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
343              /(f_star(ig,l+1)+detr_star(ig,l))
344           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
345              (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
346              /(f_star(ig,l+1)+detr_star(ig,l))
347
348        endif
349    enddo
350
351   ztemp(:)=zpspsk(:,l)*ztla(:,l)
352   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
353   do ig=1,ngrid
354      if (activetmp(ig)) then
355! on ecrit de maniere conservative (sat ou non)
356!          T = Tl +Lv/Cp ql
357           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
358           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
359           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
360!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
361           zha(ig,l) = ztva(ig,l)
362           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
363                -zqla(ig,l))-zqla(ig,l))
364           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
365           zdz=zlev(ig,l+1)-zlev(ig,l)
366           zdzbis=zlev(ig,l)-zlev(ig,l-1)
367           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
368           zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
369           zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
370           zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
371           zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
372           zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
373      endif
374   enddo
375
376   if (prt_level>=20) print*,'coucou calcul detr 460: ig, l',ig, l
377
378!===========================================================================
379! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
380!===========================================================================
381
382   nbpb=0
383   do ig=1,ngrid
384            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1)<1.e-10) then
385!               stop'On tombe sur le cas particulier de thermcell_dry'
386!               print*,'On tombe sur le cas particulier de thermcell_plume'
387                nbpb=nbpb+1
388                zw2(ig,l+1)=0.
389                linter(ig)=l+1
390            endif
391
392        if (zw2(ig,l+1)<0.) then
393           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
394                 -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
395           zw2(ig,l+1)=0.
396!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
397        elseif (f_star(ig,l+1)<0.) then
398           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
399                 -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
400           zw2(ig,l+1)=0.
401!fin CR:04/05/12
402        endif
403
404           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
405
406        if (wa_moy(ig,l+1)>wmaxa(ig)) then
407!   lmix est le niveau de la couche ou w (wa_moy) est maximum
408!on rajoute le calcul de lmix_bis
409            if (zqla(ig,l)<1.e-10) then
410               lmix_bis(ig)=l+1
411            endif
412            lmix(ig)=l+1
413            wmaxa(ig)=wa_moy(ig,l+1)
414        endif
415   enddo
416
417   if (nbpb>0) then
418   print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
419   endif
420
421!=========================================================================
422! FIN DE LA BOUCLE VERTICALE
423      enddo
424!=========================================================================
425
426!on recalcule alim_star_tot
427       do ig=1,ngrid
428          alim_star_tot(ig)=0.
429       enddo
430       do ig=1,ngrid
431          do l=1,lalim(ig)-1
432          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
433          enddo
434       enddo
435       
436
437        if (prt_level>=20) print*,'coucou calcul detr 470: ig, l', ig, l
438
439#undef wrgrads_thermcell
440#ifdef wrgrads_thermcell
441         call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta      ','esta      ')
442         call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta      ','dsta      ')
443         call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy      ','buoy      ')
444         call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt      ','dqt      ')
445         call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est     ','w_est     ')
446         call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2     ','w_es2     ')
447         call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A      ','zw2A      ')
448#endif
449
450
451 RETURN
452     end
453END MODULE lmdz_thermcell_plume
Note: See TracBrowser for help on using the repository browser.