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

Last change on this file since 2694 was 2470, checked in by yjaziri, 5 years ago

Generic GCM:
global1d and szangle for 1D simulation moved from callcorrk to callkeys
to defined a consistent 1D sza in physiq_mod used also in chemistry
+ typo
YJ
File size: 13.0 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep,                           &
5                           ztv,zhl,zqt,zql,zlev,pplev,zpopsk,                 &
6                           detr_star,entr_star,f_star,                        &
7                           ztva,zhla,zqta,zqla,zqsa,                          &
8                           zw2,lbot,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, INTENT(in) :: ngrid
40      INTEGER, INTENT(in) :: nlay
41      INTEGER, INTENT(in) :: nq
42     
43      INTEGER, INTENT(in) :: lbot(ngrid)              ! First considered layer
44     
45      REAL, INTENT(in) :: ptimestep
46      REAL, INTENT(in) :: zlev(ngrid,nlay+1)          ! Levels altitude
47      REAL, INTENT(in) :: pplev(ngrid,nlay+1)         ! Levels pressure
48      REAL, INTENT(in) :: zpopsk(ngrid,nlay)          ! Exner function
49     
50      REAL, INTENT(in) :: ztv(ngrid,nlay)             ! TRPV environment
51      REAL, INTENT(in) :: zhl(ngrid,nlay)             ! TP   environment
52      REAL, INTENT(in) :: zqt(ngrid,nlay)             ! qt   environment
53      REAL, INTENT(in) :: zql(ngrid,nlay)             ! ql   environment
54     
55!     Outputs:
56!     --------
57     
58      INTEGER, INTENT(out) :: lmin(ngrid)             ! Plume bottom level (first unstable level)
59     
60      REAL, INTENT(out) :: detr_star(ngrid,nlay)      ! Normalized detrainment
61      REAL, INTENT(out) :: entr_star(ngrid,nlay)      ! Normalized entrainment
62      REAL, INTENT(out) :: f_star(ngrid,nlay+1)       ! Normalized mass flux
63     
64      REAL, INTENT(out) :: ztva(ngrid,nlay)           ! TRPV plume (after mixing)
65      REAL, INTENT(out) :: zhla(ngrid,nlay)           ! TP   plume (after mixing)
66      REAL, INTENT(out) :: zqla(ngrid,nlay)           ! ql   plume (after mixing)
67      REAL, INTENT(out) :: zqta(ngrid,nlay)           ! qt   plume (after mixing)
68      REAL, INTENT(out) :: zqsa(ngrid,nlay)           ! qsat plume (after mixing)
69      REAL, INTENT(out) :: zw2(ngrid,nlay+1)          ! w    plume (after mixing)
70     
71!     Local:
72!     ------
73     
74      INTEGER ig, l, k
75      INTEGER l_start
76     
77      REAL ztva_est(ngrid,nlay)                       ! TRPV plume (before mixing)
78      REAL zqla_est(ngrid,nlay)                       ! ql   plume (before mixing)
79      REAL zta_est(ngrid,nlay)                        ! TR   plume (before mixing)
80      REAL zqsa_est(ngrid)                            ! qsat plume (before mixing)
81      REAL zw2_est(ngrid,nlay+1)                      ! w    plume (before mixing)
82     
83      REAL zta(ngrid,nlay)                            ! TR   plume (after mixing)
84     
85      REAL zbuoy(ngrid,nlay)                          ! Plume buoyancy
86      REAL ztemp(ngrid)                               ! Temperature to compute saturation vapor pressure
87      REAL zdz                                        ! Layers heights
88      REAL ztv2(ngrid,nlay)                           ! ztv + d_temp * Dirac(l=linf)
89     
90      REAL zdw2                                       !
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                                       !
95     
96      REAL psat                                       ! Dummy argument for Psat_water()
97     
98      LOGICAL active(ngrid)                           ! If the plume is active (speed and incoming mass flux > 0)
99      LOGICAL activetmp(ngrid)                        ! If the plume is active (active=true and outgoing mass flux > 0)
100     
101!===============================================================================
102! Initialization
103!===============================================================================
104     
105      ztva(:,:) = ztv(:,:)                                                       ! ztva is set to TPV environment
106      zhla(:,:) = zhl(:,:)                                                       ! zhla is set to TP  environment
107      zqta(:,:) = zqt(:,:)                                                       ! zqta is set to qt  environment
108      zqla(:,:) = zql(:,:)                                                       ! zqla is set to ql  environment
109      zqsa(:,:) = 0.
110      zw2(:,:)  = 0.
111     
112      ztva_est(:,:) = ztv(:,:)                                                   ! ztva_est is set to TPV environment
113      zqla_est(:,:) = zql(:,:)                                                   ! zqla_est is set to ql  environment
114      zqsa_est(:)   = 0.
115      zw2_est(:,:)  = 0.
116     
117      zbuoy(:,:) = 0.
118     
119      f_star(:,:)    = 0.
120      detr_star(:,:) = 0.
121      entr_star(:,:) = 0.
122     
123      lmin(:) = lbot(:)
124     
125      ztv2(:,:)    = ztv(:,:)
126      ztv2(:,linf) = ztv(:,linf) + d_temp
127     
128      active(:) = .false.
129     
130      l_start = nlay
131     
132!===============================================================================
133! First layer computation
134!===============================================================================
135     
136      DO ig=1,ngrid
137         l = lbot(ig)
138         l_start = MIN(l_start, lbot(ig)+1)
139         DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
140            zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
141            IF (zbuoy(ig,l) > 0.) THEN
142               lmin(ig) = l
143!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144! AB: entrainement and mass flux initial values are set to 1. The physical value
145!     will be computed thanks to the closure relation in thermcell_closure.
146!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
147               entr_star(ig,l) = 1.
148               f_star(ig,l+1) = 1.
149               zdz = zlev(ig,l+1) - zlev(ig,l)
150               zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
151               zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
152               zw2_est(ig,l+1) = exp(-zw2fact) * zdw2
153               zw2(ig,l+1) = zw2_est(ig,l+1)
154               active(ig) = .true.
155            ENDIF
156            l = l + 1
157         ENDDO
158      ENDDO
159     
160!===============================================================================
161! Thermal plumes computations
162!===============================================================================
163     
164      DO l=l_start,nlay-1
165         
166!-------------------------------------------------------------------------------
167! Is thermal plume (still) active ?
168!-------------------------------------------------------------------------------
169         
170         DO ig=1,ngrid
171            active(ig) = (zw2(ig,l) > 1.e-9).and.(f_star(ig,l) > 1.e-9)
172         ENDDO
173         
174!-------------------------------------------------------------------------------
175! Latent heat release (before mixing)
176!-------------------------------------------------------------------------------
177         
178         ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
179         
180         DO ig=1,ngrid
181            IF (active(ig)) THEN
182               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
183            ENDIF
184         ENDDO
185         
186!-------------------------------------------------------------------------------
187! Vertical speed (before mixing)
188!-------------------------------------------------------------------------------
189         
190         DO ig=1,ngrid
191            IF (active(ig)) THEN
192               zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig))              ! zqla_est is set to ql   plume
193               zta_est(ig,l)  = zhla(ig,l-1) * zpopsk(ig,l)                   &  ! zta_est  is set to TR   plume
194               &              + RLvCp * zqla_est(ig,l)
195               ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l)                  &  ! ztva_est is set to TRPV plume
196               &              * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
197               
198               zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
199               zdz = zlev(ig,l+1) - zlev(ig,l)
200               
201               zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
202               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
203               zw2_est(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2)
204            ENDIF
205         ENDDO
206         
207!-------------------------------------------------------------------------------
208! Mass flux, entrainment and detrainment
209!-------------------------------------------------------------------------------
210         
211         DO ig=1,ngrid
212            IF (active(ig)) THEN
213               
214               zdz = zlev(ig,l+1) - zlev(ig,l)
215               zw2m = (zw2_est(ig,l+1) + zw2(ig,l)) / 2.
216               gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
217               
218               IF (zw2m > 0.) THEN
219                  test = gamma / zw2m - nu
220               ELSE
221                  test = 0.
222                  print *, 'WARNING: vertical speed is negative while plume is active!'
223                  print *, 'ig,l', ig, l
224                  print *, 'zw2m', zw2m
225               ENDIF
226               
227               IF (test > 0.) THEN
228                  detr_star(ig,l) = zdz * f_star(ig,l) * nu
229                  entr_star(ig,l) = zdz * f_star(ig,l) * (betalpha * gamma / zw2m + nu) / (betalpha + 1)
230               ELSE
231                  detr_star(ig,l) = zdz * f_star(ig,l) * ((betalpha + 1) * nu - betalpha * gamma / zw2m)
232                  entr_star(ig,l) = zdz * f_star(ig,l) * nu
233               ENDIF
234               
235               f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
236               
237            ENDIF
238         ENDDO
239         
240!-------------------------------------------------------------------------------
241! Mixing between thermal plume and environment
242!-------------------------------------------------------------------------------
243         
244         activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-9)
245         
246         DO ig=1,ngrid
247            IF (activetmp(ig)) THEN
248               zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1)                      &  ! zhla is set to TP in plume (mixed)
249               &          + entr_star(ig,l) * zhl(ig,l))                      &
250               &          / (f_star(ig,l+1) + detr_star(ig,l))
251               zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1)                      &  ! zqta is set to qt in plume (mixed)
252               &          + entr_star(ig,l) * zqt(ig,l))                      &
253               &          / (f_star(ig,l+1) + detr_star(ig,l))
254            ENDIF
255         ENDDO
256         
257!-------------------------------------------------------------------------------
258! Latent heat release (after mixing)
259!-------------------------------------------------------------------------------
260         
261         ztemp(:) = zpopsk(:,l) * zhla(:,l)
262         
263         DO ig=1,ngrid
264            IF (activetmp(ig)) THEN
265               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
266            ENDIF
267         ENDDO
268         
269!-------------------------------------------------------------------------------
270! Vertical speed (after mixing)
271!-------------------------------------------------------------------------------
272         
273         DO ig=1,ngrid
274            IF (activetmp(ig)) THEN
275               zqla(ig,l) = MAX(0.,zqta(ig,l) - zqsa(ig,l))                      ! zqla is set to ql   plume (mixed)
276               zta(ig,l)  = zhla(ig,l) * zpopsk(ig,l)                         &  ! ztva is set to TR   plume (mixed)
277               &          + RLvCp * zqla(ig,l)
278               ztva(ig,l) = zta(ig,l) / zpopsk(ig,l)                          &  ! ztva is set to TRPV plume (mixed)
279               &          * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
280               
281               zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
282               zdz = zlev(ig,l+1) - zlev(ig,l)
283               
284               zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
285               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
286               zw2(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2)
287            ENDIF
288         ENDDO
289         
290      ENDDO
291     
292     
293RETURN
294END
Note: See TracBrowser for help on using the repository browser.