Context Navigation

thermcell_main.F90 @ 93

Last change on this file since 93 was 93, checked in by lfita, 10 years ago
Including checks for wth3, NaN source?!
File size: 48.5 KB

Line
1	!
2	! $Id: thermcell_main.F90 1795 2013-07-18 08:20:28Z emillour $
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	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
11	& ,zmax0, f0,zw2,fraca,ztv &
12	& ,zpspsk,ztla,zthl &
13	!!! nrlmd le 10/04/2012
14	& ,pbl_tke,pctsrf,omega,airephy &
15	& ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
16	& ,n2,s2,ale_bl_stat &
17	& ,therm_tke_max,env_tke_max &
18	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
19	& ,alp_bl_conv,alp_bl_stat &
20	!!! fin nrlmd le 10/04/2012
21	& ,ztva )
22
23	USE dimphy
24	USE ioipsl
25	USE comgeomphy , ONLY:rlond,rlatd
26	USE indice_sol_mod
27	IMPLICIT NONE
28
29	!=======================================================================
30	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
31	! Version du 09.02.07
32	! Calcul du transport vertical dans la couche limite en presence
33	! de "thermiques" explicitement representes avec processus nuageux
34	!
35	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
36	!
37	! le thermique est suppose homogene et dissipe par melange avec
38	! son environnement. la longueur l_mix controle l'efficacite du
39	! melange
40	!
41	! Le calcul du transport des differentes especes se fait en prenant
42	! en compte:
43	! 1. un flux de masse montant
44	! 2. un flux de masse descendant
45	! 3. un entrainement
46	! 4. un detrainement
47	!
48	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
49	! Introduction of an implicit computation of vertical advection in
50	! the environment of thermal plumes in thermcell_dq
51	! impl = 0 : explicit, 1 : implicit, -1 : old version
52	! controled by iflag_thermals =
53	! 15, 16 run with impl=-1 : numerical convergence with NPv3
54	! 17, 18 run with impl=1 : more stable
55	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
56	!
57	!=======================================================================
58
59
60	!-----------------------------------------------------------------------
61	! declarations:
62	! -------------
63
64	#include "dimensions.h"
65	#include "YOMCST.h"
66	#include "YOETHF.h"
67	#include "FCTTRE.h"
68	#include "iniprint.h"
69	#include "thermcell.h"
70
71	! arguments:
72	! ----------
73
74	!IM 140508
75	INTEGER itap
76
77	INTEGER ngrid,nlay
78	real ptimestep
79	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
80	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
81	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
82	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
83	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
84	real pphi(ngrid,nlay)
85
86	! local:
87	! ------
88
89	integer icount
90
91	integer, save :: dvdq=1,dqimpl=-1
92	!$OMP THREADPRIVATE(dvdq,dqimpl)
93	data icount/0/
94	save icount
95	!$OMP THREADPRIVATE(icount)
96
97	integer,save :: igout=1
98	!$OMP THREADPRIVATE(igout)
99	integer,save :: lunout1=6
100	!$OMP THREADPRIVATE(lunout1)
101	integer,save :: lev_out=10
102	!$OMP THREADPRIVATE(lev_out)
103
104	REAL susqr2pi, Reuler
105
106	INTEGER ig,k,l,ll,ierr
107	real zsortie1d(ngrid)
108	INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid)
109	INTEGER lmix(ngrid)
110	INTEGER lmix_bis(ngrid)
111	real linter(ngrid)
112	real zmix(ngrid)
113	real zmax(ngrid),zw2(ngrid,nlay+1),ztva(ngrid,nlay),zw_est(ngrid,nlay+1),ztva_est(ngrid,nlay)
114	! real fraca(ngrid,nlay)
115
116	real zmax_sec(ngrid)
117	!on garde le zmax du pas de temps precedent
118	real zmax0(ngrid)
119	!FH/IM save zmax0
120
121	real lambda
122
123	real zlev(ngrid,nlay+1),zlay(ngrid,nlay)
124	real deltaz(ngrid,nlay)
125	REAL zh(ngrid,nlay)
126	real zthl(ngrid,nlay),zdthladj(ngrid,nlay)
127	REAL ztv(ngrid,nlay)
128	real zu(ngrid,nlay),zv(ngrid,nlay),zo(ngrid,nlay)
129	real zl(ngrid,nlay)
130	real zsortie(ngrid,nlay)
131	real zva(ngrid,nlay)
132	real zua(ngrid,nlay)
133	real zoa(ngrid,nlay)
134
135	real zta(ngrid,nlay)
136	real zha(ngrid,nlay)
137	real fraca(ngrid,nlay+1)
138	real zf,zf2
139	real thetath2(ngrid,nlay),wth2(ngrid,nlay),wth3(ngrid,nlay)
140	real q2(ngrid,nlay)
141	! FH probleme de dimensionnement avec l'allocation dynamique
142	! common/comtherm/thetath2,wth2
143	real wq(ngrid,nlay)
144	real wthl(ngrid,nlay)
145	real wthv(ngrid,nlay)
146
147	real ratqscth(ngrid,nlay)
148	real var
149	real vardiff
150	real ratqsdiff(ngrid,nlay)
151
152	logical sorties
153	real rho(ngrid,nlay),rhobarz(ngrid,nlay),masse(ngrid,nlay)
154	real zpspsk(ngrid,nlay)
155
156	real wmax(ngrid)
157	real wmax_tmp(ngrid)
158	real wmax_sec(ngrid)
159	real fm0(ngrid,nlay+1),entr0(ngrid,nlay),detr0(ngrid,nlay)
160	real fm(ngrid,nlay+1),entr(ngrid,nlay),detr(ngrid,nlay)
161
162	real ztla(ngrid,nlay),zqla(ngrid,nlay),zqta(ngrid,nlay)
163	!niveau de condensation
164	integer nivcon(ngrid)
165	real zcon(ngrid)
166	REAL CHI
167	real zcon2(ngrid)
168	real pcon(ngrid)
169	real zqsat(ngrid,nlay)
170	real zqsatth(ngrid,nlay)
171
172	real f_star(ngrid,nlay+1),entr_star(ngrid,nlay)
173	real detr_star(ngrid,nlay)
174	real alim_star_tot(ngrid)
175	real alim_star(ngrid,nlay)
176	real alim_star_clos(ngrid,nlay)
177	real f(ngrid), f0(ngrid)
178	!FH/IM save f0
179	real zlevinter(ngrid)
180	logical debut
181	real seuil
182	real csc(ngrid,nlay)
183
184	!!! nrlmd le 10/04/2012
185
186	!------Entrées
187	real pbl_tke(ngrid,nlay+1,nbsrf)
188	real pctsrf(ngrid,nbsrf)
189	real omega(ngrid,nlay)
190	real airephy(ngrid)
191	!------Sorties
192	real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
193	real therm_tke_max0(ngrid),env_tke_max0(ngrid)
194	real n2(ngrid),s2(ngrid)
195	real ale_bl_stat(ngrid)
196	real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay)
197	real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid)
198	!------Local
199	integer nsrf
200	real rhobarz0(ngrid) ! Densité au LCL
201	logical ok_lcl(ngrid) ! Existence du LCL des thermiques
202	integer klcl(ngrid) ! Niveau du LCL
203	real interp(ngrid) ! Coef d'interpolation pour le LCL
204	!--Triggering
205	real Su ! Surface unité: celle d'un updraft élémentaire
206	parameter(Su=4e4)
207	real hcoef ! Coefficient directeur pour le calcul de s2
208	parameter(hcoef=1)
209	real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
210	parameter(hmincoef=0.3)
211	real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1s1/(fraca0Sd)
212	parameter(eps1=0.3)
213	real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2
214	real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
215	real zmax_moy_coef
216	parameter(zmax_moy_coef=0.33)
217	real depth(ngrid) ! Epaisseur moyenne du cumulus
218	real w_max(ngrid) ! Vitesse max statistique
219	real s_max(ngrid)
220	!--Closure
221	real pbl_tke_max(ngrid,nlay) ! Profil de TKE moyenne
222	real pbl_tke_max0(ngrid) ! TKE moyenne au LCL
223	real w_ls(ngrid,nlay) ! Vitesse verticale grande échelle (m/s)
224	real coef_m ! On considère un rendement pour alp_bl_fluct_m
225	parameter(coef_m=1.)
226	real coef_tke ! On considère un rendement pour alp_bl_fluct_tke
227	parameter(coef_tke=1.)
228
229	!!! fin nrlmd le 10/04/2012
230
231	!
232	!nouvelles variables pour la convection
233	real Ale_bl(ngrid)
234	real Alp_bl(ngrid)
235	real alp_int(ngrid),dp_int(ngrid),zdp
236	real ale_int(ngrid)
237	integer n_int(ngrid)
238	real fm_tot(ngrid)
239	real wght_th(ngrid,nlay)
240	integer lalim_conv(ngrid)
241	!v1d logical therm
242	!v1d save therm
243
244	character*2 str2
245	character*10 str10
246
247	character (len=20) :: modname='thermcell_main'
248	character (len=80) :: abort_message
249
250	EXTERNAL SCOPY
251	!
252
253	! Lluis
254	INTEGER :: llp
255	CHARACTER(LEN=50) :: lvarname, lfname
256	REAL :: largest
257
258	llp = 734
259	lfname = 'physiq'
260	largest = 10.e5
261
262	! L. Fita, LMD July 2014. Initializing variables
263	zdthladj = 0.
264	pbl_tke_max0 = 0.
265
266	!-----------------------------------------------------------------------
267	! initialisation:
268	! ---------------
269	!
270
271	seuil=0.25
272
273	if (debut) then
274	! call getin('dvdq',dvdq)
275	! call getin('dqimpl',dqimpl)
276
277	if (iflag_thermals==15.or.iflag_thermals==16) then
278	dvdq=0
279	dqimpl=-1
280	else
281	dvdq=1
282	dqimpl=1
283	endif
284
285	fm0=0.
286	entr0=0.
287	detr0=0.
288
289
290	#undef wrgrads_thermcell
291	#ifdef wrgrads_thermcell
292	! Initialisation des sorties grads pour les thermiques.
293	! Pour l'instant en 1D sur le point igout.
294	! Utilise par thermcell_out3d.h
295	str10='therm'
296	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
297	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
298	& ptimestep,str10,'therm ')
299	#endif
300
301
302
303	endif
304
305	fm=0. ; entr=0. ; detr=0.
306
307
308	icount=icount+1
309
310	!IM 090508 beg
311	!print*,'====================================================================='
312	!print*,'====================================================================='
313	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
314	!print*,'====================================================================='
315	!print*,'====================================================================='
316	!IM 090508 end
317
318	if (prt_level.ge.1) print*,'thermcell_main V4'
319
320	sorties=.true.
321	IF(ngrid.NE.ngrid) THEN
322	PRINT*
323	PRINT*,'STOP dans convadj'
324	PRINT*,'ngrid =',ngrid
325	PRINT*,'ngrid =',ngrid
326	ENDIF
327	!
328	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
329	do ig=1,ngrid
330	f0(ig)=max(f0(ig),1.e-2)
331	zmax0(ig)=max(zmax0(ig),40.)
332	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
333	enddo
334
335	if (prt_level.ge.20) then
336	do ig=1,ngrid
337	print*,'th_main ig f0',ig,f0(ig)
338	enddo
339	endif
340	!-----------------------------------------------------------------------
341	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
342	! --------------------------------------------------------------------
343	!
344	lfname='before thermcell_env'
345	lvarname = 'pt'
346	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
347	lvarname = 'pdtadj'
348	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
349
350	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
351	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
352
353	lfname='after thermcell_env'
354	lvarname = 'pt'
355	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
356	lvarname = 'pdtadj'
357	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
358
359	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
360
361	!------------------------------------------------------------------------
362	! --------------------
363	!
364	!
365	! + + + + + + + + + + +
366	!
367	!
368	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
369	! wh,wt,wo ...
370	!
371	! + + + + + + + + + + + zh,zu,zv,zo,rho
372	!
373	!
374	! -------------------- zlev(1)
375	! \\\\\\\\\\\\\\\\\\\\
376	!
377	!
378
379	!-----------------------------------------------------------------------
380	! Calcul des altitudes des couches
381	!-----------------------------------------------------------------------
382
383	do l=2,nlay
384	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
385	enddo
386	zlev(:,1)=0.
387	zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
388	do l=1,nlay
389	zlay(:,l)=pphi(:,l)/RG
390	enddo
391	!calcul de l epaisseur des couches
392	do l=1,nlay
393	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
394	enddo
395
396	! print*,'2 OK convect8'
397	!-----------------------------------------------------------------------
398	! Calcul des densites
399	!-----------------------------------------------------------------------
400
401	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
402
403	if (prt_level.ge.10)write(lunout,*) &
404	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
405	rhobarz(:,1)=rho(:,1)
406
407	do l=2,nlay
408	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
409	enddo
410
411	!calcul de la masse
412	do l=1,nlay
413	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
414	enddo
415
416	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
417
418	!------------------------------------------------------------------
419	!
420	! /\|\
421	! -------- \| F_k+1 -------
422	! ----> D_k
423	! /\|\ <---- E_k , A_k
424	! -------- \| F_k ---------
425	! ----> D_k-1
426	! <---- E_k-1 , A_k-1
427	!
428	!
429	!
430	!
431	!
432	! ---------------------------
433	!
434	! ----- F_lmax+1=0 ---------- \
435	! lmax (zmax) \|
436	! --------------------------- \|
437	! \|
438	! --------------------------- \|
439	! \|
440	! --------------------------- \|
441	! \|
442	! --------------------------- \|
443	! \|
444	! --------------------------- \|
445	! \| E
446	! --------------------------- \| D
447	! \|
448	! --------------------------- \|
449	! \|
450	! --------------------------- \ \|
451	! lalim \| \|
452	! --------------------------- \| \|
453	! \| \|
454	! --------------------------- \| \|
455	! \| A \|
456	! --------------------------- \| \|
457	! \| \|
458	! --------------------------- \| \|
459	! lmin (=1 pour le moment) \| \|
460	! ----- F_lmin=0 ------------ / /
461	!
462	! ---------------------------
463	! //////////////////////////
464	!
465	!
466	!=============================================================================
467	! Calculs initiaux ne faisant pas intervenir les changements de phase
468	!=============================================================================
469
470	!------------------------------------------------------------------
471	! 1. alim_star est le profil vertical de l'alimentation a la base du
472	! panache thermique, calcule a partir de la flotabilite de l'air sec
473	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
474	!------------------------------------------------------------------
475	!
476	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
477	lmin=1
478
479	!-----------------------------------------------------------------------------
480	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
481	! panache sec conservatif (e=d=0) alimente selon alim_star
482	! Il s'agit d'un calcul de type CAPE
483	! zmax_sec est utilise pour determiner la geometrie du thermique.
484	!------------------------------------------------------------------------------
485	!---------------------------------------------------------------------------------
486	!calcul du melange et des variables dans le thermique
487	!--------------------------------------------------------------------------------
488	!
489	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
490	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
491
492	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
493	! de la variable iflag_thermals
494
495	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
496	if (iflag_thermals_ed<=9) then
497	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
498	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
499	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
500	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
501	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
502	& ,lev_out,lunout1,igout)
503
504	elseif (iflag_thermals_ed>9) then
505	! print*,'THERM RIO et al 2010, version d Arnaud'
506	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
507	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
508	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
509	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
510	& ,lev_out,lunout1,igout)
511
512	endif
513	lfname='after thermcell_plume'
514	lvarname = 'pt'
515	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
516	lvarname = 'pdtadj'
517	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
518
519	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
520
521	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
522	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
523
524	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
525	if (prt_level.ge.10) then
526	write(lunout1,*) 'Dans thermcell_main 2'
527	write(lunout1,*) 'lmin ',lmin(igout)
528	write(lunout1,*) 'lalim ',lalim(igout)
529	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
530	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
531	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
532	endif
533
534	!-------------------------------------------------------------------------------
535	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
536	!-------------------------------------------------------------------------------
537	!
538	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
539	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
540	! Attention, w2 est transforme en sa racine carree dans cette routine
541	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
542	wmax_tmp=0.
543	do l=1,nlay
544	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
545	enddo
546	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
547
548
549
550	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
551	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
552	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
553	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
554
555	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
556
557	!-------------------------------------------------------------------------------
558	! Fermeture,determination de f
559	!-------------------------------------------------------------------------------
560	!
561	!
562	!! write(lunout,*)'THERM NOUVEAU XXXXX'
563	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
564	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
565
566	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
567	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
568
569	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
570	if (prt_level.ge.10) then
571	write(lunout1,*) 'Dans thermcell_main 1b'
572	write(lunout1,*) 'lmin ',lmin(igout)
573	write(lunout1,*) 'lalim ',lalim(igout)
574	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
575	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
576	& ,l=1,lalim(igout)+4)
577	endif
578
579
580
581
582	! Choix de la fonction d'alimentation utilisee pour la fermeture.
583	! Apparemment sans importance
584	alim_star_clos(:,:)=alim_star(:,:)
585	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
586
587	! Appel avec la version seche
588	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
589	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
590
591	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
592	! Appel avec les zmax et wmax tenant compte de la condensation
593	! Semble moins bien marcher
594	! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
595	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
596	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
597
598	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
599
600	if (tau_thermals>1.) then
601	lambda=exp(-ptimestep/tau_thermals)
602	f0=(1.-lambda)f+lambdaf0
603	else
604	f0=f
605	endif
606
607	! Test valable seulement en 1D mais pas genant
608	if (.not. (f0(1).ge.0.) ) then
609	abort_message = '.not. (f0(1).ge.0.)'
610	CALL abort_gcm (modname,abort_message,1)
611	endif
612
613	!-------------------------------------------------------------------------------
614	!deduction des flux
615	!-------------------------------------------------------------------------------
616	lfname='before flux'
617	lvarname = 'fm'
618	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
619	lvarname = 'fm'
620	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
621	lvarname = 'entr'
622	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
623	lvarname = 'detr'
624	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
625	lvarname = 'zgla'
626	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
627
628	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
629	& lalim,lmax,alim_star, &
630	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
631	& detr,zqla,lev_out,lunout1,igout)
632	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
633	lfname='after flux'
634	lvarname = 'fm'
635	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
636
637	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
638	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
639	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
640
641	!------------------------------------------------------------------
642	! On ne prend pas directement les profils issus des calculs precedents
643	! mais on s'autorise genereusement une relaxation vers ceci avec
644	! une constante de temps tau_thermals (typiquement 1800s).
645	!------------------------------------------------------------------
646
647	if (tau_thermals>1.) then
648	lambda=exp(-ptimestep/tau_thermals)
649	fm0=(1.-lambda)fm+lambdafm0
650	entr0=(1.-lambda)entr+lambdaentr0
651	detr0=(1.-lambda)detr+lambdadetr0
652	else
653	fm0=fm
654	entr0=entr
655	detr0=detr
656	endif
657
658	!c------------------------------------------------------------------
659	! calcul du transport vertical
660	!------------------------------------------------------------------
661
662	lfname='before transport_vertical'
663	lvarname = 'zdthladj'
664	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
665
666	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
667	& zthl,zdthladj,zta,lev_out)
668	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
669	& po,pdoadj,zoa,lev_out)
670	lfname='after transport_vertical'
671	lvarname = 'zdthladj'
672	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
673
674	!------------------------------------------------------------------
675	! Calcul de la fraction de l'ascendance
676	!------------------------------------------------------------------
677	do ig=1,ngrid
678	fraca(ig,1)=0.
679	fraca(ig,nlay+1)=0.
680	enddo
681	do l=2,nlay
682	do ig=1,ngrid
683	if (zw2(ig,l).gt.1.e-10) then
684	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
685	else
686	fraca(ig,l)=0.
687	endif
688	enddo
689	enddo
690
691	!------------------------------------------------------------------
692	! calcul du transport vertical du moment horizontal
693	!------------------------------------------------------------------
694	lfname='before thermcell_dv2'
695	lvarname = 'pt'
696	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
697	lvarname = 'pdtadj'
698	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
699
700	!IM 090508
701	if (dvdq == 0 ) then
702
703	! Calcul du transport de V tenant compte d'echange par gradient
704	! de pression horizontal avec l'environnement
705
706	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
707	! & ,fraca*dvdq,zmax &
708	& ,fraca,zmax &
709	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
710
711	else
712
713	! calcul purement conservatif pour le transport de V
714	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
715	& ,zu,pduadj,zua,lev_out)
716	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
717	& ,zv,pdvadj,zva,lev_out)
718
719	endif
720
721	! print*,'13 OK convect8'
722	do l=1,nlay
723	do ig=1,ngrid
724	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
725	enddo
726	enddo
727	lfname='after thermcell_dv2'
728	lvarname = 'pt'
729	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
730	lvarname = 'pdtadj'
731	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
732	lvarname = 'zdthladj'
733	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
734	lvarname = 'zpspsk'
735	CALL check_var3D(lfname, lvarname, zpspsk, ngrid, nlay, largest, .FALSE.)
736
737	if (prt_level.ge.1) print*,'14 OK convect8'
738	!------------------------------------------------------------------
739	! Calculs de diagnostiques pour les sorties
740	!------------------------------------------------------------------
741	!calcul de fraca pour les sorties
742
743	if (sorties) then
744	if (prt_level.ge.1) print*,'14a OK convect8'
745	! calcul du niveau de condensation
746	! initialisation
747	do ig=1,ngrid
748	nivcon(ig)=0
749	zcon(ig)=0.
750	enddo
751	!nouveau calcul
752	do ig=1,ngrid
753	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
754	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
755	enddo
756	!IM do k=1,nlay
757	do k=1,nlay-1
758	do ig=1,ngrid
759	if ((pcon(ig).le.pplay(ig,k)) &
760	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
761	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
762	endif
763	enddo
764	enddo
765	!IM
766	ierr=0
767	do ig=1,ngrid
768	if (pcon(ig).le.pplay(ig,nlay)) then
769	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
770	ierr=1
771	endif
772	enddo
773	if (ierr==1) then
774	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
775	CALL abort_gcm (modname,abort_message,1)
776	endif
777
778	if (prt_level.ge.1) print*,'14b OK convect8'
779	do k=nlay,1,-1
780	do ig=1,ngrid
781	if (zqla(ig,k).gt.1e-10) then
782	nivcon(ig)=k
783	zcon(ig)=zlev(ig,k)
784	endif
785	enddo
786	enddo
787	if (prt_level.ge.1) print*,'14c OK convect8'
788	!calcul des moments
789	!initialisation
790	do l=1,nlay
791	do ig=1,ngrid
792	q2(ig,l)=0.
793	wth2(ig,l)=0.
794	wth3(ig,l)=0.
795	ratqscth(ig,l)=0.
796	ratqsdiff(ig,l)=0.
797	enddo
798	enddo
799	if (prt_level.ge.1) print*,'14d OK convect8'
800	if (prt_level.ge.10)write(lunout,*) &
801	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
802	do l=1,nlay
803	do ig=1,ngrid
804	zf=fraca(ig,l)
805	zf2=zf/(1.-zf)
806	!
807	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
808	if(zw2(ig,l).gt.1.e-10) then
809	wth2(ig,l)=zf2(zw2(ig,l))*2
810	else
811	wth2(ig,l)=0.
812	endif
813	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
814	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
815	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
816	!test: on calcul q2/po=ratqsc
817	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
818	enddo
819	enddo
820	lfname='thermcell_main calculation of wth3'
821	lvarname = 'wth3'
822	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
823	lvarname = 'fraca'
824	CALL check_var3D(lfname, lvarname, fraca, ngrid, nlay, largest, .FALSE.)
825	lvarname = 'zw2'
826	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay, largest, .FALSE.)
827	lvarname = '1-fraca'
828	CALL check_var3D(lfname, lvarname, 1./(1.-fraca), ngrid, nlay, largest, .FALSE.)
829	!calcul des flux: q, thetal et thetav
830	do l=1,nlay
831	do ig=1,ngrid
832	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
833	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
834	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
835	enddo
836	enddo
837	!
838
839	!!! nrlmd le 10/04/2012
840
841	!------------Test sur le LCL des thermiques
842	do ig=1,ngrid
843	ok_lcl(ig)=.false.
844	if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
845	enddo
846
847	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
848	do l=1,nlay-1
849	do ig=1,ngrid
850	if (ok_lcl(ig)) then
851	if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
852	klcl(ig)=l
853	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
854	endif
855	endif
856	enddo
857	enddo
858
859	!------------Hauteur des thermiques
860	!!jyg le 27/04/2012
861	!! do ig =1,ngrid
862	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
863	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
864	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
865	!! zmax(ig)=pphi(ig,lmax(ig))/rg
866	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
867	!! enddo
868	do ig =1,ngrid
869	zmax(ig)=pphi(ig,lmax(ig))/rg
870	if (ok_lcl(ig)) then
871	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
872	& -rhobarz(ig,klcl(ig)))*interp(ig)
873	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
874	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
875	else
876	rhobarz0(ig)=0.
877	zlcl(ig)=zmax(ig)
878	endif
879	enddo
880	!!jyg fin
881
882	!------------Calcul des propriétés du thermique au LCL
883	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
884
885	!-----Initialisation de la TKE moyenne
886	do l=1,nlay
887	do ig=1,ngrid
888	pbl_tke_max(ig,l)=0.
889	enddo
890	enddo
891
892	!-----Calcul de la TKE moyenne
893	do nsrf=1,nbsrf
894	do l=1,nlay
895	do ig=1,ngrid
896	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
897	enddo
898	enddo
899	enddo
900
901	!-----Initialisations des TKE dans et hors des thermiques
902	do l=1,nlay
903	do ig=1,ngrid
904	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
905	env_tke_max(ig,l)=pbl_tke_max(ig,l)
906	enddo
907	enddo
908
909	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
910	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
911	& rg,pplev,therm_tke_max)
912	! print *,' thermcell_tke_transport -> ' !!jyg
913
914	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
915	do l=1,nlay
916	do ig=1,ngrid
917	pbl_tke_max(ig,l)=fraca(ig,l)therm_tke_max(ig,l)+(1.-fraca(ig,l))env_tke_max(ig,l) ! Recalcul de TKE moyenne aprés transport de TKE_TH
918	env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l)) ! Recalcul de TKE dans l'environnement aprés transport de TKE_TH
919	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
920	enddo
921	enddo
922	! print *,' apres w_ls = ' !!jyg
923
924	do ig=1,ngrid
925	if (ok_lcl(ig)) then
926	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
927	& -fraca(ig,klcl(ig)))*interp(ig)
928	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
929	& -zw2(ig,klcl(ig)))*interp(ig)
930	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
931	& -w_ls(ig,klcl(ig)))*interp(ig)
932	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
933	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
934	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
935	& -env_tke_max(ig,klcl(ig)))*interp(ig)
936	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
937	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
938	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
939	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
940	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
941	else
942	fraca0(ig)=0.
943	w0(ig)=0.
944	!!jyg le 27/04/2012
945	!! zlcl(ig)=0.
946	!!
947	endif
948	enddo
949
950	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
951	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
952
953	!------------Triggering------------------
954	IF (iflag_trig_bl.ge.1) THEN
955
956	!-----Initialisations
957	depth(:)=0.
958	n2(:)=0.
959	s2(:)=0.
960	s_max(:)=0.
961
962	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
963	do ig=1,ngrid
964	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
965	depth(ig)=zmax_moy(ig)-zlcl(ig)
966	hmin(ig)=hmincoef*zlcl(ig)
967	if (depth(ig).ge.10.) then
968	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
969	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
970	!!
971	!!jyg le 27/04/2012
972	!! s_max(ig)=s2(ig)*log(n2(ig))
973	!! if (n2(ig) .lt. 1) s_max(ig)=0.
974	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
975	!!fin jyg
976	else
977	s2(ig)=0.
978	n2(ig)=0.
979	s_max(ig)=0.
980	endif
981	enddo
982	! print *,'avant Calcul de Wmax ' !!jyg
983
984	!-----Calcul de Wmax et ALE_BL_STAT associée
985	!!jyg le 30/04/2012
986	!! do ig=1,ngrid
987	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
988	!! w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/su)-log(2.3.14)-log(2.log(s_max(ig)/su)-log(2.*3.14))))
989	!! ale_bl_stat(ig)=0.5w_max(ig)*2
990	!! else
991	!! w_max(ig)=0.
992	!! ale_bl_stat(ig)=0.
993	!! endif
994	!! enddo
995	susqr2pi=susqrt(2.Rpi)
996	Reuler=exp(1.)
997	do ig=1,ngrid
998	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
999	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
1000	ale_bl_stat(ig)=0.5w_max(ig)*2
1001	else
1002	w_max(ig)=0.
1003	ale_bl_stat(ig)=0.
1004	endif
1005	enddo
1006
1007	ENDIF ! iflag_trig_bl
1008	! print *,'ENDIF iflag_trig_bl' !!jyg
1009	lfname='thermcell_main before closure'
1010	lvarname = 'pt'
1011	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1012	lvarname = 'pdtadj'
1013	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1014	lvarname = 'alp_bl_det'
1015	CALL check_var(lfname, lvarname, alp_bl_det,ngrid, largest, .FALSE.)
1016	lvarname = 'rhobarz0'
1017	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1018	lvarname = 'fraca0'
1019	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1020	lvarname = 'w_conv'
1021	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1022	lvarname = 'w0'
1023	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1024	lvarname = 'alp_bl_fluct_m'
1025	CALL check_var(lfname, lvarname, alp_bl_fluct_m,ngrid, largest, .FALSE.)
1026	lvarname = 'alp_bl_fluct_tke'
1027	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,ngrid, largest, .FALSE.)
1028	lvarname = 'therm_tke_max0'
1029	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1030	lvarname = 'env_tke_max0'
1031	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1032	lvarname = 'pbl_tke_max0'
1033	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1034	lvarname = 'alp_bl_conv'
1035	CALL check_var(lfname, lvarname, alp_bl_conv,ngrid, largest, .FALSE.)
1036	lvarname = 'alp_bl_stat'
1037	CALL check_var(lfname, lvarname, alp_bl_stat,ngrid, largest, .FALSE.)
1038	!------------Closure------------------
1039
1040	IF (iflag_clos_bl.ge.1) THEN
1041
1042	!-----Calcul de ALP_BL_STAT
1043	do ig=1,ngrid
1044	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
1045	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
1046	& (w0(ig)**2)
1047	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
1048	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
1049	if (iflag_clos_bl.ge.2) then
1050	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
1051	& (w0(ig)**2)
1052	else
1053	alp_bl_conv(ig)=0.
1054	endif
1055	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
1056	enddo
1057
1058	!-----Sécurité ALP infinie
1059	do ig=1,ngrid
1060	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
1061	enddo
1062
1063	ENDIF ! (iflag_clos_bl.ge.1)
1064	lfname='after closure'
1065	lvarname = 'pt'
1066	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1067	lvarname = 'pdtadj'
1068	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1069	lvarname = 'alp_bl_det'
1070	CALL check_var(lfname, lvarname, alp_bl_det,ngrid, largest, .FALSE.)
1071	lvarname = 'rhobarz0'
1072	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1073	lvarname = 'fraca0'
1074	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1075	lvarname = 'w_conv'
1076	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1077	lvarname = 'w0'
1078	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1079	lvarname = 'alp_bl_fluct_m'
1080	CALL check_var(lfname, lvarname, alp_bl_fluct_m,ngrid, largest, .FALSE.)
1081	lvarname = 'alp_bl_fluct_tke'
1082	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,ngrid, largest, .FALSE.)
1083	lvarname = 'therm_tke_max0'
1084	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1085	lvarname = 'env_tke_max0'
1086	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1087	lvarname = 'pbl_tke_max0'
1088	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1089	lvarname = 'alp_bl_conv'
1090	CALL check_var(lfname, lvarname, alp_bl_conv,ngrid, largest, .FALSE.)
1091	lvarname = 'alp_bl_stat'
1092	CALL check_var(lfname, lvarname, alp_bl_stat,ngrid, largest, .FALSE.)
1093
1094	!!! fin nrlmd le 10/04/2012
1095
1096	if (prt_level.ge.10) then
1097	ig=igout
1098	do l=1,nlay
1099	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
1100	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
1101	enddo
1102	endif
1103
1104	! print*,'avant calcul ale et alp'
1105	!calcul de ALE et ALP pour la convection
1106	Alp_bl(:)=0.
1107	Ale_bl(:)=0.
1108	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
1109	do l=1,nlay
1110	do ig=1,ngrid
1111	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
1112	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
1113	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
1114	enddo
1115	enddo
1116
1117	!test:calcul de la ponderation des couches pour KE
1118	!initialisations
1119
1120	fm_tot(:)=0.
1121	wght_th(:,:)=1.
1122	lalim_conv(:)=lalim(:)
1123
1124	do k=1,nlay
1125	do ig=1,ngrid
1126	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
1127	enddo
1128	enddo
1129	lfname='after calculation of Al[p/e]_bl'
1130	lvarname = 'fm'
1131	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
1132	lvarname = 'rhobarz'
1133	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
1134	lvarname = 'wth3'
1135	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
1136	lvarname = 'Alp_bl'
1137	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1138	lvarname = 'Ale_bl'
1139	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1140
1141	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
1142	! plus petit que 1.e-10, on prend wght_th=1.
1143	do k=1,nlay
1144	do ig=1,ngrid
1145	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1146	wght_th(ig,k)=alim_star(ig,k)
1147	endif
1148	enddo
1149	enddo
1150
1151	! print*,'apres wght_th'
1152	!test pour prolonger la convection
1153	do ig=1,ngrid
1154	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1155	if ((alim_star(ig,1).lt.1.e-10)) then
1156	lalim_conv(ig)=1
1157	wght_th(ig,1)=1.
1158	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1159	endif
1160	enddo
1161
1162	!------------------------------------------------------------------------
1163	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1164	! Integrale verticale de ALP.
1165	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1166	! couches
1167	!------------------------------------------------------------------------
1168
1169	alp_int(:)=0.
1170	dp_int(:)=0.
1171	do l=2,nlay
1172	do ig=1,ngrid
1173	if(l.LE.lmax(ig)) THEN
1174	zdp=pplay(ig,l-1)-pplay(ig,l)
1175	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1176	dp_int(ig)=dp_int(ig)+zdp
1177	endif
1178	enddo
1179	enddo
1180
1181	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1182	do ig=1,ngrid
1183	!valeur integree de alp_bl * 0.5:
1184	if (dp_int(ig)>0.) then
1185	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1186	endif
1187	enddo!
1188	endif
1189
1190
1191	! Facteur multiplicatif sur Alp_bl
1192	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1193
1194	lfname='thermcell_main last computations on Alp_bl'
1195	lvarname = 'pplay'
1196	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest*10.e4, .FALSE.)
1197	lvarname = 'alp_int'
1198	CALL check_var(lfname, lvarname, alp_int, ngrid, largest*10.e4, .FALSE.)
1199	lvarname = 'dp_int'
1200	CALL check_var(lfname, lvarname, dp_int, ngrid, largest, .FALSE.)
1201	lvarname = 'Alp_bl'
1202	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1203
1204	!------------------------------------------------------------------------
1205
1206
1207	!calcul du ratqscdiff
1208	if (prt_level.ge.1) print*,'14e OK convect8'
1209	var=0.
1210	vardiff=0.
1211	ratqsdiff(:,:)=0.
1212
1213	do l=1,nlay
1214	do ig=1,ngrid
1215	if (l<=lalim(ig)) then
1216	var=var+alim_star(ig,l)zqta(ig,l)1000.
1217	endif
1218	enddo
1219	enddo
1220
1221	if (prt_level.ge.1) print*,'14f OK convect8'
1222
1223	do l=1,nlay
1224	do ig=1,ngrid
1225	if (l<=lalim(ig)) then
1226	zf=fraca(ig,l)
1227	zf2=zf/(1.-zf)
1228	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1229	endif
1230	enddo
1231	enddo
1232
1233	if (prt_level.ge.1) print*,'14g OK convect8'
1234	do l=1,nlay
1235	do ig=1,ngrid
1236	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1237	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1238	enddo
1239	enddo
1240	!--------------------------------------------------------------------
1241	!
1242	!ecriture des fichiers sortie
1243	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1244
1245	#ifdef wrgrads_thermcell
1246	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1247	#include "thermcell_out3d.h"
1248	#endif
1249
1250	endif
1251
1252	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1253
1254	lfname='before leaving thermcell_main'
1255	lvarname = 'pt'
1256	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1257	lvarname = 'pdtadj'
1258	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1259	lvarname = 'Alp_bl'
1260	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1261	lvarname = 'Ale_bl'
1262	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1263
1264
1265	return
1266	end
1267
1268	!-----------------------------------------------------------------------------
1269
1270	subroutine test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1271	IMPLICIT NONE
1272	#include "iniprint.h"
1273
1274	integer i, k, ngrid,nlay
1275	real pplev(ngrid,nlay+1),pplay(ngrid,nlay)
1276	real ztv(ngrid,nlay)
1277	real po(ngrid,nlay)
1278	real ztva(ngrid,nlay)
1279	real zqla(ngrid,nlay)
1280	real f_star(ngrid,nlay)
1281	real zw2(ngrid,nlay)
1282	integer long(ngrid)
1283	real seuil
1284	character*21 comment
1285
1286	if (prt_level.ge.1) THEN
1287	print*,'WARNING !!! TEST ',comment
1288	endif
1289	return
1290
1291	! test sur la hauteur des thermiques ...
1292	do i=1,ngrid
1293	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1294	if (prt_level.ge.10) then
1295	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1296	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1297	do k=1,nlay
1298	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)
1299	enddo
1300	endif
1301	enddo
1302
1303
1304	return
1305	end
1306
1307	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1308	! On transporte pbl_tke pour donner therm_tke
1309	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1310	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1311	& rg,pplev,therm_tke_max)
1312	implicit none
1313
1314	#include "iniprint.h"
1315	!=======================================================================
1316	!
1317	! Calcul du transport verticale dans la couche limite en presence
1318	! de "thermiques" explicitement representes
1319	! calcul du dq/dt une fois qu'on connait les ascendances
1320	!
1321	!=======================================================================
1322
1323	integer ngrid,nlay,nsrf
1324
1325	real ptimestep
1326	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1327	real entr0(ngrid,nlay),rg
1328	real therm_tke_max(ngrid,nlay)
1329	real detr0(ngrid,nlay)
1330
1331
1332	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1333	real entr(ngrid,nlay)
1334	real q(ngrid,nlay)
1335	integer lev_out ! niveau pour les print
1336
1337	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1338
1339	real zzm
1340
1341	integer ig,k
1342	integer isrf
1343
1344
1345	lev_out=0
1346
1347
1348	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1349
1350	! calcul du detrainement
1351	do k=1,nlay
1352	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1353	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1354	enddo
1355
1356
1357	! Decalage vertical des entrainements et detrainements.
1358	masse(:,1)=0.5*masse0(:,1)
1359	entr(:,1)=0.5*entr0(:,1)
1360	detr(:,1)=0.5*detr0(:,1)
1361	fm(:,1)=0.
1362	do k=1,nlay-1
1363	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1364	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1365	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1366	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1367	enddo
1368	fm(:,nlay+1)=0.
1369
1370	!!! nrlmd le 16/09/2010
1371	! calcul de la valeur dans les ascendances
1372	! do ig=1,ngrid
1373	! qa(ig,1)=q(ig,1)
1374	! enddo
1375	!!!
1376
1377	!do isrf=1,nsrf
1378
1379	! q(:,:)=therm_tke(:,:,isrf)
1380	q(:,:)=therm_tke_max(:,:)
1381	!!! nrlmd le 16/09/2010
1382	do ig=1,ngrid
1383	qa(ig,1)=q(ig,1)
1384	enddo
1385	!!!
1386
1387	if (1==1) then
1388	do k=2,nlay
1389	do ig=1,ngrid
1390	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1391	& 1.e-5*masse(ig,k)) then
1392	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1393	& /(fm(ig,k+1)+detr(ig,k))
1394	else
1395	qa(ig,k)=q(ig,k)
1396	endif
1397	if (qa(ig,k).lt.0.) then
1398	! print*,'qa<0!!!'
1399	endif
1400	if (q(ig,k).lt.0.) then
1401	! print*,'q<0!!!'
1402	endif
1403	enddo
1404	enddo
1405
1406	! Calcul du flux subsident
1407
1408	do k=2,nlay
1409	do ig=1,ngrid
1410	wqd(ig,k)=fm(ig,k)*q(ig,k)
1411	if (wqd(ig,k).lt.0.) then
1412	! print*,'wqd<0!!!'
1413	endif
1414	enddo
1415	enddo
1416	do ig=1,ngrid
1417	wqd(ig,1)=0.
1418	wqd(ig,nlay+1)=0.
1419	enddo
1420
1421	! Calcul des tendances
1422	do k=1,nlay
1423	do ig=1,ngrid
1424	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1425	& -wqd(ig,k)+wqd(ig,k+1)) &
1426	& *ptimestep/masse(ig,k)
1427	enddo
1428	enddo
1429
1430	endif
1431
1432	therm_tke_max(:,:)=q(:,:)
1433
1434	return
1435	!!! fin nrlmd le 10/04/2012
1436	end
1437

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

Download in other formats:

Original Format