Context Navigation

thermcell_main.F90 @ 4413

Last change on this file since 4413 was 4413, checked in by evignon, 16 months ago
pour les mini-projets: controle des coefficients pour le detrainement et l'entrainement des downdrafts dans les .def
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.3 KB

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

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

Download in other formats:

Original Format