Context Navigation

thermcell_main.F90 @ 89

Last change on this file since 89 was 89, checked in by lfita, 10 years ago
Fixing bad dimension size in check_var3D Initializing 'pbl_tke_max0' in thermcell_main
File size: 47.7 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(klon)
108	INTEGER lmax(klon),lmin(klon),lalim(klon)
109	INTEGER lmix(klon)
110	INTEGER lmix_bis(klon)
111	real linter(klon)
112	real zmix(klon)
113	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
114	! real fraca(klon,klev)
115
116	real zmax_sec(klon)
117	!on garde le zmax du pas de temps precedent
118	real zmax0(klon)
119	!FH/IM save zmax0
120
121	real lambda
122
123	real zlev(klon,klev+1),zlay(klon,klev)
124	real deltaz(klon,klev)
125	REAL zh(klon,klev)
126	real zthl(klon,klev),zdthladj(klon,klev)
127	REAL ztv(klon,klev)
128	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
129	real zl(klon,klev)
130	real zsortie(klon,klev)
131	real zva(klon,klev)
132	real zua(klon,klev)
133	real zoa(klon,klev)
134
135	real zta(klon,klev)
136	real zha(klon,klev)
137	real fraca(klon,klev+1)
138	real zf,zf2
139	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
140	real q2(klon,klev)
141	! FH probleme de dimensionnement avec l'allocation dynamique
142	! common/comtherm/thetath2,wth2
143	real wq(klon,klev)
144	real wthl(klon,klev)
145	real wthv(klon,klev)
146
147	real ratqscth(klon,klev)
148	real var
149	real vardiff
150	real ratqsdiff(klon,klev)
151
152	logical sorties
153	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
154	real zpspsk(klon,klev)
155
156	real wmax(klon)
157	real wmax_tmp(klon)
158	real wmax_sec(klon)
159	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
160	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
161
162	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
163	!niveau de condensation
164	integer nivcon(klon)
165	real zcon(klon)
166	REAL CHI
167	real zcon2(klon)
168	real pcon(klon)
169	real zqsat(klon,klev)
170	real zqsatth(klon,klev)
171
172	real f_star(klon,klev+1),entr_star(klon,klev)
173	real detr_star(klon,klev)
174	real alim_star_tot(klon)
175	real alim_star(klon,klev)
176	real alim_star_clos(klon,klev)
177	real f(klon), f0(klon)
178	!FH/IM save f0
179	real zlevinter(klon)
180	logical debut
181	real seuil
182	real csc(klon,klev)
183
184	!!! nrlmd le 10/04/2012
185
186	!------Entrées
187	real pbl_tke(klon,klev+1,nbsrf)
188	real pctsrf(klon,nbsrf)
189	real omega(klon,klev)
190	real airephy(klon)
191	!------Sorties
192	real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
193	real therm_tke_max0(klon),env_tke_max0(klon)
194	real n2(klon),s2(klon)
195	real ale_bl_stat(klon)
196	real therm_tke_max(klon,klev),env_tke_max(klon,klev)
197	real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
198	!------Local
199	integer nsrf
200	real rhobarz0(klon) ! Densité au LCL
201	logical ok_lcl(klon) ! Existence du LCL des thermiques
202	integer klcl(klon) ! Niveau du LCL
203	real interp(klon) ! 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(klon) ! Epaisseur moyenne du cumulus
218	real w_max(klon) ! Vitesse max statistique
219	real s_max(klon)
220	!--Closure
221	real pbl_tke_max(klon,klev) ! Profil de TKE moyenne
222	real pbl_tke_max0(klon) ! TKE moyenne au LCL
223	real w_ls(klon,klev) ! 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(klon)
234	real Alp_bl(klon)
235	real alp_int(klon),dp_int(klon),zdp
236	real ale_int(klon)
237	integer n_int(klon)
238	real fm_tot(klon)
239	real wght_th(klon,klev)
240	integer lalim_conv(klon)
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.klon) THEN
322	PRINT*
323	PRINT*,'STOP dans convadj'
324	PRINT*,'ngrid =',ngrid
325	PRINT*,'klon =',klon
326	ENDIF
327	!
328	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
329	do ig=1,klon
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, klon, klev, largest, .FALSE.)
347	lvarname = 'pdtadj'
348	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, 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, klon, klev, largest, .FALSE.)
356	lvarname = 'pdtadj'
357	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, 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(:,klev)-pphi(:,klev-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, klon, klev, largest, .FALSE.)
516	lvarname = 'pdtadj'
517	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, 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, klon, klev, 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, klon, klev, 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, klon, klev, largest, .FALSE.)
673
674	!------------------------------------------------------------------
675	! Calcul de la fraction de l'ascendance
676	!------------------------------------------------------------------
677	do ig=1,klon
678	fraca(ig,1)=0.
679	fraca(ig,nlay+1)=0.
680	enddo
681	do l=2,nlay
682	do ig=1,klon
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, klon, klev, largest, .FALSE.)
697	lvarname = 'pdtadj'
698	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, 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, klon, klev, largest, .FALSE.)
730	lvarname = 'pdtadj'
731	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, 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	!calcul des flux: q, thetal et thetav
821	do l=1,nlay
822	do ig=1,ngrid
823	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
824	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
825	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
826	enddo
827	enddo
828	!
829
830	!!! nrlmd le 10/04/2012
831
832	!------------Test sur le LCL des thermiques
833	do ig=1,ngrid
834	ok_lcl(ig)=.false.
835	if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
836	enddo
837
838	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
839	do l=1,nlay-1
840	do ig=1,ngrid
841	if (ok_lcl(ig)) then
842	if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
843	klcl(ig)=l
844	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
845	endif
846	endif
847	enddo
848	enddo
849
850	!------------Hauteur des thermiques
851	!!jyg le 27/04/2012
852	!! do ig =1,ngrid
853	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
854	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
855	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
856	!! zmax(ig)=pphi(ig,lmax(ig))/rg
857	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
858	!! enddo
859	do ig =1,ngrid
860	zmax(ig)=pphi(ig,lmax(ig))/rg
861	if (ok_lcl(ig)) then
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	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
866	else
867	rhobarz0(ig)=0.
868	zlcl(ig)=zmax(ig)
869	endif
870	enddo
871	!!jyg fin
872
873	!------------Calcul des propriétés du thermique au LCL
874	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
875
876	!-----Initialisation de la TKE moyenne
877	do l=1,nlay
878	do ig=1,ngrid
879	pbl_tke_max(ig,l)=0.
880	enddo
881	enddo
882
883	!-----Calcul de la TKE moyenne
884	do nsrf=1,nbsrf
885	do l=1,nlay
886	do ig=1,ngrid
887	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
888	enddo
889	enddo
890	enddo
891
892	!-----Initialisations des TKE dans et hors des thermiques
893	do l=1,nlay
894	do ig=1,ngrid
895	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
896	env_tke_max(ig,l)=pbl_tke_max(ig,l)
897	enddo
898	enddo
899
900	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
901	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
902	& rg,pplev,therm_tke_max)
903	! print *,' thermcell_tke_transport -> ' !!jyg
904
905	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
906	do l=1,nlay
907	do ig=1,ngrid
908	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
909	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
910	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
911	enddo
912	enddo
913	! print *,' apres w_ls = ' !!jyg
914
915	do ig=1,ngrid
916	if (ok_lcl(ig)) then
917	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
918	& -fraca(ig,klcl(ig)))*interp(ig)
919	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
920	& -zw2(ig,klcl(ig)))*interp(ig)
921	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
922	& -w_ls(ig,klcl(ig)))*interp(ig)
923	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
924	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
925	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
926	& -env_tke_max(ig,klcl(ig)))*interp(ig)
927	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
928	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
929	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
930	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
931	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
932	else
933	fraca0(ig)=0.
934	w0(ig)=0.
935	!!jyg le 27/04/2012
936	!! zlcl(ig)=0.
937	!!
938	endif
939	enddo
940
941	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
942	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
943
944	!------------Triggering------------------
945	IF (iflag_trig_bl.ge.1) THEN
946
947	!-----Initialisations
948	depth(:)=0.
949	n2(:)=0.
950	s2(:)=0.
951	s_max(:)=0.
952
953	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
954	do ig=1,ngrid
955	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
956	depth(ig)=zmax_moy(ig)-zlcl(ig)
957	hmin(ig)=hmincoef*zlcl(ig)
958	if (depth(ig).ge.10.) then
959	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
960	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
961	!!
962	!!jyg le 27/04/2012
963	!! s_max(ig)=s2(ig)*log(n2(ig))
964	!! if (n2(ig) .lt. 1) s_max(ig)=0.
965	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
966	!!fin jyg
967	else
968	s2(ig)=0.
969	n2(ig)=0.
970	s_max(ig)=0.
971	endif
972	enddo
973	! print *,'avant Calcul de Wmax ' !!jyg
974
975	!-----Calcul de Wmax et ALE_BL_STAT associée
976	!!jyg le 30/04/2012
977	!! do ig=1,ngrid
978	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
979	!! 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))))
980	!! ale_bl_stat(ig)=0.5w_max(ig)*2
981	!! else
982	!! w_max(ig)=0.
983	!! ale_bl_stat(ig)=0.
984	!! endif
985	!! enddo
986	susqr2pi=susqrt(2.Rpi)
987	Reuler=exp(1.)
988	do ig=1,ngrid
989	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
990	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
991	ale_bl_stat(ig)=0.5w_max(ig)*2
992	else
993	w_max(ig)=0.
994	ale_bl_stat(ig)=0.
995	endif
996	enddo
997
998	ENDIF ! iflag_trig_bl
999	! print *,'ENDIF iflag_trig_bl' !!jyg
1000	lfname='thermcell_main before closure'
1001	lvarname = 'pt'
1002	CALL check_var3D(lfname, lvarname, pt, klon, klev, largest, .FALSE.)
1003	lvarname = 'pdtadj'
1004	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, largest, .FALSE.)
1005	lvarname = 'alp_bl_det'
1006	CALL check_var(lfname, lvarname, alp_bl_det,klon, largest, .FALSE.)
1007	lvarname = 'rhobarz0'
1008	CALL check_var(lfname, lvarname, rhobarz0,klon, largest, .FALSE.)
1009	lvarname = 'fraca0'
1010	CALL check_var(lfname, lvarname, fraca0,klon, largest, .FALSE.)
1011	lvarname = 'w_conv'
1012	CALL check_var(lfname, lvarname, w_conv,klon, largest, .FALSE.)
1013	lvarname = 'w0'
1014	CALL check_var(lfname, lvarname, w0,klon, largest, .FALSE.)
1015	lvarname = 'alp_bl_fluct_m'
1016	CALL check_var(lfname, lvarname, alp_bl_fluct_m,klon, largest, .FALSE.)
1017	lvarname = 'alp_bl_fluct_tke'
1018	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,klon, largest, .FALSE.)
1019	lvarname = 'therm_tke_max0'
1020	CALL check_var(lfname, lvarname, therm_tke_max0,klon, largest, .FALSE.)
1021	lvarname = 'env_tke_max0'
1022	CALL check_var(lfname, lvarname, env_tke_max0,klon, largest, .FALSE.)
1023	lvarname = 'pbl_tke_max0'
1024	CALL check_var(lfname, lvarname, pbl_tke_max0,klon, largest, .FALSE.)
1025	lvarname = 'alp_bl_conv'
1026	CALL check_var(lfname, lvarname, alp_bl_conv,klon, largest, .FALSE.)
1027	lvarname = 'alp_bl_stat'
1028	CALL check_var(lfname, lvarname, alp_bl_stat,klon, largest, .FALSE.)
1029	!------------Closure------------------
1030
1031	IF (iflag_clos_bl.ge.1) THEN
1032
1033	!-----Calcul de ALP_BL_STAT
1034	do ig=1,ngrid
1035	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
1036	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
1037	& (w0(ig)**2)
1038	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
1039	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
1040	if (iflag_clos_bl.ge.2) then
1041	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
1042	& (w0(ig)**2)
1043	else
1044	alp_bl_conv(ig)=0.
1045	endif
1046	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
1047	enddo
1048
1049	!-----Sécurité ALP infinie
1050	do ig=1,ngrid
1051	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
1052	enddo
1053
1054	ENDIF ! (iflag_clos_bl.ge.1)
1055	lfname='after closure'
1056	lvarname = 'pt'
1057	CALL check_var3D(lfname, lvarname, pt, klon, klev, largest, .FALSE.)
1058	lvarname = 'pdtadj'
1059	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, largest, .FALSE.)
1060	lvarname = 'alp_bl_det'
1061	CALL check_var(lfname, lvarname, alp_bl_det,klon, largest, .FALSE.)
1062	lvarname = 'rhobarz0'
1063	CALL check_var(lfname, lvarname, rhobarz0,klon, largest, .FALSE.)
1064	lvarname = 'fraca0'
1065	CALL check_var(lfname, lvarname, fraca0,klon, largest, .FALSE.)
1066	lvarname = 'w_conv'
1067	CALL check_var(lfname, lvarname, w_conv,klon, largest, .FALSE.)
1068	lvarname = 'w0'
1069	CALL check_var(lfname, lvarname, w0,klon, largest, .FALSE.)
1070	lvarname = 'alp_bl_fluct_m'
1071	CALL check_var(lfname, lvarname, alp_bl_fluct_m,klon, largest, .FALSE.)
1072	lvarname = 'alp_bl_fluct_tke'
1073	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,klon, largest, .FALSE.)
1074	lvarname = 'therm_tke_max0'
1075	CALL check_var(lfname, lvarname, therm_tke_max0,klon, largest, .FALSE.)
1076	lvarname = 'env_tke_max0'
1077	CALL check_var(lfname, lvarname, env_tke_max0,klon, largest, .FALSE.)
1078	lvarname = 'pbl_tke_max0'
1079	CALL check_var(lfname, lvarname, pbl_tke_max0,klon, largest, .FALSE.)
1080	lvarname = 'alp_bl_conv'
1081	CALL check_var(lfname, lvarname, alp_bl_conv,klon, largest, .FALSE.)
1082	lvarname = 'alp_bl_stat'
1083	CALL check_var(lfname, lvarname, alp_bl_stat,klon, largest, .FALSE.)
1084
1085	!!! fin nrlmd le 10/04/2012
1086
1087	if (prt_level.ge.10) then
1088	ig=igout
1089	do l=1,nlay
1090	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
1091	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
1092	enddo
1093	endif
1094
1095	! print*,'avant calcul ale et alp'
1096	!calcul de ALE et ALP pour la convection
1097	Alp_bl(:)=0.
1098	Ale_bl(:)=0.
1099	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
1100	do l=1,nlay
1101	do ig=1,ngrid
1102	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
1103	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
1104	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
1105	enddo
1106	enddo
1107
1108	!test:calcul de la ponderation des couches pour KE
1109	!initialisations
1110
1111	fm_tot(:)=0.
1112	wght_th(:,:)=1.
1113	lalim_conv(:)=lalim(:)
1114
1115	do k=1,klev
1116	do ig=1,ngrid
1117	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
1118	enddo
1119	enddo
1120	lfname='after calculation of Al[p/e]_bl'
1121	lvarname = 'fm'
1122	CALL check_var3D(lfname, lvarname, fm, klon, klev, largest, .FALSE.)
1123	lvarname = 'rhobarz'
1124	CALL check_var3D(lfname, lvarname, rhobarz, klon, klev, largest, .FALSE.)
1125	lvarname = 'Alp_bl'
1126	CALL check_var(lfname, lvarname, Alp_bl,klon, largest, .FALSE.)
1127	lvarname = 'Ale_bl'
1128	CALL check_var(lfname, lvarname, Ale_bl,klon, largest, .FALSE.)
1129
1130	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
1131	! plus petit que 1.e-10, on prend wght_th=1.
1132	do k=1,klev
1133	do ig=1,ngrid
1134	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1135	wght_th(ig,k)=alim_star(ig,k)
1136	endif
1137	enddo
1138	enddo
1139
1140	! print*,'apres wght_th'
1141	!test pour prolonger la convection
1142	do ig=1,ngrid
1143	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1144	if ((alim_star(ig,1).lt.1.e-10)) then
1145	lalim_conv(ig)=1
1146	wght_th(ig,1)=1.
1147	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1148	endif
1149	enddo
1150
1151	!------------------------------------------------------------------------
1152	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1153	! Integrale verticale de ALP.
1154	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1155	! couches
1156	!------------------------------------------------------------------------
1157
1158	alp_int(:)=0.
1159	dp_int(:)=0.
1160	do l=2,nlay
1161	do ig=1,ngrid
1162	if(l.LE.lmax(ig)) THEN
1163	zdp=pplay(ig,l-1)-pplay(ig,l)
1164	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1165	dp_int(ig)=dp_int(ig)+zdp
1166	endif
1167	enddo
1168	enddo
1169
1170	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1171	do ig=1,ngrid
1172	!valeur integree de alp_bl * 0.5:
1173	if (dp_int(ig)>0.) then
1174	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1175	endif
1176	enddo!
1177	endif
1178
1179
1180	! Facteur multiplicatif sur Alp_bl
1181	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1182
1183	lfname='thermcell_main last computations on Alp_bl'
1184	lvarname = 'alp_int'
1185	CALL check_var(lfname, lvarname, alp_int, klon, largest*10.e4, .FALSE.)
1186	lvarname = 'dp_int'
1187	CALL check_var3D(lfname, lvarname, dp_int, klon, largest, .FALSE.)
1188	lvarname = 'Alp_bl'
1189	CALL check_var(lfname, lvarname, Alp_bl, klon, klev, largest, .FALSE.)
1190
1191	!------------------------------------------------------------------------
1192
1193
1194	!calcul du ratqscdiff
1195	if (prt_level.ge.1) print*,'14e OK convect8'
1196	var=0.
1197	vardiff=0.
1198	ratqsdiff(:,:)=0.
1199
1200	do l=1,klev
1201	do ig=1,ngrid
1202	if (l<=lalim(ig)) then
1203	var=var+alim_star(ig,l)zqta(ig,l)1000.
1204	endif
1205	enddo
1206	enddo
1207
1208	if (prt_level.ge.1) print*,'14f OK convect8'
1209
1210	do l=1,klev
1211	do ig=1,ngrid
1212	if (l<=lalim(ig)) then
1213	zf=fraca(ig,l)
1214	zf2=zf/(1.-zf)
1215	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1216	endif
1217	enddo
1218	enddo
1219
1220	if (prt_level.ge.1) print*,'14g OK convect8'
1221	do l=1,nlay
1222	do ig=1,ngrid
1223	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1224	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1225	enddo
1226	enddo
1227	!--------------------------------------------------------------------
1228	!
1229	!ecriture des fichiers sortie
1230	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1231
1232	#ifdef wrgrads_thermcell
1233	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1234	#include "thermcell_out3d.h"
1235	#endif
1236
1237	endif
1238
1239	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1240
1241	lfname='before leaving thermcell_main'
1242	lvarname = 'pt'
1243	CALL check_var3D(lfname, lvarname, pt, klon, klev, largest, .FALSE.)
1244	lvarname = 'pdtadj'
1245	CALL check_var3D(lfname, lvarname, pdtadj, klon, klev, largest, .FALSE.)
1246	lvarname = 'Alp_bl'
1247	CALL check_var(lfname, lvarname, Alp_bl, klon, klev, largest, .FALSE.)
1248	lvarname = 'Ale_bl'
1249	CALL check_var(lfname, lvarname, Ale_bl, klon, klev, largest, .FALSE.)
1250
1251
1252	return
1253	end
1254
1255	!-----------------------------------------------------------------------------
1256
1257	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1258	IMPLICIT NONE
1259	#include "iniprint.h"
1260
1261	integer i, k, klon,klev
1262	real pplev(klon,klev+1),pplay(klon,klev)
1263	real ztv(klon,klev)
1264	real po(klon,klev)
1265	real ztva(klon,klev)
1266	real zqla(klon,klev)
1267	real f_star(klon,klev)
1268	real zw2(klon,klev)
1269	integer long(klon)
1270	real seuil
1271	character*21 comment
1272
1273	if (prt_level.ge.1) THEN
1274	print*,'WARNING !!! TEST ',comment
1275	endif
1276	return
1277
1278	! test sur la hauteur des thermiques ...
1279	do i=1,klon
1280	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1281	if (prt_level.ge.10) then
1282	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1283	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1284	do k=1,klev
1285	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)
1286	enddo
1287	endif
1288	enddo
1289
1290
1291	return
1292	end
1293
1294	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1295	! On transporte pbl_tke pour donner therm_tke
1296	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1297	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1298	& rg,pplev,therm_tke_max)
1299	implicit none
1300
1301	#include "iniprint.h"
1302	!=======================================================================
1303	!
1304	! Calcul du transport verticale dans la couche limite en presence
1305	! de "thermiques" explicitement representes
1306	! calcul du dq/dt une fois qu'on connait les ascendances
1307	!
1308	!=======================================================================
1309
1310	integer ngrid,nlay,nsrf
1311
1312	real ptimestep
1313	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1314	real entr0(ngrid,nlay),rg
1315	real therm_tke_max(ngrid,nlay)
1316	real detr0(ngrid,nlay)
1317
1318
1319	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1320	real entr(ngrid,nlay)
1321	real q(ngrid,nlay)
1322	integer lev_out ! niveau pour les print
1323
1324	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1325
1326	real zzm
1327
1328	integer ig,k
1329	integer isrf
1330
1331
1332	lev_out=0
1333
1334
1335	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1336
1337	! calcul du detrainement
1338	do k=1,nlay
1339	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1340	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1341	enddo
1342
1343
1344	! Decalage vertical des entrainements et detrainements.
1345	masse(:,1)=0.5*masse0(:,1)
1346	entr(:,1)=0.5*entr0(:,1)
1347	detr(:,1)=0.5*detr0(:,1)
1348	fm(:,1)=0.
1349	do k=1,nlay-1
1350	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1351	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1352	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1353	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1354	enddo
1355	fm(:,nlay+1)=0.
1356
1357	!!! nrlmd le 16/09/2010
1358	! calcul de la valeur dans les ascendances
1359	! do ig=1,ngrid
1360	! qa(ig,1)=q(ig,1)
1361	! enddo
1362	!!!
1363
1364	!do isrf=1,nsrf
1365
1366	! q(:,:)=therm_tke(:,:,isrf)
1367	q(:,:)=therm_tke_max(:,:)
1368	!!! nrlmd le 16/09/2010
1369	do ig=1,ngrid
1370	qa(ig,1)=q(ig,1)
1371	enddo
1372	!!!
1373
1374	if (1==1) then
1375	do k=2,nlay
1376	do ig=1,ngrid
1377	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1378	& 1.e-5*masse(ig,k)) then
1379	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1380	& /(fm(ig,k+1)+detr(ig,k))
1381	else
1382	qa(ig,k)=q(ig,k)
1383	endif
1384	if (qa(ig,k).lt.0.) then
1385	! print*,'qa<0!!!'
1386	endif
1387	if (q(ig,k).lt.0.) then
1388	! print*,'q<0!!!'
1389	endif
1390	enddo
1391	enddo
1392
1393	! Calcul du flux subsident
1394
1395	do k=2,nlay
1396	do ig=1,ngrid
1397	wqd(ig,k)=fm(ig,k)*q(ig,k)
1398	if (wqd(ig,k).lt.0.) then
1399	! print*,'wqd<0!!!'
1400	endif
1401	enddo
1402	enddo
1403	do ig=1,ngrid
1404	wqd(ig,1)=0.
1405	wqd(ig,nlay+1)=0.
1406	enddo
1407
1408	! Calcul des tendances
1409	do k=1,nlay
1410	do ig=1,ngrid
1411	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1412	& -wqd(ig,k)+wqd(ig,k+1)) &
1413	& *ptimestep/masse(ig,k)
1414	enddo
1415	enddo
1416
1417	endif
1418
1419	therm_tke_max(:,:)=q(:,:)
1420
1421	return
1422	!!! fin nrlmd le 10/04/2012
1423	end
1424

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

Download in other formats:

Original Format