Context Navigation

lmdz_thermcell_main.F90 @ 5006

Last change on this file since 5006 was 4843, checked in by crio, 10 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