source: trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90 @ 2130

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

More satisfiying initialisations in thermcell_env.F90
Add new informations about some variables in thermcell_plume.F90
File size: 15.2 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
5                           zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
6                           detr_star,entr_star,f_star,                        &
7                           ztva,zhla,zqla,zqta,zta,                           &
8                           zw2,zqsa,lmix,lmin)
9     
10     
11!===============================================================================
12!  Purpose: calcule les valeurs de qt, thetal et w dans l ascendance
13! 
14!  Nota Bene
15!     ql   means "non-gaseous water mass mixing ratio" (liquid and solid)
16!     qv   means "vapor mass mixing ratio"
17!     qt   means "total water mass mixing ratio"
18!     TP   means "potential temperature"
19!     TRPV means "virtual potential temperature with latent heat release" 
20!     TPV  means "virtual potential temperature"
21!     TR   means "temperature with latent heat release"
22!===============================================================================
23     
24      USE print_control_mod, ONLY: prt_level
25      USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
26      USE tracer_h, ONLY: igcm_h2o_vap
27      USE thermcell_mod
28     
29      IMPLICIT NONE
30     
31     
32!===============================================================================
33! Declaration
34!===============================================================================
35     
36!     Inputs:
37!     -------
38     
39      INTEGER ngrid, nlay, nq
40     
41      REAL ptimestep
42      REAL rhobarz(ngrid,nlay)                  ! Levels density
43      REAL zlev(ngrid,nlay+1)                   ! Levels altitude
44      REAL pplev(ngrid,nlay+1)                  ! Levels pressure
45      REAL pphi(ngrid,nlay)                     ! Geopotential
46      REAL zpopsk(ngrid,nlay)                   ! Exner function
47     
48      REAL ztv(ngrid,nlay)                      ! TRPV environment
49      REAL zhl(ngrid,nlay)                      ! TP   environment
50      REAL zqt(ngrid,nlay)                      ! qt   environment
51      REAL zql(ngrid,nlay)                      ! ql   environment
52     
53!     Outputs:
54!     --------
55     
56      INTEGER lmin(ngrid)                       ! plume base level (first unstable level)
57      INTEGER lmix(ngrid)                       ! maximum vertical speed level
58     
59      REAL detr_star(ngrid,nlay)                ! normalized detrainment
60      REAL entr_star(ngrid,nlay)                ! normalized entrainment
61      REAL f_star(ngrid,nlay+1)                 ! normalized mass flux
62     
63      REAL ztva(ngrid,nlay)                     ! TRPV plume (after mixing)
64      REAL zhla(ngrid,nlay)                     ! TP   plume (after mixing)
65      REAL zqla(ngrid,nlay)                     ! ql   plume (after mixing)
66      REAL zqta(ngrid,nlay)                     ! qt   plume (after mixing)
67      REAL zqsa(ngrid,nlay)                     ! qsat plume (after mixing)
68      REAL zw2(ngrid,nlay+1)                    ! w2   plume (after mixing)
69     
70!     Local:
71!     ------
72     
73      INTEGER ig, l, k
74     
75      REAL ztva_est(ngrid,nlay)                 ! TRPV plume (before mixing)
76      REAL zqla_est(ngrid,nlay)                 ! ql   plume (before mixing)
77      REAL zta_est(ngrid,nlay)                  ! TR   plume (before mixing)
78      REAL zqsa_est(ngrid)                      ! qsat plume (before mixing)
79      REAL zw2_est(ngrid,nlay+1)                ! w2   plume (before mixing)
80     
81      REAL zta(ngrid,nlay)                      ! TR   plume (after mixing)
82     
83      REAL zbuoy(ngrid,nlay)                    ! Plume buoyancy
84      REAL ztemp(ngrid)                         ! Temperature for saturation vapor pressure computation in plume
85      REAL zdz                                  ! Layers heights
86      REAL ztv2(ngrid,nlay)                     ! ztv + d_temp * Dirac(l=linf)
87     
88      REAL zbetalpha                            !
89      REAL zdw2                                 !
90      REAL zdw2bis                              !
91      REAL zw2fact                              !
92      REAL zw2m                                 ! Average vertical velocity between two successive levels
93      REAL gamma                                ! Plume acceleration term (to compute vertical velocity)
94      REAL test                                 ! Test to know how to compute entrainment and detrainment
95       
96      REAL psat                                 ! Dummy argument for Psat_water()
97     
98      LOGICAL active(ngrid)                     ! If the plume is active at ig (speed and incoming mass flux > 0 or l=lmin)
99      LOGICAL activetmp(ngrid)                  ! If the plume is active at ig (active=true and outgoing mass flux > 0)
100     
101!===============================================================================
102! Initialization
103!===============================================================================
104     
105      zbetalpha = betalpha / (1. + betalpha)
106     
107      ztva(:,:)        = ztv(:,:)                                                ! ztva     is set to TPV environment
108      zhla(:,:)        = zhl(:,:)                                                ! zhla     is set to TP  environment
109      zqta(:,:)        = zqt(:,:)                                                ! zqta     is set to qt  environment
110      zqla(:,:)        = zql(:,:)                                                ! zqla     is set to ql  environment
111     
112      zqsa_est(:)      = 0.
113      zqsa(:,:)        = 0.
114     
115      zw2_est(:,:)     = 0.
116      zw2(:,:)         = 0.
117     
118      zbuoy(:,:)       = 0.
119     
120      f_star(:,:)      = 0.
121      detr_star(:,:)   = 0.
122      entr_star(:,:)   = 0.
123     
124      lmix(:)          = 1
125      lmin(:)          = 1
126     
127      ztv2(:,:)        = ztv(:,:)
128      ztv2(:,linf)     = ztv(:,linf) + d_temp
129     
130!===============================================================================
131! First layer computation
132!===============================================================================
133     
134      DO ig=1,ngrid
135         active(ig) = .false.
136         l = linf
137         DO WHILE (.not.active(ig).and.(pplev(ig,l+1).GT.pres_limit).and.(l.LT.nlay))
138            zdz = (zlev(ig,l+1) - zlev(ig,l))
139            zw2(ig,l+1) = 2. * afact * RG * zdz                               &  ! Do we have to divide by 1+betalpha (consider entrainment) ?
140            &           * (ztv2(ig,l) - ztv2(ig,l+1))                         &
141            &           / (ztv2(ig,l+1) * (1. + betalpha))
142            zw2m = zw2(ig,l+1) / 2.
143            entr_star(ig,l) = zdz * zbetalpha * (afact * RG * (ztv2(ig,l)     &
144            &               - ztv2(ig,l+1)) / ztv2(ig,l+1) / zw2m - fact_epsilon)
145            IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(entr_star(ig,l).GT.0.)) THEN
146               active(ig) = .true.
147               lmin(ig) = l
148               f_star(ig,l+1) = entr_star(ig,l)
149               zw2_est(ig,l+1) = zw2(ig,l+1)
150            ELSE
151               zw2(ig,l+1) = 0.
152               entr_star(ig,l) = 0.
153            ENDIF
154            l = l + 1
155         ENDDO
156      ENDDO
157     
158!===============================================================================
159! Thermal plumes computations
160!===============================================================================
161     
162      DO l=2,nlay-1
163         
164!-------------------------------------------------------------------------------
165! Check if thermal plume is (still) active
166!-------------------------------------------------------------------------------
167         
168!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
169! AB: we decide here if the plume is still active or not. When the plume's
170!     first level is reached, we set active to "true". Otherwise, it is given
171!     by zw2, f_star and entr_star.
172!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
173         DO ig=1,ngrid
174            IF (l==lmin(ig)+1) THEN
175               active(ig) = .true.
176            ENDIF
177           
178            active(ig) = active(ig)                                           &
179            &          .and. (zw2(ig,l).GT.1.e-10)                            &
180            &          .and. (f_star(ig,l) + entr_star(ig,l)).GT.1.e-10
181         ENDDO
182         
183         ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
184         
185         DO ig=1,ngrid
186            CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
187         ENDDO
188         
189!-------------------------------------------------------------------------------
190! Vertical speed (before mixing between plume and environment)
191!-------------------------------------------------------------------------------
192         
193         DO ig=1,ngrid
194            IF (active(ig)) THEN
195               zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsa_est(ig))                ! zqla_est set to ql   plume
196               zta_est(ig,l)  = zhla(ig,l-1) * zpopsk(ig,l)                   &  ! zta_est  set to TR   plume
197               &              + RLvCp * zqla_est(ig,l)
198               ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l)                  &  ! ztva_est set to TRPV plume
199               &              * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
200               
201               zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
202               zdz = zlev(ig,l+1) - zlev(ig,l)
203               
204!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
205! AB: initial formulae
206!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207!               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
208!               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
209!               zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
210!               zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2)
211!               zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2)
212!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213! AB: own derivation for zw2_est (Rio et al. 2010)
214!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
215!               zw2fact = 2. * fact_epsilon * zdz
216!               zdw2 = 2. * afact * zbuoy(ig,l) * zdz
217               zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
218               zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
219               zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2)
220            ENDIF
221         ENDDO
222         
223!-------------------------------------------------------------------------------
224! Mass flux, entrainment and detrainment
225!-------------------------------------------------------------------------------
226         
227         DO ig=1,ngrid
228            IF (active(ig)) THEN
229               
230               zdz = zlev(ig,l+1) - zlev(ig,l)
231               zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.                     ! AB: est-ce la bonne vitesse a utiliser ?
232               gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
233               
234               IF (zw2_est(ig,l).GT.0.) THEN
235                  test = gamma / zw2m - nu
236               ELSE
237                  print *, 'ERROR: zw2_est is negative while plume is active!'
238                  print *, 'ig,l', ig, l
239                  print *, 'zw2_est', zw2_est(ig,l)
240                  call abort
241               ENDIF
242               
243               IF (test.gt.0.) THEN
244                  detr_star(ig,l) = zdz * nu
245                  entr_star(ig,l) = zdz * (zbetalpha * gamma / zw2m + nu)
246               ELSE
247                  detr_star(ig,l) = zdz * (nu - betalpha * gamma / zw2m)
248                  entr_star(ig,l) = zdz * nu
249               ENDIF
250               
251!               IF (detr_star(ig,l).lt.0.) THEN
252!                  print *, 'WARNING: detrainment is negative!'
253!                  print *, 'l,detr', l, detr_star(ig,l)
254!               ENDIF
255               
256!               IF (entr_star(ig,l).lt.0.) THEN
257!                  print *, 'WARNING: entrainment is negative!'
258!                  print *, 'l,entr', l, entr_star(ig,l)
259!               ENDIF
260               
261               f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
262               
263!               IF (f_star(ig,l+1).le.0.) THEN
264!                  print *, 'WARNING: mass flux is negative!'
265!                  print *, 'l,f_star', l+1, f_star(ig,l+1)
266!               ENDIF
267               
268            ENDIF
269         ENDDO
270         
271!-------------------------------------------------------------------------------
272! Mixing between thermal plume and environment
273!-------------------------------------------------------------------------------
274         
275         activetmp(:) = active(:) .and. (f_star(:,l+1).GT.1.e-10)
276         
277         DO ig=1,ngrid
278            IF (activetmp(ig)) THEN
279               zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1)                      &  ! zhla is set to TP in plume (mixed)
280               &          + entr_star(ig,l) * zhl(ig,l))                      &
281               &          / (f_star(ig,l+1) + detr_star(ig,l))
282               zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) +                    &  ! zqta is set to qt in plume (mixed)
283               &          + entr_star(ig,l) * zqt(ig,l))                      &
284               &          / (f_star(ig,l+1) + detr_star(ig,l))
285            ENDIF
286         ENDDO
287         
288         ztemp(:) = zpopsk(:,l) * zhla(:,l)
289         
290         DO ig=1,ngrid
291            IF (activetmp(ig)) THEN
292               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
293            ENDIF
294         ENDDO
295         
296!-------------------------------------------------------------------------------
297! Vertical speed (after mixing between plume and environment)
298!-------------------------------------------------------------------------------
299         
300         DO ig=1,ngrid
301            IF (activetmp(ig)) THEN
302               zqla(ig,l) = max(0.,zqta(ig,l)-zqsa(ig,l))                        ! zqla is set to ql   plume (mixed)
303               
304               zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) + RLvCp * zqla(ig,l)        ! ztva is set to TR   plume (mixed)
305               ztva(ig,l) = zta(ig,l) / zpopsk(ig,l)                          &  ! ztva is set to TRPV plume (mixed)           
306               &          * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
307               
308               zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
309               zdz = zlev(ig,l+1) - zlev(ig,l)
310               
311!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
312! AB: initial formula
313!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
314!               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
315!               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
316!               zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
317!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318! AB: own derivation for zw2 (Rio et al. 2010)
319!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
320!               zw2fact = 2. * (fact_epsilon * zdz  + entr_star(ig,l) / f_star(ig,l))
321!               zdw2 = 2. * afact * zbuoy(ig,l) * zdz
322               zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
323               zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
324               zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)
325            ENDIF
326         ENDDO
327         
328      ENDDO
329     
330     
331RETURN
332END
Note: See TracBrowser for help on using the repository browser.