Context Navigation

lmdz_thermcell_main.F90 @ 5218

Last change on this file since 5218 was 4843, checked in by crio, 9 months ago
Nouvelle formulation du strig et correction thermiques montent trop haut
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.1 KB

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

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

Download in other formats:

Original Format