Context Navigation

thermcell_plume.F90 @ 3638

Last change on this file since 3638 was 3435, checked in by alesaux, 14 months ago
Generic PCM: Bug fix in the Thermal Plumes Model for generic tracers. Error in computation of potential temperature. ALS
File size: 14.3 KB

Line
1	!
2	!
3	!
4	SUBROUTINE 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	USE comcstfi_mod, ONLY: r, cpp, mugaz
29	USE callkeys_mod, ONLY: water, generic_condensation
30	USE generic_cloud_common_h, ONLY: Psat_generic, epsi_generic, RLVTT_generic
31
32	IMPLICIT NONE
33
34
35	!===============================================================================
36	! Declaration
37	!===============================================================================
38
39	! Inputs:
40	! -------
41
42	INTEGER, INTENT(in) :: ngrid
43	INTEGER, INTENT(in) :: nlay
44	INTEGER, INTENT(in) :: nq
45
46	INTEGER, INTENT(in) :: lbot(ngrid) ! First considered layer
47
48	REAL, INTENT(in) :: ptimestep
49	REAL, INTENT(in) :: zlev(ngrid,nlay+1) ! Levels altitude
50	REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Levels pressure
51	REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function
52
53	REAL, INTENT(in) :: ztv(ngrid,nlay) ! TRPV environment
54	REAL, INTENT(in) :: zhl(ngrid,nlay) ! TP environment
55	REAL, INTENT(in) :: zqt(ngrid,nlay) ! qt environment
56	REAL, INTENT(in) :: zql(ngrid,nlay) ! ql environment
57
58	! Outputs:
59	! --------
60
61	INTEGER, INTENT(out) :: lmin(ngrid) ! Plume bottom level (first unstable level)
62
63	REAL, INTENT(out) :: detr_star(ngrid,nlay) ! Normalized detrainment
64	REAL, INTENT(out) :: entr_star(ngrid,nlay) ! Normalized entrainment
65	REAL, INTENT(out) :: f_star(ngrid,nlay+1) ! Normalized mass flux
66
67	REAL, INTENT(out) :: ztva(ngrid,nlay) ! TRPV plume (after mixing)
68	REAL, INTENT(out) :: zhla(ngrid,nlay) ! TP plume (after mixing)
69	REAL, INTENT(out) :: zqla(ngrid,nlay) ! ql plume (after mixing)
70	REAL, INTENT(out) :: zqta(ngrid,nlay) ! qt plume (after mixing)
71	REAL, INTENT(out) :: zqsa(ngrid,nlay) ! qsat plume (after mixing)
72	REAL, INTENT(out) :: zw2(ngrid,nlay+1) ! w plume (after mixing)
73
74	! Local:
75	! ------
76
77	INTEGER ig, l, k
78	INTEGER l_start
79
80	REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing)
81	REAL zqla_est(ngrid,nlay) ! ql plume (before mixing)
82	REAL zta_est(ngrid,nlay) ! TR plume (before mixing)
83	REAL zqsa_est(ngrid) ! qsat plume (before mixing)
84	REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing)
85
86	REAL zta(ngrid,nlay) ! TR plume (after mixing)
87
88	REAL zbuoy(ngrid,nlay) ! Plume buoyancy
89	REAL ztemp(ngrid) ! Temperature to compute saturation vapor pressure
90	REAL zdz ! Layers heights
91	REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf)
92
93	REAL zdw2 !
94	REAL zw2fact !
95	REAL zw2m ! Average vertical velocity between two successive levels
96	REAL gamma ! Plume acceleration term (to compute vertical velocity)
97	REAL test !
98
99	REAL psat ! Dummy argument for Psat_water()
100
101	LOGICAL active(ngrid) ! If the plume is active (speed and incoming mass flux > 0)
102	LOGICAL activetmp(ngrid) ! If the plume is active (active=true and outgoing mass flux > 0)
103
104	REAL, SAVE :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
105	REAL :: RV_generic
106	REAL :: RETV_comp, RLvCp_comp !values used for computation (depends if water or generic tracer)
107
108	!===============================================================================
109	! Initialization
110	!===============================================================================
111
112	ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment
113	zhla(:,:) = zhl(:,:) ! zhla is set to TP environment
114	zqta(:,:) = zqt(:,:) ! zqta is set to qt environment
115	zqla(:,:) = zql(:,:) ! zqla is set to ql environment
116	zqsa(:,:) = 0.
117	zw2(:,:) = 0.
118
119	ztva_est(:,:) = ztv(:,:) ! ztva_est is set to TPV environment
120	zqla_est(:,:) = zql(:,:) ! zqla_est is set to ql environment
121	zqsa_est(:) = 0.
122	zw2_est(:,:) = 0.
123
124	zbuoy(:,:) = 0.
125
126	f_star(:,:) = 0.
127	detr_star(:,:) = 0.
128	entr_star(:,:) = 0.
129
130	lmin(:) = lbot(:)
131
132	ztv2(:,:) = ztv(:,:)
133	ztv2(:,linf) = ztv(:,linf) + d_temp
134
135	active(:) = .false.
136
137	l_start = nlay
138
139	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
140
141	! ALS24 for thermal plume model with generic tracer
142	IF (water) THEN
143	RETV_comp = RETV
144	RLvCp_comp = RLvCp
145	ELSEIF (generic_condensation .AND. .NOT. water ) THEN
146	RV_generic = (8.3145111000.)/(epsi_genericmugaz)
147	RETV_comp = RV_generic/r-1.
148	RLvCp_comp = RLVTT_generic/cpp
149	ENDIF
150
151	!===============================================================================
152	! First layer computation
153	!===============================================================================
154
155	DO ig=1,ngrid
156	l = lbot(ig)
157	l_start = MIN(l_start, lbot(ig)+1)
158	DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
159	zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
160	IF (zbuoy(ig,l) > 0.) THEN
161	lmin(ig) = l
162	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
163	! AB: entrainement and mass flux initial values are set to 1. The physical value
164	! will be computed thanks to the closure relation in thermcell_closure.
165	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
166	entr_star(ig,l) = 1.
167	f_star(ig,l+1) = 1.
168	zdz = zlev(ig,l+1) - zlev(ig,l)
169	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
170	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
171	zw2_est(ig,l+1) = exp(-zw2fact) * zdw2
172	zw2(ig,l+1) = zw2_est(ig,l+1)
173	active(ig) = .true.
174	ENDIF
175	l = l + 1
176	ENDDO
177	ENDDO
178
179	!===============================================================================
180	! Thermal plumes computations
181	!===============================================================================
182
183	DO l=l_start,nlay-1
184
185	!-------------------------------------------------------------------------------
186	! Is thermal plume (still) active ?
187	!-------------------------------------------------------------------------------
188
189	DO ig=1,ngrid
190	active(ig) = (zw2(ig,l) > 1.e-9).and.(f_star(ig,l) > 1.e-9)
191	ENDDO
192
193	!-------------------------------------------------------------------------------
194	! Latent heat release (before mixing)
195	!-------------------------------------------------------------------------------
196
197	ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
198
199	DO ig=1,ngrid
200	IF (active(ig)) THEN
201	IF (water) THEN
202	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
203	ELSEIF (generic_condensation .AND. .NOT. water) THEN
204	CALL Psat_generic(ztemp(ig),pplev(ig,l),metallicity,psat,zqsa_est(ig))
205	ENDIF
206	ENDIF
207	ENDDO
208
209	!-------------------------------------------------------------------------------
210	! Vertical speed (before mixing)
211	!-------------------------------------------------------------------------------
212
213	DO ig=1,ngrid
214	IF (active(ig)) THEN
215	zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est is set to ql plume
216	zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est is set to TR plume
217	& + RLvCp_comp * zqla_est(ig,l)
218	ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est is set to TRPV plume
219	& * (1. + RETV_comp * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
220
221	zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
222	zdz = zlev(ig,l+1) - zlev(ig,l)
223
224	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
225	zdw2 = afact * zbuoy(ig,l) / fact_epsilon
226	zw2_est(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2)
227	ENDIF
228	ENDDO
229
230	!-------------------------------------------------------------------------------
231	! Mass flux, entrainment and detrainment
232	!-------------------------------------------------------------------------------
233
234	DO ig=1,ngrid
235	IF (active(ig)) THEN
236
237	zdz = zlev(ig,l+1) - zlev(ig,l)
238	zw2m = (zw2_est(ig,l+1) + zw2(ig,l)) / 2.
239	gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
240
241	IF (zw2m > 0.) THEN
242	test = gamma / zw2m - nu
243	ELSE
244	test = 0.
245	print *, 'WARNING: vertical speed is negative while plume is active!'
246	print *, 'ig,l', ig, l
247	print *, 'zw2m', zw2m
248	ENDIF
249
250	IF (test > 0.) THEN
251	detr_star(ig,l) = zdz * f_star(ig,l) * nu
252	entr_star(ig,l) = zdz * f_star(ig,l) * (betalpha * gamma / zw2m + nu) / (betalpha + 1)
253	ELSE
254	detr_star(ig,l) = zdz * f_star(ig,l) * ((betalpha + 1) * nu - betalpha * gamma / zw2m)
255	entr_star(ig,l) = zdz * f_star(ig,l) * nu
256	ENDIF
257
258	f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
259
260	ENDIF
261	ENDDO
262
263	!-------------------------------------------------------------------------------
264	! Mixing between thermal plume and environment
265	!-------------------------------------------------------------------------------
266
267	activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-9)
268
269	DO ig=1,ngrid
270	IF (activetmp(ig)) THEN
271	zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1) & ! zhla is set to TP in plume (mixed)
272	& + entr_star(ig,l) * zhl(ig,l)) &
273	& / (f_star(ig,l+1) + detr_star(ig,l))
274	zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) & ! zqta is set to qt in plume (mixed)
275	& + entr_star(ig,l) * zqt(ig,l)) &
276	& / (f_star(ig,l+1) + detr_star(ig,l))
277	ENDIF
278	ENDDO
279
280	!-------------------------------------------------------------------------------
281	! Latent heat release (after mixing)
282	!-------------------------------------------------------------------------------
283
284	ztemp(:) = zpopsk(:,l) * zhla(:,l)
285
286	DO ig=1,ngrid
287	IF (activetmp(ig)) THEN
288	IF (water) THEN
289	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
290	ELSEIF (generic_condensation .AND. .NOT. water ) THEN
291	CALL Psat_generic(ztemp(ig),pplev(ig,l),metallicity,psat,zqsa(ig,l))
292	ENDIF
293	ENDIF
294	ENDDO
295
296	!-------------------------------------------------------------------------------
297	! Vertical speed (after mixing)
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	zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed)
304	& + RLvCp_comp * zqla(ig,l)
305	ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed)
306	& * (1. + RETV_comp*(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	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
312	zdw2 = afact * zbuoy(ig,l) / fact_epsilon
313	zw2(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2)
314	ENDIF
315	ENDDO
316
317	ENDDO
318
319
320	RETURN
321	END

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

Download in other formats:

Original Format