Context Navigation

lmdz_thermcell_main.F90 @ 5877

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

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

Download in other formats:

Original Format