Context Navigation

thermcell_main.F90 @ 2229

Last change on this file since 2229 was 2229, checked in by aboissinot, 5 years ago

A bug in thermcell_dq is fixed. Now zqt is correctly initialized when tracer h2o_vap
is missing (consistency with flag water is assumed).

Useless arguments in thermcell_dq subroutine are removed (lmin, lmax)

File size: 18.1 KB

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

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

Download in other formats:

Original Format