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

Last change on this file since 3496 was 3451, checked in by fhourdin, 6 years ago

Saving old version of thermcell_plume in thermcell_plume_6A.F90
It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
to the 5B and 6A versions used for CMIP5 and CMIP6.
The latest was previously named thermcellV1_plume.
The new thermcell_plume is a clean version (removing obsolete
options) of thermcell_plume_6A.
The 3 versions are controled by
flag_thermals_ed <= 9 thermcell_plume_6A

<= 19 thermcell_plume_5B
else thermcell_plume (default 20 for convergence with 6A)

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