Context Navigation

lmdz_thermcell_main.F90 @ 5119

Last change on this file since 5119 was 5117, checked in by abarral, 4 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 36.5 KB

Line
1	MODULE lmdz_thermcell_main
2	! $Id: lmdz_thermcell_main.F90 5117 2024-07-24 14:23:34Z abarral $
3
4	! A REGARDER !!!!!!!!!!!!!!!!!
5	! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)
6	! ATTENTION : dans thermcell_env, on condense potentiellement de l'eau. Mais comme on ne mélange pas l'eau liquide supposant qu'il n'y en n'a pas, c'est potentiellement un souci
7	CONTAINS
8
9	SUBROUTINE thermcell_main
10	(itap, ngrid, nlay, ptimestep &
11	, pplay,pplev, pphi, debut &
12	, puwind, pvwind,ptemp, p_o, ptemp_env, po_env &
13	, pduadj,pdvadj, pdtadj, pdoadj &
14	, fm0, entr0,detr0, zqta, zqla, lmax &
15	, ratqscth,ratqsdiff, zqsatth &
16	, zmax0, f0, zw2,fraca, ztv &
17	, zpspsk, ztla, zthl,ztva &
18	, pcon, rhobarz, wth3, wmax_sec,lalim, fm, alim_star, zmax, zcong &
19	#ifdef ISO
20	,xtpo,xtpdoadj &
21	#endif
22	)
23
24
25
26	! USE necessaires pour les lignes importees de thermcell_env
27	USE lmdz_thermcell_ini, ONLY: thermcell_ini, dqimpl, dvdq, prt_level, lunout, prt_level
28	USE lmdz_thermcell_ini, ONLY: iflag_thermals_closure, iflag_thermals_ed, tau_thermals, r_aspect_thermals
29	USE lmdz_thermcell_ini, ONLY: iflag_thermals_down, fact_thermals_down
30	USE lmdz_thermcell_ini, ONLY: iflag_thermals_tenv
31	USE lmdz_thermcell_ini, ONLY: RD, RG
32
33	USE lmdz_thermcell_down, ONLY: thermcell_updown_dq
34	USE lmdz_thermcell_closure, ONLY: thermcell_closure
35	USE lmdz_thermcell_dq, ONLY: thermcell_dq
36	USE lmdz_thermcell_dry, ONLY: thermcell_dry
37	USE lmdz_thermcell_dv2, ONLY: thermcell_dv2
38	USE lmdz_thermcell_env, ONLY: thermcell_env
39	USE lmdz_thermcell_flux2, ONLY: thermcell_flux2
40	USE lmdz_thermcell_height, ONLY: thermcell_height
41	USE lmdz_thermcell_plume, ONLY: thermcell_plume
42	USE lmdz_thermcell_plume_6A, ONLY: thermcell_plume_6A, thermcell_plume_5B
43
44	! USE necessaires pour les lignes importees de thermcell_env
45	USE lmdz_thermcell_ini, ONLY: RLvCp, RKAPPA, RETV
46	USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
47	USE lmdz_abort_physic, ONLY: abort_physic
48
49
50	#ifdef ISO
51	USE infotrac_phy, ONLY: ntiso
52	#ifdef ISOVERIF
53	USE isotopes_mod, ONLY: iso_eau,iso_HDO
54	USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
55	iso_verif_aberrant_encadre
56	#endif
57	#endif
58
59
60	IMPLICIT NONE
61
62	!=======================================================================
63	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
64	! Version du 09.02.07
65	! Calcul du transport vertical dans la couche limite en presence
66	! de "thermiques" explicitement representes avec processus nuageux
67
68	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
69
70	! le thermique est suppose homogene et dissipe par melange avec
71	! son environnement. la longueur l_mix controle l'efficacite du
72	! melange
73
74	! Le calcul du transport des differentes especes se fait en prenant
75	! en compte:
76	! 1. un flux de masse montant
77	! 2. un flux de masse descendant
78	! 3. un entrainement
79	! 4. un detrainement
80
81	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
82	! Introduction of an implicit computation of vertical advection in
83	! the environment of thermal plumes in thermcell_dq
84	! impl = 0 : explicit, 1 : implicit, -1 : old version
85	! controled by iflag_thermals =
86	! 15, 16 run with impl=-1 : numerical convergence with NPv3
87	! 17, 18 run with impl=1 : more stable
88	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
89
90	! Using
91	! abort_physic
92	! iso_verif_aberrant_encadre
93	! iso_verif_egalite
94	! test_ltherm
95	! thermcell_closure
96	! thermcell_dq
97	! thermcell_dry
98	! thermcell_dv2
99	! thermcell_env
100	! thermcell_flux2
101	! thermcell_height
102	! thermcell_plume
103	! thermcell_plume_5B
104	! thermcell_plume_6A
105
106	!=======================================================================
107
108
109	!-----------------------------------------------------------------------
110	! declarations:
111	! -------------
112
113
114	! arguments:
115	! ----------
116	INTEGER, INTENT(IN) :: itap, ngrid, nlay
117	REAL, INTENT(IN) :: ptimestep
118	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: ptemp, puwind, pvwind, pplay, pphi, ptemp_env, po_env
119	! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)
120	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: p_o
121	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zpspsk
122	REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: pplev
123	INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmax
124	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: pdtadj, pduadj, pdvadj, pdoadj, entr0, detr0
125	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztla, zqla, zqta, zqsatth, zthl
126	REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: fm0, zw2, fraca
127	REAL, INTENT(INOUT), DIMENSION(ngrid) :: zmax0, f0
128	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva, ztv
129	logical, INTENT(IN) :: debut
130	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ratqscth, ratqsdiff
131
132	REAL, INTENT(OUT), DIMENSION(ngrid) :: pcon
133	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: rhobarz, wth3
134	REAL, INTENT(OUT), DIMENSION(ngrid) :: wmax_sec
135	INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lalim
136	REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: fm
137	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: alim_star
138	REAL, INTENT(OUT), DIMENSION(ngrid) :: zmax, zcong
139
140	! local:
141	! ------
142
143	INTEGER, save :: igout = 1
144	!$OMP THREADPRIVATE(igout)
145	INTEGER, save :: lunout1 = 6
146	!$OMP THREADPRIVATE(lunout1)
147	INTEGER, save :: lev_out = 10
148	!$OMP THREADPRIVATE(lev_out)
149
150	REAL lambda, zf, zf2, var, vardiff, CHI
151	INTEGER ig, k, l, ierr, ll
152	LOGICAL sorties
153	REAL, DIMENSION(ngrid) :: linter, zmix, zmax_sec, lintercong
154	INTEGER, DIMENSION(ngrid) :: lmin, lmix, lmix_bis, nivcon, lcong
155	REAL, DIMENSION(ngrid, nlay) :: ztva_est
156	REAL, DIMENSION(ngrid, nlay) :: deltaz, zlay, zdthladj, zu, zv, z_o, zl, zva, zua, z_oa
157	REAL, DIMENSION(ngrid, nlay) :: ztemp_env ! temperarure liquide de l'environnement
158	REAL, DIMENSION(ngrid, nlay) :: zta, zha, q2, wq, wthl, wthv, thetath2, wth2
159	REAL, DIMENSION(ngrid, nlay) :: rho, masse
160	REAL, DIMENSION(ngrid, nlay + 1) :: zw_est, zlev
161	REAL, DIMENSION(ngrid) :: wmax, wmax_tmp
162	REAL, DIMENSION(ngrid, nlay + 1) :: f_star
163	REAL, DIMENSION(ngrid, nlay) :: entr, detr, entr_star, detr_star, alim_star_clos
164	REAL, DIMENSION(ngrid, nlay) :: zqsat, csc
165	REAL, DIMENSION(ngrid) :: zcon, zcon2, alim_star_tot, f
166	REAL, DIMENSION(ngrid, nlay) :: entrdn, detrdn
167	logical, DIMENSION(ngrid, nlay) :: mask
168
169	CHARACTER (LEN = 20) :: modname = 'thermcell_main'
170	CHARACTER (LEN = 80) :: abort_message
171
172
173	#ifdef ISO
174	REAL xtpo(ntiso,ngrid,nlay),xtpdoadj(ntiso,ngrid,nlay)
175	REAL xtzo(ntiso,ngrid,nlay)
176	REAL xtpdoadj_tmp(ngrid,nlay)
177	REAL xtpo_tmp(ngrid,nlay)
178	REAL xtzo_tmp(ngrid,nlay)
179	INTEGER ixt
180	#endif
181
182	!-----------------------------------------------------------------------
183	! initialisation:
184	! ---------------
185
186	fm = 0. ; entr = 0. ; detr = 0.
187
188	IF (prt_level>=1) PRINT*, 'thermcell_main V4'
189
190	sorties = .TRUE.
191	IF(ngrid/=ngrid) THEN
192	PRINT*
193	PRINT*, 'STOP dans convadj'
194	PRINT*, 'ngrid =', ngrid
195	PRINT*, 'ngrid =', ngrid
196	ENDIF
197
198	!PRINT*,'thermcell_main debut'
199	! WRITE(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
200	do ig = 1, ngrid
201	f0(ig) = max(f0(ig), 1.e-2)
202	zmax0(ig) = max(zmax0(ig), 40.)
203	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
204	enddo
205
206	IF (prt_level>=20) THEN
207	do ig = 1, ngrid
208	PRINT*, 'th_main ig f0', ig, f0(ig)
209	enddo
210	endif
211
212	!-----------------------------------------------------------------------
213	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
214	! --------------------------------------------------------------------
215
216	! On condense l'eau liquide si besoin.
217	! En fait on arrive ici d'habitude (jusque 6A) après réévaporation
218	! Dans une nouvelle mouture, on passe les profiles
219	! avant la couche limite : iflag_thermals_tenv=1
220	! dés le début de la physique : iflag_thermals_tenv=2
221	! Mais même pour 2) on ne veut sans doute pas réévaporer.
222	! On veut comparer thetav dans le thermique, après condensation,
223	! avec le theta_v effectif de l'environnement.
224
225	IF (iflag_thermals_tenv - 10 * (iflag_thermals_tenv / 10) == 0) THEN
226	CALL thermcell_env(ngrid, nlay, p_o, ptemp_env, puwind, pvwind, pplay, &
227	pplev, z_o, ztemp_env, zl, ztv, zthl, zu, zv, zpspsk, zqsat, lcong, lintercong, lev_out)
228
229	else
230
231	! Chantier en cours : ne pas effacer (Fredho). 15 septembre 2023
232	! Dans la version originale de thermcell_env, on condense l'eau de l'environnement
233	! pour calculer une temperature potentielle liquide.
234	! On en déduit un Theta v.
235
236	! ...
237	! contenu de thermcell_env
238	! SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
239	! & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
240	! contenu thermcell_env : CALL thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
241	! contenu thermcell_env : do ll=1,nlay
242	! contenu thermcell_env : do ig=1,ngrid
243	! contenu thermcell_env : zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
244	! contenu thermcell_env : zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll) ! T = Tl + Lv/Cp ql
245	! contenu thermcell_env : zo(ig,ll) = po(ig,ll)-zl(ig,ll)
246	! contenu thermcell_env : enddo
247	! contenu thermcell_env : enddo
248	! contenu thermcell_env : do ll=1,nlay
249	! contenu thermcell_env : do ig=1,ngrid
250	! contenu thermcell_env : zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
251	! contenu thermcell_env : zu(ig,ll)=pu(ig,ll)
252	! contenu thermcell_env : zv(ig,ll)=pv(ig,ll)
253	! contenu thermcell_env : ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll)
254	! contenu thermcell_env : ztv(ig,ll)=ztv(ig,ll)(1.+RETV(zo(ig,ll))-zl(ig,ll))
255	! contenu thermcell_env : zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll)
256	! contenu thermcell_env : enddo
257	! contenu thermcell_env : enddo
258
259	do l = 1, nlay
260	do ig = 1, ngrid
261	zl(ig, l) = 0.
262	zu(ig, l) = puwind(ig, l)
263	zv(ig, l) = pvwind(ig, l)
264	ztemp_env(ig, l) = ptemp_env(ig, l)
265	zpspsk(ig, l) = (pplay(ig, l) / 100000.)**RKAPPA
266	ztv(ig, l) = ztemp_env(ig, l) / zpspsk(ig, l)
267	ztv(ig, l) = ztv(ig, l) * (1. + RETV * po_env(ig, l))
268	zthl(ig, l) = ptemp(ig, l) / zpspsk(ig, l)
269	mask(ig, l) = .TRUE.
270	enddo
271	enddo
272	CALL thermcell_qsat(ngrid * nlay, mask, pplev, ptemp_env, p_o, zqsat)
273
274	endif
275
276	IF (prt_level>=1) PRINT*, 'thermcell_main apres thermcell_env'
277
278	!------------------------------------------------------------------------
279	! --------------------
280
281
282	! + + + + + + + + + + +
283
284
285	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
286	! wh,wt,wo ...
287
288	! + + + + + + + + + + + zh,zu,zv,z_o,rho
289
290
291	! -------------------- zlev(1)
292	! \\\\\\\\\\\\\\\\\\\\
293
294
295	!-----------------------------------------------------------------------
296	! Calcul des altitudes des couches
297	!-----------------------------------------------------------------------
298
299	do l = 2, nlay
300	zlev(:, l) = 0.5 * (pphi(:, l) + pphi(:, l - 1)) / RG
301	enddo
302	zlev(:, 1) = 0.
303	zlev(:, nlay + 1) = (2. * pphi(:, nlay) - pphi(:, nlay - 1)) / RG
304	do l = 1, nlay
305	zlay(:, l) = pphi(:, l) / RG
306	enddo
307	do l = 1, nlay
308	deltaz(:, l) = zlev(:, l + 1) - zlev(:, l)
309	enddo
310
311	!-----------------------------------------------------------------------
312	! Calcul des densites et masses
313	!-----------------------------------------------------------------------
314
315	rho(:, :) = pplay(:, :) / (zpspsk(:, :) * RD * ztv(:, :))
316	IF (prt_level>=10) WRITE(lunout, *) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
317	rhobarz(:, 1) = rho(:, 1)
318	do l = 2, nlay
319	rhobarz(:, l) = 0.5 * (rho(:, l) + rho(:, l - 1))
320	enddo
321	do l = 1, nlay
322	masse(:, l) = (pplev(:, l) - pplev(:, l + 1)) / RG
323	enddo
324	IF (prt_level>=1) PRINT*, 'thermcell_main apres initialisation'
325
326	!------------------------------------------------------------------
327
328	! /\|\
329	! -------- \| F_k+1 -------
330	! ----> D_k
331	! /\|\ <---- E_k , A_k
332	! -------- \| F_k ---------
333	! ----> D_k-1
334	! <---- E_k-1 , A_k-1
335
336
337
338
339
340	! ---------------------------
341
342	! ----- F_lmax+1=0 ---------- \
343	! lmax (zmax) \|
344	! --------------------------- \|
345	! \|
346	! --------------------------- \|
347	! \|
348	! --------------------------- \|
349	! \|
350	! --------------------------- \|
351	! \|
352	! --------------------------- \|
353	! \| E
354	! --------------------------- \| D
355	! \|
356	! --------------------------- \|
357	! \|
358	! --------------------------- \ \|
359	! lalim \| \|
360	! --------------------------- \| \|
361	! \| \|
362	! --------------------------- \| \|
363	! \| A \|
364	! --------------------------- \| \|
365	! \| \|
366	! --------------------------- \| \|
367	! lmin (=1 pour le moment) \| \|
368	! ----- F_lmin=0 ------------ / /
369
370	! ---------------------------
371	! //////////////////////////
372
373
374	!=============================================================================
375	! Calculs initiaux ne faisant pas intervenir les changements de phase
376	!=============================================================================
377
378	!------------------------------------------------------------------
379	! 1. alim_star est le profil vertical de l'alimentation a la base du
380	! panache thermique, calcule a partir de la flotabilite de l'air sec
381	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
382	!------------------------------------------------------------------
383
384	entr_star = 0. ; detr_star = 0. ; alim_star = 0. ; alim_star_tot = 0.
385	lmin = 1
386
387	!-----------------------------------------------------------------------------
388	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
389	! panache sec conservatif (e=d=0) alimente selon alim_star
390	! Il s'agit d'un calcul de type CAPE
391	! zmax_sec est utilise pour determiner la geometrie du thermique.
392	!------------------------------------------------------------------------------
393	!---------------------------------------------------------------------------------
394	!calcul du melange et des variables dans le thermique
395	!--------------------------------------------------------------------------------
396
397	IF (prt_level>=1) PRINT*, 'avant thermcell_plume ', lev_out
398
399	!=====================================================================
400	! Old version of thermcell_plume in thermcell_plume_6A.F90
401	! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
402	! to the 5B and 6A versions used for CMIP5 and CMIP6.
403	! The latest was previously named thermcellV1_plume.
404	! The new thermcell_plume is a clean version (removing obsolete
405	! options) of thermcell_plume_6A.
406	! The 3 versions are controled by
407	! flag_thermals_ed <= 9 thermcell_plume_6A
408	! <= 19 thermcell_plume_5B
409	! else thermcell_plume (default 20 for convergence with 6A)
410	! Fredho
411	!=====================================================================
412
413	IF (iflag_thermals_ed<=9) THEN
414	! PRINT*,'THERM NOUVELLE/NOUVELLE Arnaud'
415	CALL thermcell_plume_6A(itap, ngrid, nlay, ptimestep, ztv, zthl, p_o, zl, rhobarz, &
416	zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
417	lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
418	ztla, zqla, zqta, zha, zw2, zw_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
419	, lev_out, lunout1, igout)
420
421	elseif (iflag_thermals_ed<=19) THEN
422	! PRINT*,'THERM RIO et al 2010, version d Arnaud'
423	CALL thermcell_plume_5B(itap, ngrid, nlay, ptimestep, ztv, zthl, p_o, zl, rhobarz, &
424	zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
425	lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
426	ztla, zqla, zqta, zha, zw2, zw_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
427	, lev_out, lunout1, igout)
428	else
429	CALL thermcell_plume(itap, ngrid, nlay, ptimestep, ztv, zthl, p_o, zl, rhobarz, &
430	zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
431	lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
432	ztla, zqla, zqta, zha, zw2, zw_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
433	, lev_out, lunout1, igout)
434	endif
435
436	IF (prt_level>=1) PRINT*, 'apres thermcell_plume ', lev_out
437
438	CALL test_ltherm(ngrid, nlay, pplay, lalim, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_plum lalim ')
439	CALL test_ltherm(ngrid, nlay, pplay, lmix, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_plum lmix ')
440
441	IF (prt_level>=1) PRINT*, 'thermcell_main apres thermcell_plume'
442	IF (prt_level>=10) THEN
443	WRITE(lunout1, *) 'Dans thermcell_main 2'
444	WRITE(lunout1, *) 'lmin ', lmin(igout)
445	WRITE(lunout1, *) 'lalim ', lalim(igout)
446	WRITE(lunout1, *) ' ig l alim_star entr_star detr_star f_star '
447	WRITE(lunout1, '(i6,i4,4e15.5)') (igout, l, alim_star(igout, l), entr_star(igout, l), detr_star(igout, l) &
448	, f_star(igout, l + 1), l = 1, nint(linter(igout)) + 5)
449	endif
450
451	!-------------------------------------------------------------------------------
452	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
453	!-------------------------------------------------------------------------------
454
455	CALL thermcell_height(ngrid, nlay, lalim, lmin, linter, lcong, lintercong, lmix, zw2, &
456	zlev, lmax, zmax, zmax0, zmix, wmax, zcong)
457	! Attention, w2 est transforme en sa racine carree dans cette routine
458	! Le probleme vient du fait que linter et lmix sont souvent egaux a 1.
459	wmax_tmp = 0.
460	do l = 1, nlay
461	wmax_tmp(:) = max(wmax_tmp(:), zw2(:, l))
462	enddo
463	! PRINT*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
464
465	CALL test_ltherm(ngrid, nlay, pplay, lalim, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_heig lalim ')
466	CALL test_ltherm(ngrid, nlay, pplay, lmin, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_heig lmin ')
467	CALL test_ltherm(ngrid, nlay, pplay, lmix, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_heig lmix ')
468	CALL test_ltherm(ngrid, nlay, pplay, lmax, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_heig lmax ')
469
470	IF (prt_level>=1) PRINT*, 'thermcell_main apres thermcell_height'
471
472	!-------------------------------------------------------------------------------
473	! Fermeture,determination de f
474	!-------------------------------------------------------------------------------
475
476	CALL thermcell_dry(ngrid, nlay, zlev, pphi, ztv, alim_star, &
477	lalim, lmin, zmax_sec, wmax_sec)
478
479	CALL test_ltherm(ngrid, nlay, pplay, lmin, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_dry lmin ')
480	CALL test_ltherm(ngrid, nlay, pplay, lalim, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_dry lalim ')
481
482	IF (prt_level>=1) PRINT*, 'thermcell_main apres thermcell_dry'
483	IF (prt_level>=10) THEN
484	WRITE(lunout1, *) 'Dans thermcell_main 1b'
485	WRITE(lunout1, *) 'lmin ', lmin(igout)
486	WRITE(lunout1, *) 'lalim ', lalim(igout)
487	WRITE(lunout1, *) ' ig l alim_star entr_star detr_star f_star '
488	WRITE(lunout1, '(i6,i4,e15.5)') (igout, l, alim_star(igout, l) &
489	, l = 1, lalim(igout) + 4)
490	endif
491
492
493
494
495	! Choix de la fonction d'alimentation utilisee pour la fermeture.
496	! Apparemment sans importance
497	alim_star_clos(:, :) = alim_star(:, :)
498	alim_star_clos(:, :) = entr_star(:, :) + alim_star(:, :)
499
500	!CR Appel de la fermeture seche
501	IF (iflag_thermals_closure==1) THEN
502	CALL thermcell_closure(ngrid, nlay, r_aspect_thermals, ptimestep, rho, &
503	zlev, lalim, alim_star_clos, zmax_sec, wmax_sec, f)
504
505	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
506	! Appel avec les zmax et wmax tenant compte de la condensation
507	! Semble moins bien marcher
508	ELSE IF (iflag_thermals_closure==2) THEN
509	CALL thermcell_closure(ngrid, nlay, r_aspect_thermals, ptimestep, rho, &
510	zlev, lalim, alim_star, zmax, wmax, f)
511
512	endif
513
514	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
515
516	IF(prt_level>=1)PRINT*, 'thermcell_closure apres thermcell_closure'
517
518	IF (tau_thermals>1.) THEN
519	lambda = exp(-ptimestep / tau_thermals)
520	f0 = (1. - lambda) * f + lambda * f0
521	else
522	f0 = f
523	endif
524
525	! Test valable seulement en 1D mais pas genant
526	IF (.NOT. (f0(1)>=0.)) THEN
527	abort_message = '.NOT. (f0(1).ge.0.)'
528	CALL abort_physic (modname, abort_message, 1)
529	endif
530
531	!-------------------------------------------------------------------------------
532	!deduction des flux
533
534	CALL thermcell_flux2(ngrid, nlay, ptimestep, masse, &
535	lalim, lmax, alim_star, &
536	entr_star, detr_star, f, rhobarz, zlev, zw2, fm, entr, &
537	detr, zqla, lev_out, lunout1, igout)
538
539	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
540
541	IF (prt_level>=1) PRINT*, 'thermcell_main apres thermcell_flux'
542	CALL test_ltherm(ngrid, nlay, pplay, lalim, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_flux lalim ')
543	CALL test_ltherm(ngrid, nlay, pplay, lmax, ztv, p_o, ztva, zqla, f_star, zw2, 'thermcell_flux lmax ')
544
545	!------------------------------------------------------------------
546	! On ne prend pas directement les profils issus des calculs precedents
547	! mais on s'autorise genereusement une relaxation vers ceci avec
548	! une constante de temps tau_thermals (typiquement 1800s).
549	!------------------------------------------------------------------
550
551	IF (tau_thermals>1.) THEN
552	lambda = exp(-ptimestep / tau_thermals)
553	fm0 = (1. - lambda) * fm + lambda * fm0
554	entr0 = (1. - lambda) * entr + lambda * entr0
555	detr0 = (1. - lambda) * detr + lambda * detr0
556	else
557	fm0 = fm
558	entr0 = entr
559	detr0 = detr
560	endif
561
562	!------------------------------------------------------------------
563	! Calcul de la fraction de l'ascendance
564	!------------------------------------------------------------------
565	do ig = 1, ngrid
566	fraca(ig, 1) = 0.
567	fraca(ig, nlay + 1) = 0.
568	enddo
569	do l = 2, nlay
570	do ig = 1, ngrid
571	IF (zw2(ig, l)>1.e-10) THEN
572	fraca(ig, l) = fm(ig, l) / (rhobarz(ig, l) * zw2(ig, l))
573	else
574	fraca(ig, l) = 0.
575	endif
576	enddo
577	enddo
578
579	!c------------------------------------------------------------------
580	! calcul du transport vertical
581	!------------------------------------------------------------------
582	IF (iflag_thermals_down > 0) THEN
583	IF (debut) PRINT*, 'WARNING !!! routine thermcell_down en cours de developpement'
584	entrdn = fact_thermals_down * detr0
585	detrdn = fact_thermals_down * entr0
586	! we want to transport potential temperature, total water and momentum
587	CALL thermcell_updown_dq(ngrid, nlay, ptimestep, lmax, entr0, detr0, entrdn, detrdn, masse, zthl, zdthladj)
588	CALL thermcell_updown_dq(ngrid, nlay, ptimestep, lmax, entr0, detr0, entrdn, detrdn, masse, p_o, pdoadj)
589	CALL thermcell_updown_dq(ngrid, nlay, ptimestep, lmax, entr0, detr0, entrdn, detrdn, masse, zu, pduadj)
590	CALL thermcell_updown_dq(ngrid, nlay, ptimestep, lmax, entr0, detr0, entrdn, detrdn, masse, zv, pdvadj)
591	ELSE
592	!--------------------------------------------------------------
593
594	! Temperature potentielle liquide effectivement mélangée par les thermiques
595	do ll = 1, nlay
596	do ig = 1, ngrid
597	zthl(ig, ll) = ptemp(ig, ll) / zpspsk(ig, ll)
598	enddo
599	enddo
600	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, &
601	zthl, zdthladj, zta, lev_out)
602
603	do ll = 1, nlay
604	do ig = 1, ngrid
605	z_o(ig, ll) = p_o(ig, ll)
606	enddo
607	enddo
608	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, &
609	z_o, pdoadj, z_oa, lev_out)
610
611	#ifdef ISO
612	! C Risi: on utilise directement la meme routine
613	do ixt=1,ntiso
614	do ll=1,nlay
615	DO ig=1,ngrid
616	xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
617	xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
618	enddo
619	enddo
620	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
621	xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
622	do ll=1,nlay
623	DO ig=1,ngrid
624	xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
625	enddo
626	enddo
627	enddo
628	#endif
629
630	#ifdef ISO
631	#ifdef ISOVERIF
632	DO ll=1,nlay
633	DO ig=1,ngrid
634	IF (iso_eau.gt.0) THEN
635	CALL iso_verif_egalite(xtpo(iso_eau,ig,ll), &
636	p_o(ig,ll),'thermcell_main 594')
637	CALL iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
638	pdoadj(ig,ll),'thermcell_main 596')
639	endif
640	IF (iso_HDO.gt.0) THEN
641	CALL iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
642	/p_o(ig,ll),'thermcell_main 610')
643	endif
644	enddo
645	enddo !DO ll=1,nlay
646	WRITE(,) 'thermcell_main 600 tmp: apres thermcell_dq'
647	#endif
648	#endif
649
650
651	!------------------------------------------------------------------
652	! calcul du transport vertical du moment horizontal
653	!------------------------------------------------------------------
654
655	!IM 090508
656	IF (dvdq == 0) THEN
657	! Calcul du transport de V tenant compte d'echange par gradient
658	! de pression horizontal avec l'environnement
659
660	CALL thermcell_dv2(ngrid, nlay, ptimestep, fm0, entr0, masse &
661	! & ,fraca*dvdq,zmax &
662	, fraca, zmax &
663	, zu, zv, pduadj, pdvadj, zua, zva, lev_out)
664
665	else
666
667	! calcul purement conservatif pour le transport de V
668	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse &
669	, zu, pduadj, zua, lev_out)
670	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse &
671	, zv, pdvadj, zva, lev_out)
672
673	endif
674	ENDIF
675
676	! PRINT*,'13 OK convect8'
677	do l = 1, nlay
678	do ig = 1, ngrid
679	pdtadj(ig, l) = zdthladj(ig, l) * zpspsk(ig, l)
680	enddo
681	enddo
682
683	IF (prt_level>=1) PRINT*, '14 OK convect8'
684	!------------------------------------------------------------------
685	! Calculs de diagnostiques pour les sorties
686	!------------------------------------------------------------------
687	!calcul de fraca pour les sorties
688
689	IF (sorties) THEN
690	IF (prt_level>=1) PRINT*, '14a OK convect8'
691	! calcul du niveau de condensation
692	! initialisation
693	do ig = 1, ngrid
694	nivcon(ig) = 0
695	zcon(ig) = 0.
696	enddo
697	!nouveau calcul
698	do ig = 1, ngrid
699	! WARNING !!! verifier que c'est bien ztemp_env qu'on veut là
700	CHI = ztemp_env(ig, 1) / (1669.0 - 122.0 * z_o(ig, 1) / zqsat(ig, 1) - ztemp_env(ig, 1))
701	pcon(ig) = pplay(ig, 1) * (z_o(ig, 1) / zqsat(ig, 1))**CHI
702	enddo
703	!IM do k=1,nlay
704	do k = 1, nlay - 1
705	do ig = 1, ngrid
706	IF ((pcon(ig)<=pplay(ig, k)) &
707	.AND.(pcon(ig)>pplay(ig, k + 1))) THEN
708	zcon2(ig) = zlay(ig, k) - (pcon(ig) - pplay(ig, k)) / (RG * rho(ig, k)) / 100.
709	endif
710	enddo
711	enddo
712	!IM
713	ierr = 0
714	do ig = 1, ngrid
715	IF (pcon(ig)<=pplay(ig, nlay)) THEN
716	zcon2(ig) = zlay(ig, nlay) - (pcon(ig) - pplay(ig, nlay)) / (RG * rho(ig, nlay)) / 100.
717	ierr = 1
718	endif
719	enddo
720	! if (ierr==1) THEN
721	! abort_message = 'thermcellV0_main: les thermiques vont trop haut '
722	! CALL abort_physic (modname,abort_message,1)
723	! endif
724
725	IF (prt_level>=1) PRINT*, '14b OK convect8'
726	do k = nlay, 1, -1
727	do ig = 1, ngrid
728	IF (zqla(ig, k)>1e-10) THEN
729	nivcon(ig) = k
730	zcon(ig) = zlev(ig, k)
731	endif
732	enddo
733	enddo
734	IF (prt_level>=1) PRINT*, '14c OK convect8'
735	!calcul des moments
736	!initialisation
737	do l = 1, nlay
738	do ig = 1, ngrid
739	q2(ig, l) = 0.
740	wth2(ig, l) = 0.
741	wth3(ig, l) = 0.
742	ratqscth(ig, l) = 0.
743	ratqsdiff(ig, l) = 0.
744	enddo
745	enddo
746	IF (prt_level>=1) PRINT*, '14d OK convect8'
747	IF (prt_level>=10)WRITE(lunout, *) &
748	'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
749	do l = 1, nlay
750	do ig = 1, ngrid
751	zf = fraca(ig, l)
752	zf2 = zf / (1. - zf)
753
754	thetath2(ig, l) = zf2 * (ztla(ig, l) - zthl(ig, l))**2
755	IF(zw2(ig, l)>1.e-10) THEN
756	wth2(ig, l) = zf2 * (zw2(ig, l))**2
757	else
758	wth2(ig, l) = 0.
759	endif
760	wth3(ig, l) = zf2 * (1 - 2. * fraca(ig, l)) / (1 - fraca(ig, l)) &
761	* zw2(ig, l) * zw2(ig, l) * zw2(ig, l)
762	q2(ig, l) = zf2 * (zqta(ig, l) * 1000. - p_o(ig, l) * 1000.)**2
763	!test: on calcul q2/p_o=ratqsc
764	ratqscth(ig, l) = sqrt(max(q2(ig, l), 1.e-6) / (p_o(ig, l) * 1000.))
765	enddo
766	enddo
767	!calcul des flux: q, thetal et thetav
768	do l = 1, nlay
769	do ig = 1, ngrid
770	wq(ig, l) = fraca(ig, l) * zw2(ig, l) * (zqta(ig, l) * 1000. - p_o(ig, l) * 1000.)
771	wthl(ig, l) = fraca(ig, l) * zw2(ig, l) * (ztla(ig, l) - zthl(ig, l))
772	wthv(ig, l) = fraca(ig, l) * zw2(ig, l) * (ztva(ig, l) - ztv(ig, l))
773	enddo
774	enddo
775
776	!calcul du ratqscdiff
777	IF (prt_level>=1) PRINT*, '14e OK convect8'
778	var = 0.
779	vardiff = 0.
780	ratqsdiff(:, :) = 0.
781
782	do l = 1, nlay
783	do ig = 1, ngrid
784	IF (l<=lalim(ig)) THEN
785	var = var + alim_star(ig, l) * zqta(ig, l) * 1000.
786	endif
787	enddo
788	enddo
789
790	IF (prt_level>=1) PRINT*, '14f OK convect8'
791
792	do l = 1, nlay
793	do ig = 1, ngrid
794	IF (l<=lalim(ig)) THEN
795	zf = fraca(ig, l)
796	zf2 = zf / (1. - zf)
797	vardiff = vardiff + alim_star(ig, l) * (zqta(ig, l) * 1000. - var)**2
798	endif
799	enddo
800	enddo
801
802	IF (prt_level>=1) PRINT*, '14g OK convect8'
803	do l = 1, nlay
804	do ig = 1, ngrid
805	ratqsdiff(ig, l) = sqrt(vardiff) / (p_o(ig, l) * 1000.)
806	enddo
807	enddo
808	endif
809
810	IF (prt_level>=1) PRINT*, 'thermcell_main FIN OK'
811
812	!PRINT*,'thermcell_main fin'
813
814	END SUBROUTINE thermcell_main
815
816	!=============================================================================
817	!/////////////////////////////////////////////////////////////////////////////
818	!=============================================================================
819	SUBROUTINE test_ltherm(ngrid, nlay, pplay, long, ztv, p_o, ztva, & ! in
820	zqla, f_star, zw2, comment) ! in
821	!=============================================================================
822	USE lmdz_thermcell_ini, ONLY: prt_level
823	IMPLICIT NONE
824
825	INTEGER i, k, ngrid, nlay
826	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: pplay, ztv, p_o, ztva, zqla
827	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: f_star, zw2
828	INTEGER, INTENT(IN), DIMENSION(ngrid) :: long
829	REAL seuil
830	character*21 comment
831
832	seuil = 0.25
833
834	IF (prt_level>=1) THEN
835	PRINT*, 'WARNING !!! TEST ', comment
836	endif
837	return
838
839	! test sur la hauteur des thermiques ...
840	do i = 1, ngrid
841	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) THEN
842	IF (prt_level>=10) THEN
843	PRINT*, 'WARNING ', comment, ' au point ', i, ' K= ', long(i)
844	PRINT, ' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
845	do k = 1, nlay
846	WRITE(6, '(i3,7f10.3)') k, pplay(i, k), ztv(i, k), 1000 * p_o(i, k), ztva(i, k), 1000 * zqla(i, k), f_star(i, k), zw2(i, k)
847	enddo
848	endif
849	enddo
850
851	RETURN
852	end
853
854	! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
855	! On transporte pbl_tke pour donner therm_tke
856	! Copie conforme de la SUBROUTINE DTKE dans physiq.F ecrite par Frederic Hourdin
857
858	!=======================================================================
859	!///////////////////////////////////////////////////////////////////////
860	!=======================================================================
861
862	SUBROUTINE thermcell_tke_transport(&
863	ngrid, nlay, ptimestep, fm0, entr0, rg, pplev, & ! in
864	therm_tke_max) ! out
865	USE lmdz_thermcell_ini, ONLY: prt_level
866	IMPLICIT NONE
867
868	!=======================================================================
869
870	! Calcul du transport verticale dans la couche limite en presence
871	! de "thermiques" explicitement representes
872	! calcul du dq/dt une fois qu'on connait les ascendances
873
874	!=======================================================================
875
876	INTEGER ngrid, nlay
877
878	REAL, INTENT(IN) :: ptimestep
879	REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fm0, pplev
880	REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: entr0
881	REAL, INTENT(IN) :: rg
882	REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: therm_tke_max
883
884	REAL detr0(ngrid, nlay)
885	REAL masse0(ngrid, nlay)
886	REAL masse(ngrid, nlay), fm(ngrid, nlay + 1)
887	REAL entr(ngrid, nlay)
888	REAL q(ngrid, nlay)
889	INTEGER lev_out ! niveau pour les print
890
891	REAL qa(ngrid, nlay), detr(ngrid, nlay), wqd(ngrid, nlay + 1)
892	INTEGER ig, k
893
894	lev_out = 0
895
896	IF (prt_level>=1) PRINT*, 'Q2 THERMCEL_DQ 0'
897
898	! calcul du detrainement
899	do k = 1, nlay
900	detr0(:, k) = fm0(:, k) - fm0(:, k + 1) + entr0(:, k)
901	masse0(:, k) = (pplev(:, k) - pplev(:, k + 1)) / RG
902	enddo
903
904
905	! Decalage vertical des entrainements et detrainements.
906	masse(:, 1) = 0.5 * masse0(:, 1)
907	entr(:, 1) = 0.5 * entr0(:, 1)
908	detr(:, 1) = 0.5 * detr0(:, 1)
909	fm(:, 1) = 0.
910	do k = 1, nlay - 1
911	masse(:, k + 1) = 0.5 * (masse0(:, k) + masse0(:, k + 1))
912	entr(:, k + 1) = 0.5 * (entr0(:, k) + entr0(:, k + 1))
913	detr(:, k + 1) = 0.5 * (detr0(:, k) + detr0(:, k + 1))
914	fm(:, k + 1) = fm(:, k) + entr(:, k) - detr(:, k)
915	enddo
916	fm(:, nlay + 1) = 0.
917
918	q(:, :) = therm_tke_max(:, :)
919	!!! nrlmd le 16/09/2010
920	do ig = 1, ngrid
921	qa(ig, 1) = q(ig, 1)
922	enddo
923	!!!
924
925	IF (1==1) THEN
926	do k = 2, nlay
927	do ig = 1, ngrid
928	IF ((fm(ig, k + 1) + detr(ig, k)) * ptimestep> &
929	1.e-5 * masse(ig, k)) THEN
930	qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k)) &
931	/ (fm(ig, k + 1) + detr(ig, k))
932	else
933	qa(ig, k) = q(ig, k)
934	endif
935	IF (qa(ig, k)<0.) THEN
936	! PRINT*,'qa<0!!!'
937	endif
938	IF (q(ig, k)<0.) THEN
939	! PRINT*,'q<0!!!'
940	endif
941	enddo
942	enddo
943
944	! Calcul du flux subsident
945
946	do k = 2, nlay
947	do ig = 1, ngrid
948	wqd(ig, k) = fm(ig, k) * q(ig, k)
949	IF (wqd(ig, k)<0.) THEN
950	! PRINT*,'wqd<0!!!'
951	endif
952	enddo
953	enddo
954	do ig = 1, ngrid
955	wqd(ig, 1) = 0.
956	wqd(ig, nlay + 1) = 0.
957	enddo
958
959	! Calcul des tendances
960	do k = 1, nlay
961	do ig = 1, ngrid
962	q(ig, k) = q(ig, k) + (detr(ig, k) * qa(ig, k) - entr(ig, k) * q(ig, k) &
963	- wqd(ig, k) + wqd(ig, k + 1)) &
964	* ptimestep / masse(ig, k)
965	enddo
966	enddo
967
968	endif
969
970	therm_tke_max(:, :) = q(:, :)
971
972	RETURN
973	!!! fin nrlmd le 10/04/2012
974	end
975
976	END MODULE lmdz_thermcell_main

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

Download in other formats:

Original Format