source: trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90 @ 2127

Last change on this file since 2127 was 2127, checked in by aboissinot, 6 years ago

Fix in convadj.F
New version of the thermal plume model (new entrainment and detrainment formulae, alimentation is removed)

File size: 18.9 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_main(ngrid,nlay,nq,ptimestep,firstcall,                  &
5                          pplay,pplev,pphi,zpopsk,                            &
6                          pu,pv,pt,pq,                                        &
7                          pduadj,pdvadj,pdtadj,pdqadj,                        &
8                          f0,fm0,entr0,detr0,zw2,fraca,                       &
9                          zqta,zqla,ztv,ztva,zhla,zhl,zqsa,                   &
10                          lmin,lmix,lalim,lmax)
11     
12     
13!===============================================================================
14!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
15!   Version du 09.02.07
16!   Calcul du transport vertical dans la couche limite en presence
17!   de "thermiques" explicitement representes avec processus nuageux
18!
19!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
20!
21!   le thermique est suppose homogene et dissipe par melange avec
22!   son environnement. la longueur l_mix controle l'efficacite du
23!   melange
24!
25!   Le calcul du transport des differentes especes se fait en prenant
26!   en compte:
27!     1. un flux de masse montant
28!     2. un flux de masse descendant
29!     3. un entrainement
30!     4. un detrainement
31!
32! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
33!    Introduction of an implicit computation of vertical advection in
34!    the environment of thermal plumes in thermcell_dq
35!    impl = 0 : explicit ; impl = 1 : implicit ; impl =-1 : old version
36!    controled by iflag_thermals =
37!       15, 16 run with impl=-1 : numerical convergence with NPv3
38!       17, 18 run with impl=1  : more stable
39!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
40!
41! Major changes 2018-19 (AB alexandre.boissinot@lmd.jussieu.fr)
42!    New detr and entre formulae (no longer alimentation)
43!    lmin can be greater than 1
44!    Mix every tracer (EN COURS)
45!    Old version of thermcell_dq is removed
46!    Alternative version thermcell_dv2 is removed
47!
48!===============================================================================
49     
50      USE thermcell_mod
51      USE tracer_h, ONLY: igcm_h2o_vap
52      USE print_control_mod, ONLY: lunout, prt_level
53     
54      IMPLICIT NONE
55     
56     
57!===============================================================================
58! Declaration
59!===============================================================================
60     
61!     Inputs:
62!     -------
63     
64      INTEGER ngrid, nlay, nq
65     
66      REAL ptimestep
67      REAL pplay(ngrid,nlay)                    ! Layer pressure
68      REAL pplev(ngrid,nlay+1)                  ! Level pressure
69      REAL pphi(ngrid,nlay)                     ! Geopotential
70     
71      REAL pu(ngrid,nlay)                       ! Zonal wind
72      REAL pv(ngrid,nlay)                       ! Meridional wind
73      REAL pt(ngrid,nlay)                       ! Temperature
74      REAL pq(ngrid,nlay,nq)                    ! Tracers mass mixing ratio
75     
76      LOGICAL firstcall
77     
78!     Outputs:
79!     --------
80     
81      REAL pduadj(ngrid,nlay)                   ! u convective variations
82      REAL pdvadj(ngrid,nlay)                   ! v convective variations
83      REAL pdtadj(ngrid,nlay)                   ! t convective variations
84      REAL pdqadj(ngrid,nlay,nq)                ! q convective variations
85     
86      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
87      REAL fm0(ngrid,nlay+1)                    ! mass flux      (after possible time relaxation)
88      REAL entr0(ngrid,nlay)                    ! entrainment    (after possible time relaxation)
89      REAL detr0(ngrid,nlay)                    ! detrainment    (after possible time relaxation)
90     
91!     Local:
92!     ------
93     
94      INTEGER ig, k, l, iq
95      INTEGER lmax(ngrid)                       ! Highest layer reached by the plume
96      INTEGER lmix(ngrid)                       ! Layer in which plume vertical speed is maximal
97      INTEGER lmin(ngrid)                       ! First unstable layer
98     
99      REAL zlay(ngrid,nlay)                     ! Layers altitudes
100      REAL zlev(ngrid,nlay+1)                   ! Levels altitudes
101      REAL rho(ngrid,nlay)                      ! Layers densities
102      REAL rhobarz(ngrid,nlay)                  ! Levels densities
103      REAL masse(ngrid,nlay)                    ! Layers masses
104      REAL zpopsk(ngrid,nlay)                   ! Exner function
105     
106      REAL zu(ngrid,nlay)                       ! u    environment
107      REAL zv(ngrid,nlay)                       ! v    environment
108      REAL zt(ngrid,nlay)                       ! TR   environment
109      REAL zqt(ngrid,nlay)                      ! qt   environment
110      REAL zql(ngrid,nlay)                      ! ql   environment
111      REAL zhl(ngrid,nlay)                      ! TP   environment
112      REAL ztv(ngrid,nlay)                      ! TRPV environment
113      REAL zqs(ngrid,nlay)                      ! qsat environment
114     
115      REAL zua(ngrid,nlay)                      ! u    plume
116      REAL zva(ngrid,nlay)                      ! v    plume
117      REAL zta(ngrid,nlay)                      ! TR   plume
118      REAL zqla(ngrid,nlay)                     ! qv   plume
119      REAL zqta(ngrid,nlay)                     ! qt   plume
120      REAL zhla(ngrid,nlay)                     ! TP   plume
121      REAL ztva(ngrid,nlay)                     ! TRPV plume
122      REAL zqsa(ngrid,nlay)                     ! qsat plume
123     
124      REAL zqa(ngrid,nlay,nq)                   ! q    plume (ql=0, qv=qt)
125     
126      REAL linter(ngrid)                        ! Level (continuous) of maximal vertical speed
127      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
128      REAL zmax(ngrid)                          ! Maximal altitude reached by the plume
129      REAL wmax(ngrid)                          ! Maximal vertical speed
130      REAL zw2(ngrid,nlay+1)                    ! Plume vertical speed
131     
132      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
133     
134      REAL f_star(ngrid,nlay+1)                 ! Normalized mass flux
135      REAL entr_star(ngrid,nlay)                ! Normalized entrainment
136      REAL detr_star(ngrid,nlay)                ! Normalized detrainment
137     
138      REAL f(ngrid)                             ! Mass flux norm
139      REAL fm(ngrid,nlay+1)                     ! Mass flux
140      REAL entr(ngrid,nlay)                     ! Entrainment
141      REAL detr(ngrid,nlay)                     ! Detrainment
142     
143      REAL lambda                               ! Time relaxation coefficent
144     
145      REAL zdthladj(ngrid,nlay)                 ! Potential temperature variations
146      REAL dummy(ngrid,nlay)                    ! Dummy argument for thermcell_dq()
147     
148      CHARACTER (len=20) :: modname='thermcell_main'
149      CHARACTER (len=80) :: abort_message
150     
151! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
152      INTEGER lalim(ngrid)                      ! Highest alimentation level
153      REAL alim_star(ngrid,nlay)                ! Normalized alimentation
154      REAL alim_star_tot(ngrid)                 ! Integrated alimentation
155      REAL alim_star_clos(ngrid,nlay)           ! Closure alimentation
156     
157!===============================================================================
158! Initialization
159!===============================================================================
160     
161      IF (firstcall) THEN
162         fm0(:,:) = 0.
163         entr0(:,:) = 0.
164         detr0(:,:) = 0.
165      ENDIF
166     
167      f_star(:,:) = 0.
168      entr_star(:,:) = 0.
169      detr_star(:,:) = 0.
170     
171      f(:) = 0.
172     
173      fm(:,:) = 0.
174      entr(:,:) = 0.
175      detr(:,:) = 0.
176     
177      lmax(:) = 1
178      lmix(:) = 1
179      lmin(:) = 1
180     
181      pduadj(:,:) = 0.0
182      pdvadj(:,:) = 0.0
183      pdtadj(:,:) = 0.0
184      pdqadj(:,:,:) = 0.0
185     
186! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
187      alim_star(:,:) = 0.
188      alim_star_tot(:) = 0.
189      alim_star_clos(:,:) = 0.
190      lalim(:) = 1
191     
192! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model!
193      DO ig=1,ngrid
194         f0(ig) = max(f0(ig), 1.e-2)
195      ENDDO
196     
197      IF (prt_level.ge.20) then
198         DO ig=1,ngrid
199            print *, 'ig,f0', ig, f0(ig)
200         ENDDO
201      ENDIF
202     
203!===============================================================================
204! Environment settings
205!===============================================================================
206     
207!-------------------------------------------------------------------------------
208! Calcul de T,q,ql a partir de Tl et qt dans l environnement
209!-------------------------------------------------------------------------------
210     
211      CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
212      &                  zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
213     
214!-------------------------------------------------------------------------------
215! Levels and layers altitudes
216!-------------------------------------------------------------------------------
217     
218      DO l=2,nlay
219         zlev(:,l) = 0.5 * (pphi(:,l) + pphi(:,l-1)) / RG
220      ENDDO
221     
222      zlev(:,1) = 0.
223      zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
224     
225      DO l=1,nlay
226         zlay(:,l) = pphi(:,l)/RG
227      ENDDO
228     
229!-------------------------------------------------------------------------------
230! Levels and layers densities
231!-------------------------------------------------------------------------------
232     
233      rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
234     
235      IF (prt_level.ge.10) THEN
236         write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)'
237      ENDIF
238     
239      rhobarz(:,1) = rho(:,1)
240     
241      DO l=2,nlay
242         rhobarz(:,l) = 0.5 * (rho(:,l) + rho(:,l-1))
243      ENDDO
244     
245!-------------------------------------------------------------------------------
246! Layers masses
247!-------------------------------------------------------------------------------
248     
249      DO l=1,nlay
250         masse(:,l) = (pplev(:,l) - pplev(:,l+1)) / RG
251      ENDDO
252     
253!===============================================================================
254! Explicative schemes
255!===============================================================================
256
257!-------------------------------------------------------------------------------
258! Thermal plume variables
259!-------------------------------------------------------------------------------
260     
261!           top of the model     
262!     ===========================
263!                               
264!     ---------------------------
265!                                        _
266!     ----- F_lmax+1=0 ------zmax         \
267!     lmax                                 |
268!     ------F_lmax>0-------------          |
269!                                          |
270!     ---------------------------          |
271!                                          |
272!     ---------------------------          |
273!                                          |
274!     ------------------wmax,zmix          |
275!     lmix                                 |
276!     ---------------------------          |
277!                                          |
278!     ---------------------------          |
279!                                          | E, D
280!     ---------------------------          |
281!                                          |
282!     --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca
283!         zt, zu, zv, zo, rho              |
284!     ---------------------------          |
285!                                          |
286!     ---------------------------          |
287!                                          |
288!     ---------------------------          |
289!                                          |
290!     ------F_lmin+1>0-----------          |
291!     lmin                                 |
292!     ----- F_lmin=0 ------------        _/
293!                               
294!     ---------------------------
295!                               
296!     ===========================
297!         bottom of the model   
298     
299!-------------------------------------------------------------------------------
300! Zoom on layers k and k-1
301!-------------------------------------------------------------------------------
302     
303!     |     /|\                  |                          |
304!     |----  |  F_k+1 -----------|--------------------------|   level k+1
305!     |      |  w_k+1         |                             |
306!     |                     --|--> D_k                      |
307!     |                       |                             |   layer k
308!     |                    <--|--  E_k                      |
309!     |     /|\               |                             |
310!     |----  |  F_k ----------|-----------------------------|   level k
311!     |      |  w_k        |                                |
312!     |                  --|--> D_k-1                       |
313!     |                    |                                |   layer k-1
314!     |                 <--|--  E_k-1                       |
315!     |     /|\            |                                |
316!     |----  |  F_k-1 -----|--------------------------------|   level k-1
317!            |  w_k-1                                       
318!     0                  fraca                              1
319!      \__________________/ \______________________________/
320!          plume (fraca)          environment (1-fraca)   
321     
322!===============================================================================
323! Thermal plumes computation
324!===============================================================================
325     
326!-------------------------------------------------------------------------------
327! Thermal plumes speeds, fluxes, tracers and temperatures
328!-------------------------------------------------------------------------------
329     
330      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
331      &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
332      &                    detr_star,entr_star,f_star,                        &
333      &                    ztva,zhla,zqla,zqta,zta,                           &
334      &                    zw2,zqsa,lmix,lmin)
335     
336!-------------------------------------------------------------------------------
337! Thermal plumes characteristics: zmax, zmix, wmax
338!-------------------------------------------------------------------------------
339     
340! AB: Careful, zw2 became its square root in thermcell_height!
341      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
342      &                     zlev,lmax,zmax,zmix,wmax,f_star)
343     
344!===============================================================================
345! Closure and mass fluxes computation
346!===============================================================================
347     
348!-------------------------------------------------------------------------------
349! Closure
350!-------------------------------------------------------------------------------
351     
352      CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
353      &                      lmax,entr_star,zmax,wmax,f)
354     
355      IF (tau_thermals>1.) THEN
356         lambda = exp(-ptimestep/tau_thermals)
357         f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
358      ELSE
359         f0(:) = f(:)
360      ENDIF
361     
362!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
363! Test valable seulement en 1D mais pas genant
364!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
365      IF (.not. (f0(1).ge.0.) ) THEN
366         abort_message = '.not. (f0(1).ge.0.)'
367         print *, 'f0 =', f0(1)
368         CALL abort_physic(modname,abort_message,1)
369      ENDIF
370     
371!-------------------------------------------------------------------------------
372! Mass fluxes
373!-------------------------------------------------------------------------------
374     
375      CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
376! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
377!     That is not already done for thermcell_flux.
378      &                   lalim,lmin,lmax,entr_star,detr_star,                &
379      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
380     
381!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
382! On ne prend pas directement les profils issus des calculs precedents mais on
383! s'autorise genereusement une relaxation vers ceci avec une constante de temps
384! tau_thermals (typiquement 1800s sur Terre).
385!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
386     
387      IF (tau_thermals>1.) THEN
388         lambda = exp(-ptimestep/tau_thermals)
389         fm0   = (1.-lambda) * fm   + lambda * fm0
390         entr0 = (1.-lambda) * entr + lambda * entr0
391         detr0 = (1.-lambda) * detr + lambda * detr0
392      ELSE
393         fm0(:,:)   = fm(:,:)
394         entr0(:,:) = entr(:,:)
395         detr0(:,:) = detr(:,:)
396      ENDIF
397     
398!-------------------------------------------------------------------------------
399! Updraft fraction
400!-------------------------------------------------------------------------------
401     
402      DO ig=1,ngrid
403         fraca(ig,1) = 0.
404         fraca(ig,nlay+1) = 0.
405      ENDDO
406     
407      DO l=2,nlay
408         DO ig=1,ngrid
409            IF (zw2(ig,l).gt.0.) THEN
410               fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
411            ELSE
412               fraca(ig,l) = 0.
413            ENDIF
414         ENDDO
415      ENDDO
416     
417!===============================================================================
418! Transport vertical
419!===============================================================================
420     
421!-------------------------------------------------------------------------------
422! Calcul du transport vertical de la temperature potentielle
423!-------------------------------------------------------------------------------
424     
425      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
426      &                 zhl,zdthladj,dummy,lmin)
427     
428      DO l=1,nlay
429         DO ig=1,ngrid
430            pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l) 
431         ENDDO
432      ENDDO
433     
434!-------------------------------------------------------------------------------
435! Calcul du transport vertical des traceurs
436!-------------------------------------------------------------------------------
437     
438      DO iq=1,nq
439         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,        &
440         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
441      ENDDO
442     
443!-------------------------------------------------------------------------------
444! Calcul du transport vertical du moment horizontal
445!-------------------------------------------------------------------------------
446     
447      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
448      &                 zu,pduadj,zua,lmin)
449     
450      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
451      &                 zv,pdvadj,zva,lmin)
452     
453     
454RETURN
455END
Note: See TracBrowser for help on using the repository browser.