Context Navigation

thermcell_main.F90 @ 4396

Last change on this file since 4396 was 4396, checked in by evignon, 20 months ago
Attention les yeux, commission qui enterine la premiere version operationnelle des downdrafts dans LMDZ. Speciale dedicace a Fleur
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: 30.1 KB

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

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

Download in other formats:

Original Format