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

Last change on this file was 2480, checked in by aboissinot, 4 years ago

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