Context Navigation

lmdz_thermcell_main.F90 @ 5878

Last change on this file since 5878 was 5853, checked in by fhourdin, 6 weeks ago
Separation thermcell_plume 5B/6A
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: 35.2 KB

Line
1	MODULE lmdz_thermcell_main
2	! $Id: lmdz_thermcell_main.F90 5853 2025-10-10 13:59:07Z ymeurdesoif $
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
42	USE lmdz_thermcell_plume_5B, ONLY: 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
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
143	integer, parameter :: igout=1
144	integer, parameter :: lunout1=6
145	integer, parameter :: lev_out=10
146
147	real lambda, zf,zf2,var,vardiff,CHI
148	integer ig,k,l,ierr,ll
149	logical sorties
150	real, dimension(ngrid) :: linter,zmix, zmax_sec,lintercong
151	integer,dimension(ngrid) :: lmin,lmix,lmix_bis,nivcon, lcong
152	real, dimension(ngrid,nlay) :: ztva_est
153	real, dimension(ngrid,nlay) :: deltaz,zlay,zdthladj,zu,zv,z_o,zl,zva,zua,z_oa
154	real, dimension(ngrid,nlay) :: ztemp_env ! temperarure liquide de l'environnement
155	real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2
156	real, dimension(ngrid,nlay) :: rho,masse
157	real, dimension(ngrid,nlay+1) :: zw_est,zlev
158	real, dimension(ngrid) :: wmax,wmax_tmp
159	real, dimension(ngrid,nlay+1) :: f_star
160	real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos
161	real, dimension(ngrid,nlay) :: zqsat,csc
162	real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f
163	real, dimension(ngrid,nlay) :: entrdn,detrdn
164	logical, dimension(ngrid,nlay) :: mask
165
166	character (len=20), parameter :: modname='thermcell_main'
167	character (len=80) :: abort_message
168
169
170	#ifdef ISO
171	REAL xtpo(ntiso,ngrid,nlay),xtpdoadj(ntiso,ngrid,nlay)
172	REAL xtzo(ntiso,ngrid,nlay)
173	REAL xtpdoadj_tmp(ngrid,nlay)
174	REAL xtpo_tmp(ngrid,nlay)
175	REAL xtzo_tmp(ngrid,nlay)
176	integer ixt
177	#endif
178
179	!
180
181	!-----------------------------------------------------------------------
182	! initialisation:
183	! ---------------
184	!
185	fm=0. ; entr=0. ; detr=0.
186
187	if (prt_level.ge.1) print*,'thermcell_main V4'
188
189	sorties=.true.
190	IF(ngrid.NE.ngrid) THEN
191	PRINT*,'STOP dans convadj'
192	PRINT*,'ngrid =',ngrid
193	PRINT*,'ngrid =',ngrid
194	ENDIF
195	!
196	!print*,'thermcell_main debut'
197	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
198	do ig=1,ngrid
199	f0(ig)=max(f0(ig),1.e-2)
200	zmax0(ig)=max(zmax0(ig),40.)
201	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
202	enddo
203
204	if (prt_level.ge.20) then
205	do ig=1,ngrid
206	print*,'th_main ig f0',ig,f0(ig)
207	enddo
208	endif
209
210	!-----------------------------------------------------------------------
211	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
212	! --------------------------------------------------------------------
213	!
214	! On condense l'eau liquide si besoin.
215	! En fait on arrive ici d'habitude (jusque 6A) après réévaporation
216	! Dans une nouvelle mouture, on passe les profiles
217	! avant la couche limite : iflag_thermals_tenv=1
218	! dés le début de la physique : iflag_thermals_tenv=2
219	! Mais même pour 2) on ne veut sans doute pas réévaporer.
220	! On veut comparer thetav dans le thermique, après condensation,
221	! avec le theta_v effectif de l'environnement.
222
223	if (iflag_thermals_tenv - 10 * ( iflag_thermals_tenv / 10 ) == 0) then
224
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.+RETVpo_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.ge.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	!-----------------------------------------------------------------------
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(:,:)RDztv(:,:))
316	if (prt_level.ge.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.ge.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.ge.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.ge.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.ge.1) print*,'thermcell_main apres thermcell_plume'
442	if (prt_level.ge.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
466
467	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
468	call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
469	call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
470	call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
471
472	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
473
474	!-------------------------------------------------------------------------------
475	! Fermeture,determination de f
476	!-------------------------------------------------------------------------------
477	!
478	!
479	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
480	& lalim,lmin,zmax_sec,wmax_sec)
481
482
483	call test_ltherm(ngrid,nlay,pplay,lmin,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
484	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
485
486	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
487	if (prt_level.ge.10) then
488	write(lunout1,*) 'Dans thermcell_main 1b'
489	write(lunout1,*) 'lmin ',lmin(igout)
490	write(lunout1,*) 'lalim ',lalim(igout)
491	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
492	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
493	& ,l=1,lalim(igout)+4)
494	endif
495
496
497
498
499	! Choix de la fonction d'alimentation utilisee pour la fermeture.
500	! Apparemment sans importance
501	alim_star_clos(:,:)=alim_star(:,:)
502	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
503	!
504	!CR Appel de la fermeture seche
505	if (iflag_thermals_closure.eq.1) then
506
507	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
508	& zlev,lalim,alim_star_clos,zmax_sec,wmax_sec,f)
509
510	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
511	! Appel avec les zmax et wmax tenant compte de la condensation
512	! Semble moins bien marcher
513	else if (iflag_thermals_closure.eq.2) then
514
515	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
516	& zlev,lalim,alim_star,zmax,wmax,f)
517
518
519	endif
520
521	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522
523	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
524
525	if (tau_thermals>1.) then
526	lambda=exp(-ptimestep/tau_thermals)
527	f0=(1.-lambda)f+lambdaf0
528	else
529	f0=f
530	endif
531
532	! Test valable seulement en 1D mais pas genant
533	if (.not. (f0(1).ge.0.) ) then
534	abort_message = '.not. (f0(1).ge.0.)'
535	CALL abort_physic (modname,abort_message,1)
536	endif
537
538	!-------------------------------------------------------------------------------
539	!deduction des flux
540
541	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
542	& lalim,lmax,alim_star, &
543	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
544	& detr,zqla,lev_out,lunout1,igout)
545
546	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
547
548	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
549	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
550	call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
551
552	!------------------------------------------------------------------
553	! On ne prend pas directement les profils issus des calculs precedents
554	! mais on s'autorise genereusement une relaxation vers ceci avec
555	! une constante de temps tau_thermals (typiquement 1800s).
556	!------------------------------------------------------------------
557
558	if (tau_thermals>1.) then
559	lambda=exp(-ptimestep/tau_thermals)
560	fm0=(1.-lambda)fm+lambdafm0
561	entr0=(1.-lambda)entr+lambdaentr0
562	detr0=(1.-lambda)detr+lambdadetr0
563	else
564	fm0=fm
565	entr0=entr
566	detr0=detr
567	endif
568
569	!------------------------------------------------------------------
570	! Calcul de la fraction de l'ascendance
571	!------------------------------------------------------------------
572	do ig=1,ngrid
573	fraca(ig,1)=0.
574	fraca(ig,nlay+1)=0.
575	enddo
576	do l=2,nlay
577	do ig=1,ngrid
578	if (zw2(ig,l).gt.1.e-10) then
579	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
580	else
581	fraca(ig,l)=0.
582	endif
583	enddo
584	enddo
585
586	!c------------------------------------------------------------------
587	! calcul du transport vertical
588	!------------------------------------------------------------------
589	IF (iflag_thermals_down .GT. 0) THEN
590	if (debut) print*,'WARNING !!! routine thermcell_down en cours de developpement'
591	entrdn=fact_thermals_down*detr0
592	detrdn=fact_thermals_down*entr0
593	! we want to transport potential temperature, total water and momentum
594	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zthl,zdthladj)
595	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,p_o,pdoadj)
596	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zu,pduadj)
597	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zv,pdvadj)
598	ELSE
599	!--------------------------------------------------------------
600
601	! Temperature potentielle liquide effectivement mélangée par les thermiques
602	do ll=1,nlay
603	do ig=1,ngrid
604	zthl(ig,ll)=ptemp(ig,ll)/zpspsk(ig,ll)
605	enddo
606	enddo
607	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
608	& zthl,zdthladj,zta,lev_out)
609
610	do ll=1,nlay
611	do ig=1,ngrid
612	z_o(ig,ll)=p_o(ig,ll)
613	enddo
614	enddo
615	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
616	& z_o,pdoadj,z_oa,lev_out)
617
618	#ifdef ISO
619	! C Risi: on utilise directement la meme routine
620	do ixt=1,ntiso
621	do ll=1,nlay
622	DO ig=1,ngrid
623	xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
624	xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
625	enddo
626	enddo
627	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
628	& xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
629	do ll=1,nlay
630	DO ig=1,ngrid
631	xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
632	enddo
633	enddo
634	enddo
635	#endif
636
637	#ifdef ISO
638	#ifdef ISOVERIF
639	DO ll=1,nlay
640	DO ig=1,ngrid
641	if (iso_eau.gt.0) then
642	call iso_verif_egalite(xtpo(iso_eau,ig,ll), &
643	& p_o(ig,ll),'thermcell_main 594')
644	call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
645	& pdoadj(ig,ll),'thermcell_main 596')
646	endif
647	if (iso_HDO.gt.0) then
648	call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
649	& /p_o(ig,ll),'thermcell_main 610')
650	endif
651	enddo
652	enddo !DO ll=1,nlay
653	write(,) 'thermcell_main 600 tmp: apres thermcell_dq'
654	#endif
655	#endif
656
657
658	!------------------------------------------------------------------
659	! calcul du transport vertical du moment horizontal
660	!------------------------------------------------------------------
661
662	!IM 090508
663	if (dvdq == 0 ) then
664
665	! Calcul du transport de V tenant compte d'echange par gradient
666	! de pression horizontal avec l'environnement
667
668	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
669	! & ,fraca*dvdq,zmax &
670	& ,fraca,zmax &
671	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
672
673	else
674
675	! calcul purement conservatif pour le transport de V
676	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
677	& ,zu,pduadj,zua,lev_out)
678	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
679	& ,zv,pdvadj,zva,lev_out)
680
681	endif
682	ENDIF
683
684	! print*,'13 OK convect8'
685	do l=1,nlay
686	do ig=1,ngrid
687	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
688	enddo
689	enddo
690
691	if (prt_level.ge.1) print*,'14 OK convect8'
692	!------------------------------------------------------------------
693	! Calculs de diagnostiques pour les sorties
694	!------------------------------------------------------------------
695	!calcul de fraca pour les sorties
696
697	if (sorties) then
698	if (prt_level.ge.1) print*,'14a OK convect8'
699	! calcul du niveau de condensation
700	! initialisation
701	do ig=1,ngrid
702	nivcon(ig)=0
703	zcon(ig)=0.
704	enddo
705	!nouveau calcul
706	do ig=1,ngrid
707	! WARNING !!! verifier que c'est bien ztemp_env qu'on veut là
708	CHI=ztemp_env(ig,1)/(1669.0-122.0*z_o(ig,1)/zqsat(ig,1)-ztemp_env(ig,1))
709	pcon(ig)=pplay(ig,1)(z_o(ig,1)/zqsat(ig,1))*CHI
710	enddo
711	!IM do k=1,nlay
712	do k=1,nlay-1
713	do ig=1,ngrid
714	if ((pcon(ig).le.pplay(ig,k)) &
715	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
716	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
717	endif
718	enddo
719	enddo
720	!IM
721	ierr=0
722	do ig=1,ngrid
723	if (pcon(ig).le.pplay(ig,nlay)) then
724	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
725	ierr=1
726	endif
727	enddo
728	! if (ierr==1) then
729	! abort_message = 'thermcellV0_main: les thermiques vont trop haut '
730	! CALL abort_physic (modname,abort_message,1)
731	! endif
732
733	if (prt_level.ge.1) print*,'14b OK convect8'
734	do k=nlay,1,-1
735	do ig=1,ngrid
736	if (zqla(ig,k).gt.1e-10) then
737	nivcon(ig)=k
738	zcon(ig)=zlev(ig,k)
739	endif
740	enddo
741	enddo
742	if (prt_level.ge.1) print*,'14c OK convect8'
743	!calcul des moments
744	!initialisation
745	do l=1,nlay
746	do ig=1,ngrid
747	q2(ig,l)=0.
748	wth2(ig,l)=0.
749	wth3(ig,l)=0.
750	ratqscth(ig,l)=0.
751	ratqsdiff(ig,l)=0.
752	enddo
753	enddo
754	if (prt_level.ge.1) print*,'14d OK convect8'
755	if (prt_level.ge.10)write(lunout,*) &
756	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
757	do l=1,nlay
758	do ig=1,ngrid
759	zf=fraca(ig,l)
760	zf2=zf/(1.-zf)
761	!
762	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
763	if(zw2(ig,l).gt.1.e-10) then
764	wth2(ig,l)=zf2(zw2(ig,l))*2
765	else
766	wth2(ig,l)=0.
767	endif
768	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
769	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
770	q2(ig,l)=zf2(zqta(ig,l)1000.-p_o(ig,l)1000.)*2
771	!test: on calcul q2/p_o=ratqsc
772	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(p_o(ig,l)*1000.))
773	enddo
774	enddo
775	!calcul des flux: q, thetal et thetav
776	do l=1,nlay
777	do ig=1,ngrid
778	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-p_o(ig,l)1000.)
779	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
780	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
781	enddo
782	enddo
783
784	!calcul du ratqscdiff
785	if (prt_level.ge.1) print*,'14e OK convect8'
786	var=0.
787	vardiff=0.
788	ratqsdiff(:,:)=0.
789
790	do l=1,nlay
791	do ig=1,ngrid
792	if (l<=lalim(ig)) then
793	var=var+alim_star(ig,l)zqta(ig,l)1000.
794	endif
795	enddo
796	enddo
797
798	if (prt_level.ge.1) print*,'14f OK convect8'
799
800	do l=1,nlay
801	do ig=1,ngrid
802	if (l<=lalim(ig)) then
803	zf=fraca(ig,l)
804	zf2=zf/(1.-zf)
805	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
806	endif
807	enddo
808	enddo
809
810	if (prt_level.ge.1) print*,'14g OK convect8'
811	do l=1,nlay
812	do ig=1,ngrid
813	ratqsdiff(ig,l)=sqrt(vardiff)/(p_o(ig,l)*1000.)
814	enddo
815	enddo
816	endif
817
818	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
819
820	!print*,'thermcell_main fin'
821	RETURN
822	end subroutine thermcell_main
823
824	!=============================================================================
825	!/////////////////////////////////////////////////////////////////////////////
826	!=============================================================================
827	subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,p_o,ztva, & ! in
828	& zqla,f_star,zw2,comment) ! in
829	!=============================================================================
830	USE lmdz_thermcell_ini, ONLY: prt_level
831	IMPLICIT NONE
832
833	integer i, k, ngrid,nlay
834	real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,p_o,ztva,zqla
835	real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
836	integer, intent(in), dimension(ngrid) :: long
837	real seuil
838	character*21 comment
839
840	seuil=0.25
841
842	if (prt_level.ge.1) THEN
843	print*,'WARNING !!! TEST ',comment
844	endif
845	return
846
847	! test sur la hauteur des thermiques ...
848	do i=1,ngrid
849	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
850	if (prt_level.ge.10) then
851	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
852	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
853	do k=1,nlay
854	write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000p_o(i,k),ztva(i,k),1000zqla(i,k),f_star(i,k),zw2(i,k)
855	enddo
856	endif
857	enddo
858
859
860	return
861	end subroutine test_ltherm
862
863	! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
864	! On transporte pbl_tke pour donner therm_tke
865	! Copie conforme de la subroutine DTKE dans physiq.F ecrite par Frederic Hourdin
866
867	!=======================================================================
868	!///////////////////////////////////////////////////////////////////////
869	!=======================================================================
870
871	subroutine thermcell_tke_transport( &
872	& ngrid,nlay,ptimestep,fm0,entr0,rg,pplev, & ! in
873	& therm_tke_max) ! out
874	USE lmdz_thermcell_ini, ONLY: prt_level
875	implicit none
876
877	!=======================================================================
878	!
879	! Calcul du transport verticale dans la couche limite en presence
880	! de "thermiques" explicitement representes
881	! calcul du dq/dt une fois qu'on connait les ascendances
882	!
883	!=======================================================================
884
885	integer ngrid,nlay
886
887	real, intent(in) :: ptimestep
888	real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
889	real, intent(in), dimension(ngrid,nlay) :: entr0
890	real, intent(in) :: rg
891	real, intent(out), dimension(ngrid,nlay) :: therm_tke_max
892
893	real detr0(ngrid,nlay)
894	real masse0(ngrid,nlay)
895	real masse(ngrid,nlay),fm(ngrid,nlay+1)
896	real entr(ngrid,nlay)
897	real q(ngrid,nlay)
898	integer lev_out ! niveau pour les print
899
900	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
901	integer ig,k
902
903
904	lev_out=0
905
906
907	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
908
909	! calcul du detrainement
910	do k=1,nlay
911	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
912	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
913	enddo
914
915
916	! Decalage vertical des entrainements et detrainements.
917	masse(:,1)=0.5*masse0(:,1)
918	entr(:,1)=0.5*entr0(:,1)
919	detr(:,1)=0.5*detr0(:,1)
920	fm(:,1)=0.
921	do k=1,nlay-1
922	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
923	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
924	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
925	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
926	enddo
927	fm(:,nlay+1)=0.
928
929
930	q(:,:)=therm_tke_max(:,:)
931	!!! nrlmd le 16/09/2010
932	do ig=1,ngrid
933	qa(ig,1)=q(ig,1)
934	enddo
935	!!!
936
937	if (1==1) then
938	do k=2,nlay
939	do ig=1,ngrid
940	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
941	& 1.e-5*masse(ig,k)) then
942	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
943	& /(fm(ig,k+1)+detr(ig,k))
944	else
945	qa(ig,k)=q(ig,k)
946	endif
947	if (qa(ig,k).lt.0.) then
948	! print*,'qa<0!!!'
949	endif
950	if (q(ig,k).lt.0.) then
951	! print*,'q<0!!!'
952	endif
953	enddo
954	enddo
955
956	! Calcul du flux subsident
957
958	do k=2,nlay
959	do ig=1,ngrid
960	wqd(ig,k)=fm(ig,k)*q(ig,k)
961	if (wqd(ig,k).lt.0.) then
962	! print*,'wqd<0!!!'
963	endif
964	enddo
965	enddo
966	do ig=1,ngrid
967	wqd(ig,1)=0.
968	wqd(ig,nlay+1)=0.
969	enddo
970
971	! Calcul des tendances
972	do k=1,nlay
973	do ig=1,ngrid
974	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
975	& -wqd(ig,k)+wqd(ig,k+1)) &
976	& *ptimestep/masse(ig,k)
977	enddo
978	enddo
979
980	endif
981
982	therm_tke_max(:,:)=q(:,:)
983
984	return
985	!!! fin nrlmd le 10/04/2012
986	end subroutine thermcell_tke_transport
987
988	END MODULE lmdz_thermcell_main

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

Download in other formats:

Original Format