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

Last change on this file since 2176 was 2144, checked in by aboissinot, 5 years ago

Thermal plume model parameters are now set in callphys.def and (re)add flag to choose between thermcell_dv2 or thermcell_dq to mix horizontal momentum.

File size: 18.4 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,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 print_control_mod, ONLY: prt_level
52     
53      IMPLICIT NONE
54     
55     
56!===============================================================================
57! Declaration
58!===============================================================================
59     
60!     Inputs:
61!     -------
62     
63      INTEGER ngrid, nlay, nq
64     
65      REAL ptimestep
66      REAL pplay(ngrid,nlay)                    ! Layer pressure
67      REAL pplev(ngrid,nlay+1)                  ! Level pressure
68      REAL pphi(ngrid,nlay)                     ! Geopotential
69     
70      REAL pu(ngrid,nlay)                       ! Zonal wind
71      REAL pv(ngrid,nlay)                       ! Meridional wind
72      REAL pt(ngrid,nlay)                       ! Temperature
73      REAL pq(ngrid,nlay,nq)                    ! Tracers mass mixing ratio
74     
75      LOGICAL firstcall
76     
77!     Outputs:
78!     --------
79     
80      REAL pduadj(ngrid,nlay)                   ! u convective variations
81      REAL pdvadj(ngrid,nlay)                   ! v convective variations
82      REAL pdtadj(ngrid,nlay)                   ! t convective variations
83      REAL pdqadj(ngrid,nlay,nq)                ! q convective variations
84     
85      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
86      REAL fm0(ngrid,nlay+1)                    ! mass flux      (after possible time relaxation)
87      REAL entr0(ngrid,nlay)                    ! entrainment    (after possible time relaxation)
88      REAL detr0(ngrid,nlay)                    ! detrainment    (after possible time relaxation)
89     
90!     Local:
91!     ------
92     
93      INTEGER ig, k, l, iq
94      INTEGER lmax(ngrid)                       ! Highest layer reached by the plume
95      INTEGER lmix(ngrid)                       ! Layer in which plume vertical speed is maximal
96      INTEGER lmin(ngrid)                       ! First unstable layer
97     
98      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
99      REAL zmax(ngrid)                          ! Maximal altitudes where plumes are active
100      REAL zmin(ngrid)                          ! Minimal altitudes where plumes are active
101     
102      REAL zlay(ngrid,nlay)                     ! Layers altitudes
103      REAL zlev(ngrid,nlay+1)                   ! Levels altitudes
104      REAL rho(ngrid,nlay)                      ! Layers densities
105      REAL rhobarz(ngrid,nlay)                  ! Levels densities
106      REAL masse(ngrid,nlay)                    ! Layers masses
107      REAL zpopsk(ngrid,nlay)                   ! Exner function
108     
109      REAL zu(ngrid,nlay)                       ! u    environment
110      REAL zv(ngrid,nlay)                       ! v    environment
111      REAL zt(ngrid,nlay)                       ! TR   environment
112      REAL zqt(ngrid,nlay)                      ! qt   environment
113      REAL zql(ngrid,nlay)                      ! ql   environment
114      REAL zhl(ngrid,nlay)                      ! TP   environment
115      REAL ztv(ngrid,nlay)                      ! TRPV environment
116      REAL zqs(ngrid,nlay)                      ! qsat environment
117     
118      REAL zua(ngrid,nlay)                      ! u    plume
119      REAL zva(ngrid,nlay)                      ! v    plume
120      REAL zta(ngrid,nlay)                      ! TR   plume
121      REAL zqla(ngrid,nlay)                     ! qv   plume
122      REAL zqta(ngrid,nlay)                     ! qt   plume
123      REAL zhla(ngrid,nlay)                     ! TP   plume
124      REAL ztva(ngrid,nlay)                     ! TRPV plume
125      REAL zqsa(ngrid,nlay)                     ! qsat plume
126     
127      REAL zqa(ngrid,nlay,nq)                   ! q    plume (ql=0, qv=qt)
128     
129      REAL f_star(ngrid,nlay+1)                 ! Normalized mass flux
130      REAL entr_star(ngrid,nlay)                ! Normalized entrainment (E* = e* dz)
131      REAL detr_star(ngrid,nlay)                ! Normalized detrainment (D* = d* dz)
132     
133      REAL fm(ngrid,nlay+1)                     ! Mass flux
134      REAL entr(ngrid,nlay)                     ! Entrainment (E = e dz)
135      REAL detr(ngrid,nlay)                     ! Detrainment (D = d dz)
136     
137      REAL f(ngrid)                             ! Mass flux norm
138      REAL lambda                               ! Time relaxation coefficent
139      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
140      REAL linter(ngrid)                        ! Level (continuous) of maximal vertical speed
141      REAL wmax(ngrid)                          ! Maximal vertical speed
142      REAL zw2(ngrid,nlay+1)                    ! Plume vertical speed
143      REAL zdthladj(ngrid,nlay)                 ! Potential temperature variations
144      REAL dummy(ngrid,nlay)                    ! Dummy argument for thermcell_dq()
145     
146!===============================================================================
147! Initialization
148!===============================================================================
149     
150      IF (firstcall) THEN
151         fm0(:,:) = 0.
152         entr0(:,:) = 0.
153         detr0(:,:) = 0.
154      ENDIF
155     
156      f_star(:,:) = 0.
157      entr_star(:,:) = 0.
158      detr_star(:,:) = 0.
159     
160      f(:) = 0.
161     
162      fm(:,:) = 0.
163      entr(:,:) = 0.
164      detr(:,:) = 0.
165     
166      lmax(:) = 1
167      lmix(:) = 1
168      lmin(:) = 1
169     
170      pduadj(:,:) = 0.0
171      pdvadj(:,:) = 0.0
172      pdtadj(:,:) = 0.0
173      pdqadj(:,:,:) = 0.0
174     
175      DO ig=1,ngrid
176! AB: Careful: Hard-coded value from Earth version!
177!         f0(ig) = max(f0(ig), 1.e-2)
178! AB: No pescribed minimal value for f0
179         f0(ig) = max(f0(ig), 0.)
180      ENDDO
181     
182!===============================================================================
183! Environment settings
184!===============================================================================
185     
186!-------------------------------------------------------------------------------
187! Calcul de T,q,ql a partir de Tl et qt dans l environnement
188!-------------------------------------------------------------------------------
189     
190      CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
191      &                  zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
192     
193!-------------------------------------------------------------------------------
194! Levels and layers altitudes
195!-------------------------------------------------------------------------------
196     
197      DO l=2,nlay
198         zlev(:,l) = 0.5 * (pphi(:,l) + pphi(:,l-1)) / RG
199      ENDDO
200     
201      zlev(:,1) = 0.
202      zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
203     
204      DO l=1,nlay
205         zlay(:,l) = pphi(:,l) / RG
206      ENDDO
207     
208!-------------------------------------------------------------------------------
209! Levels and layers densities
210!-------------------------------------------------------------------------------
211     
212      rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
213     
214      IF (prt_level.ge.10) THEN
215         print *, 'WARNING: density in the first layer is equal to density at the first level!'
216         print *, 'rhobarz(:,1)', rhobarz(:,1)
217         print *, 'rho(:,1)    ', rho(:,1)
218      ENDIF
219     
220      rhobarz(:,1) = rho(:,1)
221     
222      DO l=2,nlay
223         rhobarz(:,l) = 0.5 * (rho(:,l) + rho(:,l-1))
224      ENDDO
225     
226!-------------------------------------------------------------------------------
227! Layers masses
228!-------------------------------------------------------------------------------
229     
230      DO l=1,nlay
231         masse(:,l) = (pplev(:,l) - pplev(:,l+1)) / RG
232      ENDDO
233     
234!===============================================================================
235! Explicative schemes
236!===============================================================================
237
238!-------------------------------------------------------------------------------
239! Thermal plume variables
240!-------------------------------------------------------------------------------
241     
242!           top of the model     
243!     ===========================
244!                               
245!     ---------------------------
246!                                        _
247!     ----- F_lmax+1=0 ------zmax         \
248!     lmax                                 |
249!     ------F_lmax>0-------------          |
250!                                          |
251!     ---------------------------          |
252!                                          |
253!     ---------------------------          |
254!                                          |
255!     ------------------wmax,zmix          |
256!     lmix                                 |
257!     ---------------------------          |
258!                                          |
259!     ---------------------------          |
260!                                          | E, D
261!     ---------------------------          |
262!                                          |
263!     --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca
264!         zt, zu, zv, zo, rho              |
265!     ---------------------------          |
266!                                          |
267!     ---------------------------          |
268!                                          |
269!     ---------------------------          |
270!                                          |
271!     ------F_lmin+1>0-----------          |
272!     lmin                                 |
273!     ----- F_lmin=0 ------------        _/
274!                               
275!     ---------------------------
276!                               
277!     ===========================
278!         bottom of the model   
279     
280!-------------------------------------------------------------------------------
281! Zoom on layers k and k-1
282!-------------------------------------------------------------------------------
283     
284!     |     /|\                  |                          |
285!     |----  |  F_k+1 -----------|--------------------------|   level k+1
286!     |      |  w_k+1         |                             |
287!     |                     --|--> D_k                      |
288!     |                       |                             |   layer k
289!     |                    <--|--  E_k                      |
290!     |     /|\               |                             |
291!     |----  |  F_k ----------|-----------------------------|   level k
292!     |      |  w_k        |                                |
293!     |                  --|--> D_k-1                       |
294!     |                    |                                |   layer k-1
295!     |                 <--|--  E_k-1                       |
296!     |     /|\            |                                |
297!     |----  |  F_k-1 -----|--------------------------------|   level k-1
298!            |  w_k-1                                       
299!     0                  fraca                              1
300!      \__________________/ \______________________________/
301!          plume (fraca)          environment (1-fraca)   
302     
303!===============================================================================
304! Thermal plumes computation
305!===============================================================================
306     
307!-------------------------------------------------------------------------------
308! Thermal plumes speeds, fluxes, tracers and temperatures
309!-------------------------------------------------------------------------------
310     
311      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
312      &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
313      &                    detr_star,entr_star,f_star,                        &
314      &                    ztva,zhla,zqla,zqta,zta,zqsa,                      &
315      &                    zw2,lmix,lmin)
316     
317!-------------------------------------------------------------------------------
318! Thermal plumes characteristics: zmax, zmix, wmax
319!-------------------------------------------------------------------------------
320     
321! AB: Careful, zw2 became its square root in thermcell_height!
322      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
323      &                     zlev,zmin,zmix,zmax,wmax,f_star)
324     
325!===============================================================================
326! Closure and mass fluxes computation
327!===============================================================================
328     
329!-------------------------------------------------------------------------------
330! Closure
331!-------------------------------------------------------------------------------
332     
333      CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
334      &                      lmax,entr_star,zmin,zmax,wmax,f)
335     
336      IF (tau_thermals>1.) THEN
337         lambda = exp(-ptimestep/tau_thermals)
338         f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
339      ELSE
340         f0(:) = f(:)
341      ENDIF
342     
343!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
344! Test valable seulement en 1D mais pas genant
345!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
346      IF (.not. (f0(1).ge.0.) ) THEN
347         print *, 'ERROR: mass flux norm is not positive!'
348         print *, 'f0 =', f0(1)
349         CALL abort
350      ENDIF
351     
352!-------------------------------------------------------------------------------
353! Mass fluxes
354!-------------------------------------------------------------------------------
355     
356      CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
357      &                   lmin,lmax,entr_star,detr_star,                      &
358      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
359     
360!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
361! On ne prend pas directement les profils issus des calculs precedents mais on
362! s'autorise genereusement une relaxation vers ceci avec une constante de temps
363! tau_thermals (typiquement 1800s sur Terre).
364!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
365     
366      IF (tau_thermals>1.) THEN
367         lambda = exp(-ptimestep/tau_thermals)
368         fm0   = (1.-lambda) * fm   + lambda * fm0
369         entr0 = (1.-lambda) * entr + lambda * entr0
370         detr0 = (1.-lambda) * detr + lambda * detr0
371      ELSE
372         fm0(:,:)   = fm(:,:)
373         entr0(:,:) = entr(:,:)
374         detr0(:,:) = detr(:,:)
375      ENDIF
376     
377!-------------------------------------------------------------------------------
378! Updraft fraction
379!-------------------------------------------------------------------------------
380     
381      DO ig=1,ngrid
382         fraca(ig,1) = 0.
383         fraca(ig,nlay+1) = 0.
384      ENDDO
385     
386      DO l=2,nlay
387         DO ig=1,ngrid
388            IF (zw2(ig,l) > 0.) THEN
389               fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
390            ELSE
391               fraca(ig,l) = 0.
392            ENDIF
393         ENDDO
394      ENDDO
395     
396!===============================================================================
397! Transport vertical
398!===============================================================================
399     
400!-------------------------------------------------------------------------------
401! Calcul du transport vertical de la temperature potentielle
402!-------------------------------------------------------------------------------
403     
404      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
405      &                 zhl,zdthladj,dummy,lmin)
406     
407      DO l=1,nlay
408         DO ig=1,ngrid
409            pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l) 
410         ENDDO
411      ENDDO
412     
413!-------------------------------------------------------------------------------
414! Calcul du transport vertical des traceurs
415!-------------------------------------------------------------------------------
416     
417      DO iq=1,nq
418         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,        &
419         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
420      ENDDO
421     
422!-------------------------------------------------------------------------------
423! Calcul du transport vertical du moment horizontal
424!-------------------------------------------------------------------------------
425     
426      IF (dvimpl) THEN
427         CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,fraca,    &
428         &                  zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva)
429      ELSE
430         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
431         &                 zu,pduadj,zua,lmin)
432         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
433         &                 zv,pdvadj,zva,lmin)
434      ENDIF
435     
436     
437RETURN
438END
Note: See TracBrowser for help on using the repository browser.