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

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

Alimentation (more precisely lalim) is removed from thermcell_flux as well as corrections which need it.
Consequently, iflag_thermals_optflux flag is removed from thermcell_mod.

File size: 18.2 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 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)                          ! Plume height
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!===============================================================================
152! Initialization
153!===============================================================================
154     
155      IF (firstcall) THEN
156         fm0(:,:) = 0.
157         entr0(:,:) = 0.
158         detr0(:,:) = 0.
159      ENDIF
160     
161      f_star(:,:) = 0.
162      entr_star(:,:) = 0.
163      detr_star(:,:) = 0.
164     
165      f(:) = 0.
166     
167      fm(:,:) = 0.
168      entr(:,:) = 0.
169      detr(:,:) = 0.
170     
171      lmax(:) = 1
172      lmix(:) = 1
173      lmin(:) = 1
174     
175      pduadj(:,:) = 0.0
176      pdvadj(:,:) = 0.0
177      pdtadj(:,:) = 0.0
178      pdqadj(:,:,:) = 0.0
179     
180! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model!
181      DO ig=1,ngrid
182         f0(ig) = max(f0(ig), 1.e-2)
183      ENDDO
184     
185      IF (prt_level.ge.20) then
186         DO ig=1,ngrid
187            print *, 'ig,f0', ig, f0(ig)
188         ENDDO
189      ENDIF
190     
191!===============================================================================
192! Environment settings
193!===============================================================================
194     
195!-------------------------------------------------------------------------------
196! Calcul de T,q,ql a partir de Tl et qt dans l environnement
197!-------------------------------------------------------------------------------
198     
199      CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
200      &                  zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
201     
202!-------------------------------------------------------------------------------
203! Levels and layers altitudes
204!-------------------------------------------------------------------------------
205     
206      DO l=2,nlay
207         zlev(:,l) = 0.5 * (pphi(:,l) + pphi(:,l-1)) / RG
208      ENDDO
209     
210      zlev(:,1) = 0.
211      zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
212     
213      DO l=1,nlay
214         zlay(:,l) = pphi(:,l)/RG
215      ENDDO
216     
217!-------------------------------------------------------------------------------
218! Levels and layers densities
219!-------------------------------------------------------------------------------
220     
221      rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
222     
223      IF (prt_level.ge.10) THEN
224         write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)'
225      ENDIF
226     
227      rhobarz(:,1) = rho(:,1)
228     
229      DO l=2,nlay
230         rhobarz(:,l) = 0.5 * (rho(:,l) + rho(:,l-1))
231      ENDDO
232     
233!-------------------------------------------------------------------------------
234! Layers masses
235!-------------------------------------------------------------------------------
236     
237      DO l=1,nlay
238         masse(:,l) = (pplev(:,l) - pplev(:,l+1)) / RG
239      ENDDO
240     
241!===============================================================================
242! Explicative schemes
243!===============================================================================
244
245!-------------------------------------------------------------------------------
246! Thermal plume variables
247!-------------------------------------------------------------------------------
248     
249!           top of the model     
250!     ===========================
251!                               
252!     ---------------------------
253!                                        _
254!     ----- F_lmax+1=0 ------zmax         \
255!     lmax                                 |
256!     ------F_lmax>0-------------          |
257!                                          |
258!     ---------------------------          |
259!                                          |
260!     ---------------------------          |
261!                                          |
262!     ------------------wmax,zmix          |
263!     lmix                                 |
264!     ---------------------------          |
265!                                          |
266!     ---------------------------          |
267!                                          | E, D
268!     ---------------------------          |
269!                                          |
270!     --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca
271!         zt, zu, zv, zo, rho              |
272!     ---------------------------          |
273!                                          |
274!     ---------------------------          |
275!                                          |
276!     ---------------------------          |
277!                                          |
278!     ------F_lmin+1>0-----------          |
279!     lmin                                 |
280!     ----- F_lmin=0 ------------        _/
281!                               
282!     ---------------------------
283!                               
284!     ===========================
285!         bottom of the model   
286     
287!-------------------------------------------------------------------------------
288! Zoom on layers k and k-1
289!-------------------------------------------------------------------------------
290     
291!     |     /|\                  |                          |
292!     |----  |  F_k+1 -----------|--------------------------|   level k+1
293!     |      |  w_k+1         |                             |
294!     |                     --|--> D_k                      |
295!     |                       |                             |   layer k
296!     |                    <--|--  E_k                      |
297!     |     /|\               |                             |
298!     |----  |  F_k ----------|-----------------------------|   level k
299!     |      |  w_k        |                                |
300!     |                  --|--> D_k-1                       |
301!     |                    |                                |   layer k-1
302!     |                 <--|--  E_k-1                       |
303!     |     /|\            |                                |
304!     |----  |  F_k-1 -----|--------------------------------|   level k-1
305!            |  w_k-1                                       
306!     0                  fraca                              1
307!      \__________________/ \______________________________/
308!          plume (fraca)          environment (1-fraca)   
309     
310!===============================================================================
311! Thermal plumes computation
312!===============================================================================
313     
314!-------------------------------------------------------------------------------
315! Thermal plumes speeds, fluxes, tracers and temperatures
316!-------------------------------------------------------------------------------
317     
318      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
319      &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
320      &                    detr_star,entr_star,f_star,                        &
321      &                    ztva,zhla,zqla,zqta,zta,                           &
322      &                    zw2,zqsa,lmix,lmin)
323           
324!-------------------------------------------------------------------------------
325! Thermal plumes characteristics: zmax, zmix, wmax
326!-------------------------------------------------------------------------------
327     
328! AB: Careful, zw2 became its square root in thermcell_height!
329      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
330      &                     zlev,lmax,zmax,zmix,wmax,f_star)
331     
332!===============================================================================
333! Closure and mass fluxes computation
334!===============================================================================
335     
336!-------------------------------------------------------------------------------
337! Closure
338!-------------------------------------------------------------------------------
339     
340      CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
341      &                      lmax,entr_star,zmax,wmax,f)
342     
343      IF (tau_thermals>1.) THEN
344         lambda = exp(-ptimestep/tau_thermals)
345         f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
346      ELSE
347         f0(:) = f(:)
348      ENDIF
349     
350!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
351! Test valable seulement en 1D mais pas genant
352!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
353      IF (.not. (f0(1).ge.0.) ) THEN
354         abort_message = '.not. (f0(1).ge.0.)'
355         print *, 'f0 =', f0(1)
356         CALL abort_physic(modname,abort_message,1)
357      ENDIF
358     
359!-------------------------------------------------------------------------------
360! Mass fluxes
361!-------------------------------------------------------------------------------
362     
363      CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
364      &                   lmin,lmax,entr_star,detr_star,                      &
365      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
366     
367!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
368! On ne prend pas directement les profils issus des calculs precedents mais on
369! s'autorise genereusement une relaxation vers ceci avec une constante de temps
370! tau_thermals (typiquement 1800s sur Terre).
371!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
372     
373      IF (tau_thermals>1.) THEN
374         lambda = exp(-ptimestep/tau_thermals)
375         fm0   = (1.-lambda) * fm   + lambda * fm0
376         entr0 = (1.-lambda) * entr + lambda * entr0
377         detr0 = (1.-lambda) * detr + lambda * detr0
378      ELSE
379         fm0(:,:)   = fm(:,:)
380         entr0(:,:) = entr(:,:)
381         detr0(:,:) = detr(:,:)
382      ENDIF
383     
384!-------------------------------------------------------------------------------
385! Updraft fraction
386!-------------------------------------------------------------------------------
387     
388      DO ig=1,ngrid
389         fraca(ig,1) = 0.
390         fraca(ig,nlay+1) = 0.
391      ENDDO
392     
393      DO l=2,nlay
394         DO ig=1,ngrid
395            IF (zw2(ig,l).gt.0.) THEN
396               fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
397            ELSE
398               fraca(ig,l) = 0.
399            ENDIF
400         ENDDO
401      ENDDO
402     
403!===============================================================================
404! Transport vertical
405!===============================================================================
406     
407!-------------------------------------------------------------------------------
408! Calcul du transport vertical de la temperature potentielle
409!-------------------------------------------------------------------------------
410     
411      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
412      &                 zhl,zdthladj,dummy,lmin)
413     
414      DO l=1,nlay
415         DO ig=1,ngrid
416            pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l) 
417         ENDDO
418      ENDDO
419     
420!-------------------------------------------------------------------------------
421! Calcul du transport vertical des traceurs
422!-------------------------------------------------------------------------------
423     
424      DO iq=1,nq
425         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,        &
426         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
427      ENDDO
428     
429!-------------------------------------------------------------------------------
430! Calcul du transport vertical du moment horizontal
431!-------------------------------------------------------------------------------
432     
433      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
434      &                 zu,pduadj,zua,lmin)
435     
436      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
437      &                 zv,pdvadj,zva,lmin)
438     
439     
440RETURN
441END
Note: See TracBrowser for help on using the repository browser.