Context Navigation

thermcell_main.F90 @ 4203

Last change on this file since 4203 was 4143, checked in by dcugnet, 3 years ago
Some variables are renamed or replaced by direct equivalents: iso_indnum -> tracers(:)%iso_iName niso_possibles -> niso iqiso -> iqIsoPha ; index_trac -> itZonIso ok_iso_verif -> isoCheck ntraceurs_zone -> nzone ; ntraciso -> ntiso qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times. Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18). The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
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: 29.3 KB

Line
1
2	! $Id: thermcell_main.F90 4143 2022-05-09 10:35:40Z lguez $
3	!
4	subroutine thermcell_main(itap,ngrid,nlay,ptimestep &
5	& ,pplay,pplev,pphi,debut &
6	& ,pu,pv,pt,po &
7	& ,pduadj,pdvadj,pdtadj,pdoadj &
8	& ,fm0,entr0,detr0,zqta,zqla,lmax &
9	& ,ratqscth,ratqsdiff,zqsatth &
10	& ,zmax0, f0,zw2,fraca,ztv &
11	& ,zpspsk,ztla,zthl,ztva &
12	& ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
13	#ifdef ISO
14	& ,xtpo,xtpdoadj &
15	#endif
16	& )
17
18
19	USE thermcell_ini_mod, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level
20	USE thermcell_ini_mod, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals
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	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
458
459	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
460	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
461	call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
462
463	!------------------------------------------------------------------
464	! On ne prend pas directement les profils issus des calculs precedents
465	! mais on s'autorise genereusement une relaxation vers ceci avec
466	! une constante de temps tau_thermals (typiquement 1800s).
467	!------------------------------------------------------------------
468
469	if (tau_thermals>1.) then
470	lambda=exp(-ptimestep/tau_thermals)
471	fm0=(1.-lambda)fm+lambdafm0
472	entr0=(1.-lambda)entr+lambdaentr0
473	detr0=(1.-lambda)detr+lambdadetr0
474	else
475	fm0=fm
476	entr0=entr
477	detr0=detr
478	endif
479
480	!c------------------------------------------------------------------
481	! calcul du transport vertical
482	!------------------------------------------------------------------
483
484	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
485	& zthl,zdthladj,zta,lev_out)
486	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
487	& po,pdoadj,zoa,lev_out)
488
489	#ifdef ISO
490	! C Risi: on utilise directement la meme routine
491	do ixt=1,ntiso
492	do ll=1,nlay
493	DO ig=1,ngrid
494	xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
495	xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
496	enddo
497	enddo
498	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
499	& xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
500	do ll=1,nlay
501	DO ig=1,ngrid
502	xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
503	enddo
504	enddo
505	enddo
506	#endif
507
508	#ifdef ISO
509	#ifdef ISOVERIF
510	DO ll=1,nlay
511	DO ig=1,ngrid
512	if (iso_eau.gt.0) then
513	call iso_verif_egalite(xtpo(iso_eau,ig,ll), &
514	& po(ig,ll),'thermcell_main 594')
515	call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
516	& pdoadj(ig,ll),'thermcell_main 596')
517	endif
518	if (iso_HDO.gt.0) then
519	call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
520	& /po(ig,ll),'thermcell_main 610')
521	endif
522	enddo
523	enddo !DO ll=1,nlay
524	write(,) 'thermcell_main 600 tmp: apres thermcell_dq'
525	#endif
526	#endif
527
528
529
530	!------------------------------------------------------------------
531	! Calcul de la fraction de l'ascendance
532	!------------------------------------------------------------------
533	do ig=1,ngrid
534	fraca(ig,1)=0.
535	fraca(ig,nlay+1)=0.
536	enddo
537	do l=2,nlay
538	do ig=1,ngrid
539	if (zw2(ig,l).gt.1.e-10) then
540	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
541	else
542	fraca(ig,l)=0.
543	endif
544	enddo
545	enddo
546
547	!------------------------------------------------------------------
548	! calcul du transport vertical du moment horizontal
549	!------------------------------------------------------------------
550
551	!IM 090508
552	if (dvdq == 0 ) then
553
554	! Calcul du transport de V tenant compte d'echange par gradient
555	! de pression horizontal avec l'environnement
556
557	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
558	! & ,fraca*dvdq,zmax &
559	& ,fraca,zmax &
560	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
561
562	else
563
564	! calcul purement conservatif pour le transport de V
565	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
566	& ,zu,pduadj,zua,lev_out)
567	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
568	& ,zv,pdvadj,zva,lev_out)
569
570	endif
571
572	! print*,'13 OK convect8'
573	do l=1,nlay
574	do ig=1,ngrid
575	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
576	enddo
577	enddo
578
579	if (prt_level.ge.1) print*,'14 OK convect8'
580	!------------------------------------------------------------------
581	! Calculs de diagnostiques pour les sorties
582	!------------------------------------------------------------------
583	!calcul de fraca pour les sorties
584
585	if (sorties) then
586	if (prt_level.ge.1) print*,'14a OK convect8'
587	! calcul du niveau de condensation
588	! initialisation
589	do ig=1,ngrid
590	nivcon(ig)=0
591	zcon(ig)=0.
592	enddo
593	!nouveau calcul
594	do ig=1,ngrid
595	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
596	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
597	enddo
598	!IM do k=1,nlay
599	do k=1,nlay-1
600	do ig=1,ngrid
601	if ((pcon(ig).le.pplay(ig,k)) &
602	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
603	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
604	endif
605	enddo
606	enddo
607	!IM
608	ierr=0
609	do ig=1,ngrid
610	if (pcon(ig).le.pplay(ig,nlay)) then
611	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
612	ierr=1
613	endif
614	enddo
615	if (ierr==1) then
616	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
617	CALL abort_physic (modname,abort_message,1)
618	endif
619
620	if (prt_level.ge.1) print*,'14b OK convect8'
621	do k=nlay,1,-1
622	do ig=1,ngrid
623	if (zqla(ig,k).gt.1e-10) then
624	nivcon(ig)=k
625	zcon(ig)=zlev(ig,k)
626	endif
627	enddo
628	enddo
629	if (prt_level.ge.1) print*,'14c OK convect8'
630	!calcul des moments
631	!initialisation
632	do l=1,nlay
633	do ig=1,ngrid
634	q2(ig,l)=0.
635	wth2(ig,l)=0.
636	wth3(ig,l)=0.
637	ratqscth(ig,l)=0.
638	ratqsdiff(ig,l)=0.
639	enddo
640	enddo
641	if (prt_level.ge.1) print*,'14d OK convect8'
642	if (prt_level.ge.10)write(lunout,*) &
643	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
644	do l=1,nlay
645	do ig=1,ngrid
646	zf=fraca(ig,l)
647	zf2=zf/(1.-zf)
648	!
649	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
650	if(zw2(ig,l).gt.1.e-10) then
651	wth2(ig,l)=zf2(zw2(ig,l))*2
652	else
653	wth2(ig,l)=0.
654	endif
655	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
656	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
657	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
658	!test: on calcul q2/po=ratqsc
659	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
660	enddo
661	enddo
662	!calcul des flux: q, thetal et thetav
663	do l=1,nlay
664	do ig=1,ngrid
665	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
666	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
667	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
668	enddo
669	enddo
670
671	!calcul du ratqscdiff
672	if (prt_level.ge.1) print*,'14e OK convect8'
673	var=0.
674	vardiff=0.
675	ratqsdiff(:,:)=0.
676
677	do l=1,nlay
678	do ig=1,ngrid
679	if (l<=lalim(ig)) then
680	var=var+alim_star(ig,l)zqta(ig,l)1000.
681	endif
682	enddo
683	enddo
684
685	if (prt_level.ge.1) print*,'14f OK convect8'
686
687	do l=1,nlay
688	do ig=1,ngrid
689	if (l<=lalim(ig)) then
690	zf=fraca(ig,l)
691	zf2=zf/(1.-zf)
692	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
693	endif
694	enddo
695	enddo
696
697	if (prt_level.ge.1) print*,'14g OK convect8'
698	do l=1,nlay
699	do ig=1,ngrid
700	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
701	enddo
702	enddo
703	endif
704
705	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
706
707	RETURN
708	end subroutine thermcell_main
709
710	!=============================================================================
711	!/////////////////////////////////////////////////////////////////////////////
712	!=============================================================================
713	subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,po,ztva, & ! in
714	& zqla,f_star,zw2,comment) ! in
715	!=============================================================================
716	USE thermcell_ini_mod, ONLY: prt_level
717	IMPLICIT NONE
718
719	integer i, k, ngrid,nlay
720	real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,po,ztva,zqla
721	real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
722	integer, intent(in), dimension(ngrid) :: long
723	real seuil
724	character*21 comment
725
726	seuil=0.25
727
728	if (prt_level.ge.1) THEN
729	print*,'WARNING !!! TEST ',comment
730	endif
731	return
732
733	! test sur la hauteur des thermiques ...
734	do i=1,ngrid
735	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
736	if (prt_level.ge.10) then
737	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
738	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
739	do k=1,nlay
740	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)
741	enddo
742	endif
743	enddo
744
745
746	return
747	end
748
749	! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
750	! On transporte pbl_tke pour donner therm_tke
751	! Copie conforme de la subroutine DTKE dans physiq.F ecrite par Frederic Hourdin
752
753	!=======================================================================
754	!///////////////////////////////////////////////////////////////////////
755	!=======================================================================
756
757	subroutine thermcell_tke_transport( &
758	& ngrid,nlay,ptimestep,fm0,entr0,rg,pplev, & ! in
759	& therm_tke_max) ! out
760	USE thermcell_ini_mod, ONLY: prt_level
761	implicit none
762
763	!=======================================================================
764	!
765	! Calcul du transport verticale dans la couche limite en presence
766	! de "thermiques" explicitement representes
767	! calcul du dq/dt une fois qu'on connait les ascendances
768	!
769	!=======================================================================
770
771	integer ngrid,nlay
772
773	real, intent(in) :: ptimestep
774	real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
775	real, intent(in), dimension(ngrid,nlay) :: entr0
776	real, intent(in) :: rg
777	real, intent(out), dimension(ngrid,nlay) :: therm_tke_max
778
779	real detr0(ngrid,nlay)
780	real masse0(ngrid,nlay)
781	real masse(ngrid,nlay),fm(ngrid,nlay+1)
782	real entr(ngrid,nlay)
783	real q(ngrid,nlay)
784	integer lev_out ! niveau pour les print
785
786	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
787	integer ig,k
788
789
790	lev_out=0
791
792
793	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
794
795	! calcul du detrainement
796	do k=1,nlay
797	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
798	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
799	enddo
800
801
802	! Decalage vertical des entrainements et detrainements.
803	masse(:,1)=0.5*masse0(:,1)
804	entr(:,1)=0.5*entr0(:,1)
805	detr(:,1)=0.5*detr0(:,1)
806	fm(:,1)=0.
807	do k=1,nlay-1
808	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
809	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
810	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
811	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
812	enddo
813	fm(:,nlay+1)=0.
814
815
816	q(:,:)=therm_tke_max(:,:)
817	!!! nrlmd le 16/09/2010
818	do ig=1,ngrid
819	qa(ig,1)=q(ig,1)
820	enddo
821	!!!
822
823	if (1==1) then
824	do k=2,nlay
825	do ig=1,ngrid
826	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
827	& 1.e-5*masse(ig,k)) then
828	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
829	& /(fm(ig,k+1)+detr(ig,k))
830	else
831	qa(ig,k)=q(ig,k)
832	endif
833	if (qa(ig,k).lt.0.) then
834	! print*,'qa<0!!!'
835	endif
836	if (q(ig,k).lt.0.) then
837	! print*,'q<0!!!'
838	endif
839	enddo
840	enddo
841
842	! Calcul du flux subsident
843
844	do k=2,nlay
845	do ig=1,ngrid
846	wqd(ig,k)=fm(ig,k)*q(ig,k)
847	if (wqd(ig,k).lt.0.) then
848	! print*,'wqd<0!!!'
849	endif
850	enddo
851	enddo
852	do ig=1,ngrid
853	wqd(ig,1)=0.
854	wqd(ig,nlay+1)=0.
855	enddo
856
857	! Calcul des tendances
858	do k=1,nlay
859	do ig=1,ngrid
860	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
861	& -wqd(ig,k)+wqd(ig,k+1)) &
862	& *ptimestep/masse(ig,k)
863	enddo
864	enddo
865
866	endif
867
868	therm_tke_max(:,:)=q(:,:)
869
870	return
871	!!! fin nrlmd le 10/04/2012
872	end
873

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

Download in other formats:

Original Format