Context Navigation

thermcell_plume.F90 @ 2176

Last change on this file since 2176 was 2145, checked in by aboissinot, 5 years ago
Some temporary outputs in thermcell_closure are removed and a bug in thermcell_plume is fixed.
File size: 14.5 KB

Line
1	!
2	!
3	!
4	SUBROUTINE 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,zqsa, &
8	zw2,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 ?
65	REAL zqla(ngrid,nlay) ! ql plume (after mixing)
66	REAL zqta(ngrid,nlay) ! qt plume ?
67	REAL zqsa(ngrid,nlay) ! qsat plume (after mixing)
68	REAL zw2(ngrid,nlay+1) ! w 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) ! w 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	active(:) = .false.
131
132	!===============================================================================
133	! First layer computation
134	!===============================================================================
135
136	DO ig=1,ngrid
137	l = linf
138	DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
139	zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
140	zdz = zlev(ig,l+1) - zlev(ig,l)
141	zw2m = afact * zbuoy(ig,l) * zdz / (1. + betalpha)
142	! gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
143	! test = gamma / zw2m - nu
144	test = zbuoy(ig,l)
145	IF (test > 0.) THEN
146	lmin(ig) = l
147	! entr_star(ig,l) = zdz * f_star(ig,l) * zbetalpha * gamma / zw2m - nu ! Problem because f*(ig,l) = 0
148	! detr_star(ig,l) = f_star(ig,l) * nu ! Problem because f*(ig,l) = 0
149	! f_star(ig,l+1) = entr_star(ig,l) - detr_star(ig,l)
150	entr_star(ig,l) = 1.
151	f_star(ig,l+1) = 1.
152	zw2_est(ig,l+1) = zw2m * 2.
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=2,nlay-1
165
166	!-------------------------------------------------------------------------------
167	! Is thermal plume (still) active ?
168	!-------------------------------------------------------------------------------
169
170	DO ig=1,ngrid
171	active(ig) = (active(ig).or.(l == lmin(ig)+1)) &
172	& .and.(zw2(ig,l) > 1.e-10) &
173	& .and.(f_star(ig,l) > 1.e-10)
174	ENDDO
175
176	!-------------------------------------------------------------------------------
177	! Latent heat release (before mixing)
178	!-------------------------------------------------------------------------------
179
180	ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
181
182	DO ig=1,ngrid
183	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
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 set to ql plume
193	zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est set to TR plume
194	& + RLvCp * zqla_est(ig,l)
195	ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est 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	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202	! AB: initial formulae
203	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
204	! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
205	! zdw2 = afact * zbuoy(ig,l) / fact_epsilon
206	! zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
207	! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2)
208	! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2)
209	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
210	! AB: own derivation for zw2_est (Rio et al. 2010)
211	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
212	! zw2fact = 2. * fact_epsilon * zdz
213	! zdw2 = 2. * afact * zbuoy(ig,l) * zdz
214	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
215	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
216	zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2)
217	ENDIF
218	ENDDO
219
220	!-------------------------------------------------------------------------------
221	! Mass flux, entrainment and detrainment
222	!-------------------------------------------------------------------------------
223
224	DO ig=1,ngrid
225	IF (active(ig)) THEN
226
227	zdz = zlev(ig,l+1) - zlev(ig,l)
228	zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.
229	gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
230
231	IF (zw2_est(ig,l) > 0.) THEN
232	test = gamma / zw2m - nu
233	ELSE
234	print *, 'ERROR: zw2_est is negative while plume is active!'
235	print *, 'ig,l', ig, l
236	print *, 'zw2_est', zw2_est(ig,l)
237	call abort
238	ENDIF
239
240	IF (test > 0.) THEN
241	detr_star(ig,l) = zdz * f_star(ig,l) * nu
242	entr_star(ig,l) = zdz * f_star(ig,l) * (zbetalpha * gamma / zw2m + nu)
243	ELSE
244	detr_star(ig,l) = zdz * f_star(ig,l) * (nu - betalpha * gamma / zw2m)
245	entr_star(ig,l) = zdz * f_star(ig,l) * nu
246	ENDIF
247
248	f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
249
250	ENDIF
251	ENDDO
252
253	!-------------------------------------------------------------------------------
254	! Mixing between thermal plume and environment
255	!-------------------------------------------------------------------------------
256
257	activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-10)
258
259	DO ig=1,ngrid
260	IF (activetmp(ig)) THEN
261	zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1) & ! zhla is set to TP in plume (mixed)
262	& + entr_star(ig,l) * zhl(ig,l)) &
263	& / (f_star(ig,l+1) + detr_star(ig,l))
264	zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) + & ! zqta is set to qt in plume (mixed)
265	& + entr_star(ig,l) * zqt(ig,l)) &
266	& / (f_star(ig,l+1) + detr_star(ig,l))
267	ENDIF
268	ENDDO
269
270	!-------------------------------------------------------------------------------
271	! Latent heat release (after mixing)
272	!-------------------------------------------------------------------------------
273
274	ztemp(:) = zpopsk(:,l) * zhla(:,l)
275
276	DO ig=1,ngrid
277	IF (activetmp(ig)) THEN
278	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
279	ENDIF
280	ENDDO
281
282	!-------------------------------------------------------------------------------
283	! Vertical speed (after mixing)
284	!-------------------------------------------------------------------------------
285
286	DO ig=1,ngrid
287	IF (activetmp(ig)) THEN
288	zqla(ig,l) = MAX(0.,zqta(ig,l) - zqsa(ig,l)) ! zqla is set to ql plume (mixed)
289	zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed)
290	& + RLvCp * zqla(ig,l)
291	ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed)
292	& * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
293
294	zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
295	zdz = zlev(ig,l+1) - zlev(ig,l)
296
297	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
298	! AB: initial formula
299	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
300	! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
301	! zdw2 = afact * zbuoy(ig,l) / fact_epsilon
302	! zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
303	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
304	! AB: own derivation for zw2 (Rio et al. 2010)
305	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
306	! zw2fact = 2. * (fact_epsilon * zdz + entr_star(ig,l) / f_star(ig,l))
307	! zdw2 = 2. * afact * zbuoy(ig,l) * zdz
308	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
309	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
310	zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)
311	ENDIF
312	ENDDO
313
314	ENDDO
315
316
317	RETURN
318	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format