Context Navigation

thermcell_main.F90 @ 113

Last change on this file since 113 was 109, checked in by lfita, 11 years ago
Checking NaNs? from alim_star
File size: 51.0 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='thermcell_main 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	lvarname = 'pplev'
350	CALL check_var3D(lfname, lvarname, pplev, ngrid, nlay, largest, .FALSE.)
351
352	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
353	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
354
355	lfname='thermcell_main after thermcell_env'
356	lvarname = 'pt'
357	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
358	lvarname = 'pdtadj'
359	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
360	lvarname = 'pplev'
361	CALL check_var3D(lfname, lvarname, pplev, ngrid, nlay, largest, .FALSE.)
362
363	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
364
365	!------------------------------------------------------------------------
366	! --------------------
367	!
368	!
369	! + + + + + + + + + + +
370	!
371	!
372	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
373	! wh,wt,wo ...
374	!
375	! + + + + + + + + + + + zh,zu,zv,zo,rho
376	!
377	!
378	! -------------------- zlev(1)
379	! \\\\\\\\\\\\\\\\\\\\
380	!
381	!
382
383	!-----------------------------------------------------------------------
384	! Calcul des altitudes des couches
385	!-----------------------------------------------------------------------
386
387	do l=2,nlay
388	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
389	enddo
390	zlev(:,1)=0.
391	zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
392	do l=1,nlay
393	zlay(:,l)=pphi(:,l)/RG
394	enddo
395	!calcul de l epaisseur des couches
396	do l=1,nlay
397	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
398	enddo
399
400	! print*,'2 OK convect8'
401	!-----------------------------------------------------------------------
402	! Calcul des densites
403	!-----------------------------------------------------------------------
404
405	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
406
407	if (prt_level.ge.10)write(lunout,*) &
408	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
409	rhobarz(:,1)=rho(:,1)
410
411	do l=2,nlay
412	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
413	enddo
414
415	!calcul de la masse
416	do l=1,nlay
417	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
418	enddo
419
420	lfname='thermcell_main after initialization'
421	lvarname = 'rhobarz'
422	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
423	lvarname = 'rho'
424	CALL check_var3D(lfname, lvarname, rho, ngrid, nlay, largest, .FALSE.)
425	lvarname = 'zpspsk'
426	CALL check_var3D(lfname, lvarname, zpspsk, ngrid, nlay, largest, .FALSE.)
427	lvarname = 'pplay'
428	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest, .FALSE.)
429	lvarname = 'zw2'
430	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
431	lvarname = 'alim_star'
432	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
433
434	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
435
436	!------------------------------------------------------------------
437	!
438	! /\|\
439	! -------- \| F_k+1 -------
440	! ----> D_k
441	! /\|\ <---- E_k , A_k
442	! -------- \| F_k ---------
443	! ----> D_k-1
444	! <---- E_k-1 , A_k-1
445	!
446	!
447	!
448	!
449	!
450	! ---------------------------
451	!
452	! ----- F_lmax+1=0 ---------- \
453	! lmax (zmax) \|
454	! --------------------------- \|
455	! \|
456	! --------------------------- \|
457	! \|
458	! --------------------------- \|
459	! \|
460	! --------------------------- \|
461	! \|
462	! --------------------------- \|
463	! \| E
464	! --------------------------- \| D
465	! \|
466	! --------------------------- \|
467	! \|
468	! --------------------------- \ \|
469	! lalim \| \|
470	! --------------------------- \| \|
471	! \| \|
472	! --------------------------- \| \|
473	! \| A \|
474	! --------------------------- \| \|
475	! \| \|
476	! --------------------------- \| \|
477	! lmin (=1 pour le moment) \| \|
478	! ----- F_lmin=0 ------------ / /
479	!
480	! ---------------------------
481	! //////////////////////////
482	!
483	!
484	!=============================================================================
485	! Calculs initiaux ne faisant pas intervenir les changements de phase
486	!=============================================================================
487
488	!------------------------------------------------------------------
489	! 1. alim_star est le profil vertical de l'alimentation a la base du
490	! panache thermique, calcule a partir de la flotabilite de l'air sec
491	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
492	!------------------------------------------------------------------
493	!
494	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
495	lmin=1
496
497	!-----------------------------------------------------------------------------
498	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
499	! panache sec conservatif (e=d=0) alimente selon alim_star
500	! Il s'agit d'un calcul de type CAPE
501	! zmax_sec est utilise pour determiner la geometrie du thermique.
502	!------------------------------------------------------------------------------
503	!---------------------------------------------------------------------------------
504	!calcul du melange et des variables dans le thermique
505	!--------------------------------------------------------------------------------
506	!
507	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
508	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
509
510	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
511	! de la variable iflag_thermals
512	lfname='thermcell_main before plume'
513	lvarname = 'alim_star'
514	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
515
516	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
517	if (iflag_thermals_ed<=9) then
518	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
519	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
520	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
521	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
522	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
523	& ,lev_out,lunout1,igout)
524
525	elseif (iflag_thermals_ed>9) then
526	! print*,'THERM RIO et al 2010, version d Arnaud'
527	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
528	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
529	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
530	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
531	& ,lev_out,lunout1,igout)
532
533	endif
534	lfname='thermcell_main after thermcell_plume'
535	lvarname = 'pt'
536	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
537	lvarname = 'pdtadj'
538	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
539	lvarname = 'zw2'
540	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
541	lvarname = 'alim_star'
542	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
543
544	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
545
546	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
547	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
548
549	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
550	if (prt_level.ge.10) then
551	write(lunout1,*) 'Dans thermcell_main 2'
552	write(lunout1,*) 'lmin ',lmin(igout)
553	write(lunout1,*) 'lalim ',lalim(igout)
554	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
555	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
556	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
557	endif
558
559	!-------------------------------------------------------------------------------
560	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
561	!-------------------------------------------------------------------------------
562	!
563	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
564	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
565	! Attention, w2 est transforme en sa racine carree dans cette routine
566	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
567	wmax_tmp=0.
568	do l=1,nlay
569	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
570	enddo
571	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
572
573	lfname='thermcell_main after thermcell_height'
574	lvarname = 'zw2'
575	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
576
577	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
578	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
579	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
580	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
581
582	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
583
584	!-------------------------------------------------------------------------------
585	! Fermeture,determination de f
586	!-------------------------------------------------------------------------------
587	!
588	!
589	!! write(lunout,*)'THERM NOUVEAU XXXXX'
590	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
591	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
592
593	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
594	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
595
596	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
597	if (prt_level.ge.10) then
598	write(lunout1,*) 'Dans thermcell_main 1b'
599	write(lunout1,*) 'lmin ',lmin(igout)
600	write(lunout1,*) 'lalim ',lalim(igout)
601	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
602	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
603	& ,l=1,lalim(igout)+4)
604	endif
605
606	lfname='thermcell_main after thermcell_dry'
607	lvarname = 'alim_star'
608	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
609
610
611
612	! Choix de la fonction d'alimentation utilisee pour la fermeture.
613	! Apparemment sans importance
614	alim_star_clos(:,:)=alim_star(:,:)
615	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
616
617	! Appel avec la version seche
618	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
619	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
620
621	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
622	! Appel avec les zmax et wmax tenant compte de la condensation
623	! Semble moins bien marcher
624	! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
625	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
626	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
627
628	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
629
630	if (tau_thermals>1.) then
631	lambda=exp(-ptimestep/tau_thermals)
632	f0=(1.-lambda)f+lambdaf0
633	else
634	f0=f
635	endif
636
637	! Test valable seulement en 1D mais pas genant
638	if (.not. (f0(1).ge.0.) ) then
639	abort_message = '.not. (f0(1).ge.0.)'
640	CALL abort_gcm (modname,abort_message,1)
641	endif
642
643	!-------------------------------------------------------------------------------
644	!deduction des flux
645	!-------------------------------------------------------------------------------
646	lfname='thermcell_main before flux'
647	lvarname = 'fm'
648	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
649	lvarname = 'entr'
650	CALL check_var3D(lfname, lvarname, entr, ngrid, nlay, largest, .FALSE.)
651	lvarname = 'detr'
652	CALL check_var3D(lfname, lvarname, detr, ngrid, nlay, largest, .FALSE.)
653	lvarname = 'zqla'
654	CALL check_var3D(lfname, lvarname, zqla, ngrid, nlay, largest, .FALSE.)
655	lvarname = 'zw2'
656	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
657	lvarname = 'alim_star'
658	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
659
660	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
661	& lalim,lmax,alim_star, &
662	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
663	& detr,zqla,lev_out,lunout1,igout)
664	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
665	lfname='thermcell_main after flux'
666	lvarname = 'fm'
667	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
668	lvarname = 'zw2'
669	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
670	lvarname = 'alim_star'
671	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
672
673	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
674	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
675	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
676
677	!------------------------------------------------------------------
678	! On ne prend pas directement les profils issus des calculs precedents
679	! mais on s'autorise genereusement une relaxation vers ceci avec
680	! une constante de temps tau_thermals (typiquement 1800s).
681	!------------------------------------------------------------------
682
683	if (tau_thermals>1.) then
684	lambda=exp(-ptimestep/tau_thermals)
685	fm0=(1.-lambda)fm+lambdafm0
686	entr0=(1.-lambda)entr+lambdaentr0
687	detr0=(1.-lambda)detr+lambdadetr0
688	else
689	fm0=fm
690	entr0=entr
691	detr0=detr
692	endif
693
694	!c------------------------------------------------------------------
695	! calcul du transport vertical
696	!------------------------------------------------------------------
697
698	lfname='thermcell_main before transport_vertical'
699	lvarname = 'zdthladj'
700	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
701
702	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
703	& zthl,zdthladj,zta,lev_out)
704	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
705	& po,pdoadj,zoa,lev_out)
706	lfname='thermcell_main after transport_vertical'
707	lvarname = 'zdthladj'
708	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
709	lvarname = 'masse'
710	CALL check_var3D(lfname, lvarname, masse, ngrid, nlay, largest, .FALSE.)
711
712	lfname='thermcell_main before fraction ascendance'
713	lvarname = 'fm'
714	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
715	lvarname = 'zw2'
716	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
717
718	!------------------------------------------------------------------
719	! Calcul de la fraction de l'ascendance
720	!------------------------------------------------------------------
721	do ig=1,ngrid
722	fraca(ig,1)=0.
723	fraca(ig,nlay+1)=0.
724	enddo
725	do l=2,nlay
726	do ig=1,ngrid
727	if (zw2(ig,l).gt.1.e-10) then
728	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
729	else
730	fraca(ig,l)=0.
731	endif
732	enddo
733	enddo
734
735	!------------------------------------------------------------------
736	! calcul du transport vertical du moment horizontal
737	!------------------------------------------------------------------
738	lfname='before thermcell_dv2'
739	lvarname = 'pt'
740	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
741	lvarname = 'pdtadj'
742	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
743	lvarname = 'zw2'
744	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
745
746	!IM 090508
747	if (dvdq == 0 ) then
748
749	! Calcul du transport de V tenant compte d'echange par gradient
750	! de pression horizontal avec l'environnement
751
752	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
753	! & ,fraca*dvdq,zmax &
754	& ,fraca,zmax &
755	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
756
757	else
758
759	! calcul purement conservatif pour le transport de V
760	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
761	& ,zu,pduadj,zua,lev_out)
762	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
763	& ,zv,pdvadj,zva,lev_out)
764
765	endif
766
767	! print*,'13 OK convect8'
768	do l=1,nlay
769	do ig=1,ngrid
770	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
771	enddo
772	enddo
773	lfname='after thermcell_dv2'
774	lvarname = 'pt'
775	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
776	lvarname = 'pdtadj'
777	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
778	lvarname = 'zdthladj'
779	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
780	lvarname = 'zpspsk'
781	CALL check_var3D(lfname, lvarname, zpspsk, ngrid, nlay, largest, .FALSE.)
782
783	if (prt_level.ge.1) print*,'14 OK convect8'
784	!------------------------------------------------------------------
785	! Calculs de diagnostiques pour les sorties
786	!------------------------------------------------------------------
787	!calcul de fraca pour les sorties
788
789	if (sorties) then
790	if (prt_level.ge.1) print*,'14a OK convect8'
791	! calcul du niveau de condensation
792	! initialisation
793	do ig=1,ngrid
794	nivcon(ig)=0
795	zcon(ig)=0.
796	enddo
797	!nouveau calcul
798	do ig=1,ngrid
799	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
800	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
801	enddo
802	!IM do k=1,nlay
803	do k=1,nlay-1
804	do ig=1,ngrid
805	if ((pcon(ig).le.pplay(ig,k)) &
806	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
807	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
808	endif
809	enddo
810	enddo
811	!IM
812	ierr=0
813	do ig=1,ngrid
814	if (pcon(ig).le.pplay(ig,nlay)) then
815	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
816	ierr=1
817	endif
818	enddo
819	if (ierr==1) then
820	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
821	CALL abort_gcm (modname,abort_message,1)
822	endif
823
824	if (prt_level.ge.1) print*,'14b OK convect8'
825	do k=nlay,1,-1
826	do ig=1,ngrid
827	if (zqla(ig,k).gt.1e-10) then
828	nivcon(ig)=k
829	zcon(ig)=zlev(ig,k)
830	endif
831	enddo
832	enddo
833	if (prt_level.ge.1) print*,'14c OK convect8'
834	!calcul des moments
835	!initialisation
836	do l=1,nlay
837	do ig=1,ngrid
838	q2(ig,l)=0.
839	wth2(ig,l)=0.
840	wth3(ig,l)=0.
841	ratqscth(ig,l)=0.
842	ratqsdiff(ig,l)=0.
843	enddo
844	enddo
845	if (prt_level.ge.1) print*,'14d OK convect8'
846	if (prt_level.ge.10)write(lunout,*) &
847	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
848	do l=1,nlay
849	do ig=1,ngrid
850	zf=fraca(ig,l)
851	zf2=zf/(1.-zf)
852	!
853	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
854	if(zw2(ig,l).gt.1.e-10) then
855	wth2(ig,l)=zf2(zw2(ig,l))*2
856	else
857	wth2(ig,l)=0.
858	endif
859	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
860	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
861	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
862	!test: on calcul q2/po=ratqsc
863	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
864	enddo
865	enddo
866	lfname='thermcell_main calculation of wth3'
867	lvarname = 'wth3'
868	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
869	lvarname = 'fraca'
870	CALL check_var3D(lfname, lvarname, fraca, ngrid, nlay, largest, .FALSE.)
871	lvarname = 'zw2'
872	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay, largest, .FALSE.)
873	lvarname = '1-fraca'
874	CALL check_var3D(lfname, lvarname, 1./(1.-fraca), ngrid, nlay, largest, .FALSE.)
875	!calcul des flux: q, thetal et thetav
876	do l=1,nlay
877	do ig=1,ngrid
878	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
879	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
880	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
881	enddo
882	enddo
883	!
884
885	!!! nrlmd le 10/04/2012
886
887	!------------Test sur le LCL des thermiques
888	do ig=1,ngrid
889	ok_lcl(ig)=.false.
890	if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
891	enddo
892
893	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
894	do l=1,nlay-1
895	do ig=1,ngrid
896	if (ok_lcl(ig)) then
897	if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
898	klcl(ig)=l
899	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
900	endif
901	endif
902	enddo
903	enddo
904
905	!------------Hauteur des thermiques
906	!!jyg le 27/04/2012
907	!! do ig =1,ngrid
908	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
909	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
910	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
911	!! zmax(ig)=pphi(ig,lmax(ig))/rg
912	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
913	!! enddo
914	do ig =1,ngrid
915	zmax(ig)=pphi(ig,lmax(ig))/rg
916	if (ok_lcl(ig)) then
917	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
918	& -rhobarz(ig,klcl(ig)))*interp(ig)
919	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
920	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
921	else
922	rhobarz0(ig)=0.
923	zlcl(ig)=zmax(ig)
924	endif
925	enddo
926	!!jyg fin
927
928	!------------Calcul des propriétés du thermique au LCL
929	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
930
931	!-----Initialisation de la TKE moyenne
932	do l=1,nlay
933	do ig=1,ngrid
934	pbl_tke_max(ig,l)=0.
935	enddo
936	enddo
937
938	!-----Calcul de la TKE moyenne
939	do nsrf=1,nbsrf
940	do l=1,nlay
941	do ig=1,ngrid
942	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
943	enddo
944	enddo
945	enddo
946
947	!-----Initialisations des TKE dans et hors des thermiques
948	do l=1,nlay
949	do ig=1,ngrid
950	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
951	env_tke_max(ig,l)=pbl_tke_max(ig,l)
952	enddo
953	enddo
954
955	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
956	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
957	& rg,pplev,therm_tke_max)
958	! print *,' thermcell_tke_transport -> ' !!jyg
959
960	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
961	do l=1,nlay
962	do ig=1,ngrid
963	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
964	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
965	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
966	enddo
967	enddo
968	! print *,' apres w_ls = ' !!jyg
969
970	do ig=1,ngrid
971	if (ok_lcl(ig)) then
972	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
973	& -fraca(ig,klcl(ig)))*interp(ig)
974	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
975	& -zw2(ig,klcl(ig)))*interp(ig)
976	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
977	& -w_ls(ig,klcl(ig)))*interp(ig)
978	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
979	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
980	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
981	& -env_tke_max(ig,klcl(ig)))*interp(ig)
982	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
983	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
984	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
985	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
986	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
987	else
988	fraca0(ig)=0.
989	w0(ig)=0.
990	!!jyg le 27/04/2012
991	!! zlcl(ig)=0.
992	!!
993	endif
994	enddo
995
996	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
997	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
998
999	!------------Triggering------------------
1000	IF (iflag_trig_bl.ge.1) THEN
1001
1002	!-----Initialisations
1003	depth(:)=0.
1004	n2(:)=0.
1005	s2(:)=0.
1006	s_max(:)=0.
1007
1008	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
1009	do ig=1,ngrid
1010	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
1011	depth(ig)=zmax_moy(ig)-zlcl(ig)
1012	hmin(ig)=hmincoef*zlcl(ig)
1013	if (depth(ig).ge.10.) then
1014	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
1015	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
1016	!!
1017	!!jyg le 27/04/2012
1018	!! s_max(ig)=s2(ig)*log(n2(ig))
1019	!! if (n2(ig) .lt. 1) s_max(ig)=0.
1020	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
1021	!!fin jyg
1022	else
1023	s2(ig)=0.
1024	n2(ig)=0.
1025	s_max(ig)=0.
1026	endif
1027	enddo
1028	! print *,'avant Calcul de Wmax ' !!jyg
1029
1030	!-----Calcul de Wmax et ALE_BL_STAT associée
1031	!!jyg le 30/04/2012
1032	!! do ig=1,ngrid
1033	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
1034	!! 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))))
1035	!! ale_bl_stat(ig)=0.5w_max(ig)*2
1036	!! else
1037	!! w_max(ig)=0.
1038	!! ale_bl_stat(ig)=0.
1039	!! endif
1040	!! enddo
1041	susqr2pi=susqrt(2.Rpi)
1042	Reuler=exp(1.)
1043	do ig=1,ngrid
1044	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
1045	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
1046	ale_bl_stat(ig)=0.5w_max(ig)*2
1047	else
1048	w_max(ig)=0.
1049	ale_bl_stat(ig)=0.
1050	endif
1051	enddo
1052
1053	ENDIF ! iflag_trig_bl
1054	! print *,'ENDIF iflag_trig_bl' !!jyg
1055	lfname='thermcell_main before closure'
1056	lvarname = 'pt'
1057	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1058	lvarname = 'pdtadj'
1059	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1060	lvarname = 'alp_bl_det'
1061	CALL check_var(lfname, lvarname, alp_bl_det,ngrid, largest, .FALSE.)
1062	lvarname = 'rhobarz0'
1063	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1064	lvarname = 'fraca0'
1065	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1066	lvarname = 'w_conv'
1067	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1068	lvarname = 'w0'
1069	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1070	lvarname = 'alp_bl_fluct_m'
1071	CALL check_var(lfname, lvarname, alp_bl_fluct_m,ngrid, largest, .FALSE.)
1072	lvarname = 'alp_bl_fluct_tke'
1073	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,ngrid, largest, .FALSE.)
1074	lvarname = 'therm_tke_max0'
1075	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1076	lvarname = 'env_tke_max0'
1077	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1078	lvarname = 'pbl_tke_max0'
1079	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1080	lvarname = 'alp_bl_conv'
1081	CALL check_var(lfname, lvarname, alp_bl_conv,ngrid, largest, .FALSE.)
1082	lvarname = 'alp_bl_stat'
1083	CALL check_var(lfname, lvarname, alp_bl_stat,ngrid, largest, .FALSE.)
1084	!------------Closure------------------
1085
1086	IF (iflag_clos_bl.ge.1) THEN
1087
1088	!-----Calcul de ALP_BL_STAT
1089	do ig=1,ngrid
1090	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
1091	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
1092	& (w0(ig)**2)
1093	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
1094	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
1095	if (iflag_clos_bl.ge.2) then
1096	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
1097	& (w0(ig)**2)
1098	else
1099	alp_bl_conv(ig)=0.
1100	endif
1101	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
1102	enddo
1103
1104	!-----Sécurité ALP infinie
1105	do ig=1,ngrid
1106	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
1107	enddo
1108
1109	ENDIF ! (iflag_clos_bl.ge.1)
1110	lfname='after closure'
1111	lvarname = 'pt'
1112	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1113	lvarname = 'pdtadj'
1114	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1115	lvarname = 'alp_bl_det'
1116	CALL check_var(lfname, lvarname, alp_bl_det,ngrid, largest, .FALSE.)
1117	lvarname = 'rhobarz0'
1118	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1119	lvarname = 'fraca0'
1120	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1121	lvarname = 'w_conv'
1122	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1123	lvarname = 'w0'
1124	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1125	lvarname = 'alp_bl_fluct_m'
1126	CALL check_var(lfname, lvarname, alp_bl_fluct_m,ngrid, largest, .FALSE.)
1127	lvarname = 'alp_bl_fluct_tke'
1128	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,ngrid, largest, .FALSE.)
1129	lvarname = 'therm_tke_max0'
1130	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1131	lvarname = 'env_tke_max0'
1132	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1133	lvarname = 'pbl_tke_max0'
1134	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1135	lvarname = 'alp_bl_conv'
1136	CALL check_var(lfname, lvarname, alp_bl_conv,ngrid, largest, .FALSE.)
1137	lvarname = 'alp_bl_stat'
1138	CALL check_var(lfname, lvarname, alp_bl_stat,ngrid, largest, .FALSE.)
1139
1140	!!! fin nrlmd le 10/04/2012
1141
1142	if (prt_level.ge.10) then
1143	ig=igout
1144	do l=1,nlay
1145	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
1146	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
1147	enddo
1148	endif
1149
1150	! print*,'avant calcul ale et alp'
1151	!calcul de ALE et ALP pour la convection
1152	Alp_bl(:)=0.
1153	Ale_bl(:)=0.
1154	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
1155	do l=1,nlay
1156	do ig=1,ngrid
1157	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
1158	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
1159	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
1160	enddo
1161	enddo
1162
1163	!test:calcul de la ponderation des couches pour KE
1164	!initialisations
1165
1166	fm_tot(:)=0.
1167	wght_th(:,:)=1.
1168	lalim_conv(:)=lalim(:)
1169
1170	do k=1,nlay
1171	do ig=1,ngrid
1172	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
1173	enddo
1174	enddo
1175	lfname='after calculation of Al[p/e]_bl'
1176	lvarname = 'fm'
1177	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
1178	lvarname = 'rhobarz'
1179	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
1180	lvarname = 'wth3'
1181	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
1182	lvarname = 'Alp_bl'
1183	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1184	lvarname = 'Ale_bl'
1185	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1186
1187	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
1188	! plus petit que 1.e-10, on prend wght_th=1.
1189	do k=1,nlay
1190	do ig=1,ngrid
1191	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1192	wght_th(ig,k)=alim_star(ig,k)
1193	endif
1194	enddo
1195	enddo
1196
1197	! print*,'apres wght_th'
1198	!test pour prolonger la convection
1199	do ig=1,ngrid
1200	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1201	if ((alim_star(ig,1).lt.1.e-10)) then
1202	lalim_conv(ig)=1
1203	wght_th(ig,1)=1.
1204	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1205	endif
1206	enddo
1207
1208	!------------------------------------------------------------------------
1209	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1210	! Integrale verticale de ALP.
1211	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1212	! couches
1213	!------------------------------------------------------------------------
1214
1215	alp_int(:)=0.
1216	dp_int(:)=0.
1217	do l=2,nlay
1218	do ig=1,ngrid
1219	if(l.LE.lmax(ig)) THEN
1220	zdp=pplay(ig,l-1)-pplay(ig,l)
1221	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1222	dp_int(ig)=dp_int(ig)+zdp
1223	endif
1224	enddo
1225	enddo
1226
1227	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1228	do ig=1,ngrid
1229	!valeur integree de alp_bl * 0.5:
1230	if (dp_int(ig)>0.) then
1231	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1232	endif
1233	enddo!
1234	endif
1235
1236
1237	! Facteur multiplicatif sur Alp_bl
1238	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1239
1240	lfname='thermcell_main last computations on Alp_bl'
1241	lvarname = 'pplay'
1242	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest*10.e4, .FALSE.)
1243	lvarname = 'alp_int'
1244	CALL check_var(lfname, lvarname, alp_int, ngrid, largest*10.e4, .FALSE.)
1245	lvarname = 'dp_int'
1246	CALL check_var(lfname, lvarname, dp_int, ngrid, largest, .FALSE.)
1247	lvarname = 'Alp_bl'
1248	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1249
1250	!------------------------------------------------------------------------
1251
1252
1253	!calcul du ratqscdiff
1254	if (prt_level.ge.1) print*,'14e OK convect8'
1255	var=0.
1256	vardiff=0.
1257	ratqsdiff(:,:)=0.
1258
1259	do l=1,nlay
1260	do ig=1,ngrid
1261	if (l<=lalim(ig)) then
1262	var=var+alim_star(ig,l)zqta(ig,l)1000.
1263	endif
1264	enddo
1265	enddo
1266
1267	if (prt_level.ge.1) print*,'14f OK convect8'
1268
1269	do l=1,nlay
1270	do ig=1,ngrid
1271	if (l<=lalim(ig)) then
1272	zf=fraca(ig,l)
1273	zf2=zf/(1.-zf)
1274	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1275	endif
1276	enddo
1277	enddo
1278
1279	if (prt_level.ge.1) print*,'14g OK convect8'
1280	do l=1,nlay
1281	do ig=1,ngrid
1282	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1283	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1284	enddo
1285	enddo
1286	!--------------------------------------------------------------------
1287	!
1288	!ecriture des fichiers sortie
1289	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1290
1291	#ifdef wrgrads_thermcell
1292	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1293	#include "thermcell_out3d.h"
1294	#endif
1295
1296	endif
1297
1298	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1299
1300	lfname='before leaving thermcell_main'
1301	lvarname = 'pt'
1302	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1303	lvarname = 'pdtadj'
1304	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1305	lvarname = 'Alp_bl'
1306	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1307	lvarname = 'Ale_bl'
1308	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1309
1310
1311	return
1312	end
1313
1314	!-----------------------------------------------------------------------------
1315
1316	subroutine test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1317	IMPLICIT NONE
1318	#include "iniprint.h"
1319
1320	integer i, k, ngrid,nlay
1321	real pplev(ngrid,nlay+1),pplay(ngrid,nlay)
1322	real ztv(ngrid,nlay)
1323	real po(ngrid,nlay)
1324	real ztva(ngrid,nlay)
1325	real zqla(ngrid,nlay)
1326	real f_star(ngrid,nlay)
1327	real zw2(ngrid,nlay)
1328	integer long(ngrid)
1329	real seuil
1330	character*21 comment
1331
1332	if (prt_level.ge.1) THEN
1333	print*,'WARNING !!! TEST ',comment
1334	endif
1335	return
1336
1337	! test sur la hauteur des thermiques ...
1338	do i=1,ngrid
1339	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1340	if (prt_level.ge.10) then
1341	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1342	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1343	do k=1,nlay
1344	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)
1345	enddo
1346	endif
1347	enddo
1348
1349
1350	return
1351	end
1352
1353	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1354	! On transporte pbl_tke pour donner therm_tke
1355	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1356	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1357	& rg,pplev,therm_tke_max)
1358	implicit none
1359
1360	#include "iniprint.h"
1361	!=======================================================================
1362	!
1363	! Calcul du transport verticale dans la couche limite en presence
1364	! de "thermiques" explicitement representes
1365	! calcul du dq/dt une fois qu'on connait les ascendances
1366	!
1367	!=======================================================================
1368
1369	integer ngrid,nlay,nsrf
1370
1371	real ptimestep
1372	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1373	real entr0(ngrid,nlay),rg
1374	real therm_tke_max(ngrid,nlay)
1375	real detr0(ngrid,nlay)
1376
1377
1378	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1379	real entr(ngrid,nlay)
1380	real q(ngrid,nlay)
1381	integer lev_out ! niveau pour les print
1382
1383	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1384
1385	real zzm
1386
1387	integer ig,k
1388	integer isrf
1389
1390
1391	lev_out=0
1392
1393
1394	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1395
1396	! calcul du detrainement
1397	do k=1,nlay
1398	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1399	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1400	enddo
1401
1402
1403	! Decalage vertical des entrainements et detrainements.
1404	masse(:,1)=0.5*masse0(:,1)
1405	entr(:,1)=0.5*entr0(:,1)
1406	detr(:,1)=0.5*detr0(:,1)
1407	fm(:,1)=0.
1408	do k=1,nlay-1
1409	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1410	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1411	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1412	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1413	enddo
1414	fm(:,nlay+1)=0.
1415
1416	!!! nrlmd le 16/09/2010
1417	! calcul de la valeur dans les ascendances
1418	! do ig=1,ngrid
1419	! qa(ig,1)=q(ig,1)
1420	! enddo
1421	!!!
1422
1423	!do isrf=1,nsrf
1424
1425	! q(:,:)=therm_tke(:,:,isrf)
1426	q(:,:)=therm_tke_max(:,:)
1427	!!! nrlmd le 16/09/2010
1428	do ig=1,ngrid
1429	qa(ig,1)=q(ig,1)
1430	enddo
1431	!!!
1432
1433	if (1==1) then
1434	do k=2,nlay
1435	do ig=1,ngrid
1436	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1437	& 1.e-5*masse(ig,k)) then
1438	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1439	& /(fm(ig,k+1)+detr(ig,k))
1440	else
1441	qa(ig,k)=q(ig,k)
1442	endif
1443	if (qa(ig,k).lt.0.) then
1444	! print*,'qa<0!!!'
1445	endif
1446	if (q(ig,k).lt.0.) then
1447	! print*,'q<0!!!'
1448	endif
1449	enddo
1450	enddo
1451
1452	! Calcul du flux subsident
1453
1454	do k=2,nlay
1455	do ig=1,ngrid
1456	wqd(ig,k)=fm(ig,k)*q(ig,k)
1457	if (wqd(ig,k).lt.0.) then
1458	! print*,'wqd<0!!!'
1459	endif
1460	enddo
1461	enddo
1462	do ig=1,ngrid
1463	wqd(ig,1)=0.
1464	wqd(ig,nlay+1)=0.
1465	enddo
1466
1467	! Calcul des tendances
1468	do k=1,nlay
1469	do ig=1,ngrid
1470	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1471	& -wqd(ig,k)+wqd(ig,k+1)) &
1472	& *ptimestep/masse(ig,k)
1473	enddo
1474	enddo
1475
1476	endif
1477
1478	therm_tke_max(:,:)=q(:,:)
1479
1480	return
1481	!!! fin nrlmd le 10/04/2012
1482	end
1483

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

Download in other formats:

Original Format