Context Navigation

thermcell_main.F90 @ 2306

Last change on this file since 2306 was 2193, checked in by fhourdin, 10 years ago

New option iflag_trig_bl for triggering :
1 : ALE computed from the PDF of thermals
2 : mean ALE with moist cape
3 : mean ALE with dry cape

Catherine et Frédéric

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 40.8 KB

Line
1	!
2	! $Id: thermcell_main.F90 2193 2015-02-04 06:48:02Z jyg $
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	LOGICAL debut
86
87	! local:
88	! ------
89
90	integer icount
91
92	integer, save :: dvdq=1,dqimpl=-1
93	!$OMP THREADPRIVATE(dvdq,dqimpl)
94	data icount/0/
95	save icount
96	!$OMP THREADPRIVATE(icount)
97
98	integer,save :: igout=1
99	!$OMP THREADPRIVATE(igout)
100	integer,save :: lunout1=6
101	!$OMP THREADPRIVATE(lunout1)
102	integer,save :: lev_out=10
103	!$OMP THREADPRIVATE(lev_out)
104
105	REAL susqr2pi, Reuler
106
107	INTEGER ig,k,l,ll,ierr
108	real zsortie1d(klon)
109	INTEGER lmax(klon),lmin(klon),lalim(klon)
110	INTEGER lmix(klon)
111	INTEGER lmix_bis(klon)
112	real linter(klon)
113	real zmix(klon)
114	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
115	! real fraca(klon,klev)
116
117	real zmax_sec(klon)
118	!on garde le zmax du pas de temps precedent
119	real zmax0(klon)
120	!FH/IM save zmax0
121
122	real lambda
123
124	real zlev(klon,klev+1),zlay(klon,klev)
125	real deltaz(klon,klev)
126	REAL zh(klon,klev)
127	real zthl(klon,klev),zdthladj(klon,klev)
128	REAL ztv(klon,klev)
129	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
130	real zl(klon,klev)
131	real zsortie(klon,klev)
132	real zva(klon,klev)
133	real zua(klon,klev)
134	real zoa(klon,klev)
135
136	real zta(klon,klev)
137	real zha(klon,klev)
138	real fraca(klon,klev+1)
139	real zf,zf2
140	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
141	real q2(klon,klev)
142	! FH probleme de dimensionnement avec l'allocation dynamique
143	! common/comtherm/thetath2,wth2
144	real wq(klon,klev)
145	real wthl(klon,klev)
146	real wthv(klon,klev)
147
148	real ratqscth(klon,klev)
149	real var
150	real vardiff
151	real ratqsdiff(klon,klev)
152
153	logical sorties
154	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
155	real zpspsk(klon,klev)
156
157	real wmax(klon)
158	real wmax_tmp(klon)
159	real wmax_sec(klon)
160	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
161	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
162
163	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
164	!niveau de condensation
165	integer nivcon(klon)
166	real zcon(klon)
167	REAL CHI
168	real zcon2(klon)
169	real pcon(klon)
170	real zqsat(klon,klev)
171	real zqsatth(klon,klev)
172
173	real f_star(klon,klev+1),entr_star(klon,klev)
174	real detr_star(klon,klev)
175	real alim_star_tot(klon)
176	real alim_star(klon,klev)
177	real alim_star_clos(klon,klev)
178	real f(klon), f0(klon)
179	!FH/IM save f0
180	real zlevinter(klon)
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	!-----------------------------------------------------------------------
254	! initialisation:
255	! ---------------
256	!
257
258	seuil=0.25
259
260	if (debut) then
261	if (iflag_thermals==15.or.iflag_thermals==16) then
262	dvdq=0
263	dqimpl=-1
264	else
265	dvdq=1
266	dqimpl=1
267	endif
268
269	fm0=0.
270	entr0=0.
271	detr0=0.
272	endif
273	fm=0. ; entr=0. ; detr=0.
274	icount=icount+1
275
276	!IM 090508 beg
277	!print*,'====================================================================='
278	!print*,'====================================================================='
279	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
280	!print*,'====================================================================='
281	!print*,'====================================================================='
282	!IM 090508 end
283
284	if (prt_level.ge.1) print*,'thermcell_main V4'
285
286	sorties=.true.
287	IF(ngrid.NE.klon) THEN
288	PRINT*
289	PRINT*,'STOP dans convadj'
290	PRINT*,'ngrid =',ngrid
291	PRINT*,'klon =',klon
292	ENDIF
293	!
294	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
295	do ig=1,klon
296	f0(ig)=max(f0(ig),1.e-2)
297	zmax0(ig)=max(zmax0(ig),40.)
298	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
299	enddo
300
301	if (prt_level.ge.20) then
302	do ig=1,ngrid
303	print*,'th_main ig f0',ig,f0(ig)
304	enddo
305	endif
306	!-----------------------------------------------------------------------
307	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
308	! --------------------------------------------------------------------
309	!
310	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
311	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
312
313	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
314
315	!------------------------------------------------------------------------
316	! --------------------
317	!
318	!
319	! + + + + + + + + + + +
320	!
321	!
322	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
323	! wh,wt,wo ...
324	!
325	! + + + + + + + + + + + zh,zu,zv,zo,rho
326	!
327	!
328	! -------------------- zlev(1)
329	! \\\\\\\\\\\\\\\\\\\\
330	!
331	!
332
333	!-----------------------------------------------------------------------
334	! Calcul des altitudes des couches
335	!-----------------------------------------------------------------------
336
337	do l=2,nlay
338	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
339	enddo
340	zlev(:,1)=0.
341	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
342	do l=1,nlay
343	zlay(:,l)=pphi(:,l)/RG
344	enddo
345	!calcul de l epaisseur des couches
346	do l=1,nlay
347	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
348	enddo
349
350	! print*,'2 OK convect8'
351	!-----------------------------------------------------------------------
352	! Calcul des densites
353	!-----------------------------------------------------------------------
354
355	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
356
357	if (prt_level.ge.10)write(lunout,*) &
358	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
359	rhobarz(:,1)=rho(:,1)
360
361	do l=2,nlay
362	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
363	enddo
364
365	!calcul de la masse
366	do l=1,nlay
367	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
368	enddo
369
370	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
371
372	!------------------------------------------------------------------
373	!
374	! /\|\
375	! -------- \| F_k+1 -------
376	! ----> D_k
377	! /\|\ <---- E_k , A_k
378	! -------- \| F_k ---------
379	! ----> D_k-1
380	! <---- E_k-1 , A_k-1
381	!
382	!
383	!
384	!
385	!
386	! ---------------------------
387	!
388	! ----- F_lmax+1=0 ---------- \
389	! lmax (zmax) \|
390	! --------------------------- \|
391	! \|
392	! --------------------------- \|
393	! \|
394	! --------------------------- \|
395	! \|
396	! --------------------------- \|
397	! \|
398	! --------------------------- \|
399	! \| E
400	! --------------------------- \| D
401	! \|
402	! --------------------------- \|
403	! \|
404	! --------------------------- \ \|
405	! lalim \| \|
406	! --------------------------- \| \|
407	! \| \|
408	! --------------------------- \| \|
409	! \| A \|
410	! --------------------------- \| \|
411	! \| \|
412	! --------------------------- \| \|
413	! lmin (=1 pour le moment) \| \|
414	! ----- F_lmin=0 ------------ / /
415	!
416	! ---------------------------
417	! //////////////////////////
418	!
419	!
420	!=============================================================================
421	! Calculs initiaux ne faisant pas intervenir les changements de phase
422	!=============================================================================
423
424	!------------------------------------------------------------------
425	! 1. alim_star est le profil vertical de l'alimentation a la base du
426	! panache thermique, calcule a partir de la flotabilite de l'air sec
427	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
428	!------------------------------------------------------------------
429	!
430	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
431	lmin=1
432
433	!-----------------------------------------------------------------------------
434	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
435	! panache sec conservatif (e=d=0) alimente selon alim_star
436	! Il s'agit d'un calcul de type CAPE
437	! zmax_sec est utilise pour determiner la geometrie du thermique.
438	!------------------------------------------------------------------------------
439	!---------------------------------------------------------------------------------
440	!calcul du melange et des variables dans le thermique
441	!--------------------------------------------------------------------------------
442	!
443	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
444	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
445
446	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
447	! de la variable iflag_thermals
448
449	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
450	if (iflag_thermals_ed<=9) then
451	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
452	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
453	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
454	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
455	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
456	& ,lev_out,lunout1,igout)
457
458	elseif (iflag_thermals_ed>9) then
459	! print*,'THERM RIO et al 2010, version d Arnaud'
460	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
461	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
462	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
463	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
464	& ,lev_out,lunout1,igout)
465
466	endif
467
468	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
469
470	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
471	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
472
473	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
474	if (prt_level.ge.10) then
475	write(lunout1,*) 'Dans thermcell_main 2'
476	write(lunout1,*) 'lmin ',lmin(igout)
477	write(lunout1,*) 'lalim ',lalim(igout)
478	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
479	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
480	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
481	endif
482
483	!-------------------------------------------------------------------------------
484	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
485	!-------------------------------------------------------------------------------
486	!
487	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
488	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
489	! Attention, w2 est transforme en sa racine carree dans cette routine
490	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
491	wmax_tmp=0.
492	do l=1,nlay
493	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
494	enddo
495	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
496
497
498
499	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
500	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
501	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
502	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
503
504	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
505
506	!-------------------------------------------------------------------------------
507	! Fermeture,determination de f
508	!-------------------------------------------------------------------------------
509	!
510	!
511	!! write(lunout,*)'THERM NOUVEAU XXXXX'
512	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
513	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
514
515
516	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
517	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
518
519	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
520	if (prt_level.ge.10) then
521	write(lunout1,*) 'Dans thermcell_main 1b'
522	write(lunout1,*) 'lmin ',lmin(igout)
523	write(lunout1,*) 'lalim ',lalim(igout)
524	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
525	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
526	& ,l=1,lalim(igout)+4)
527	endif
528
529
530
531
532	! Choix de la fonction d'alimentation utilisee pour la fermeture.
533	! Apparemment sans importance
534	alim_star_clos(:,:)=alim_star(:,:)
535	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
536	!
537	!CR Appel de la fermeture seche
538	if (iflag_thermals_closure.eq.1) then
539
540	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
541	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
542
543	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
544	! Appel avec les zmax et wmax tenant compte de la condensation
545	! Semble moins bien marcher
546	else if (iflag_thermals_closure.eq.2) then
547
548	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
549	& zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
550
551	endif
552
553	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
554
555	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
556
557	if (tau_thermals>1.) then
558	lambda=exp(-ptimestep/tau_thermals)
559	f0=(1.-lambda)f+lambdaf0
560	else
561	f0=f
562	endif
563
564	! Test valable seulement en 1D mais pas genant
565	if (.not. (f0(1).ge.0.) ) then
566	abort_message = '.not. (f0(1).ge.0.)'
567	CALL abort_gcm (modname,abort_message,1)
568	endif
569
570	!-------------------------------------------------------------------------------
571	!deduction des flux
572	!-------------------------------------------------------------------------------
573
574	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
575	& lalim,lmax,alim_star, &
576	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
577	& detr,zqla,lev_out,lunout1,igout)
578	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
579
580	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
581	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
582	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
583
584	!------------------------------------------------------------------
585	! On ne prend pas directement les profils issus des calculs precedents
586	! mais on s'autorise genereusement une relaxation vers ceci avec
587	! une constante de temps tau_thermals (typiquement 1800s).
588	!------------------------------------------------------------------
589
590	if (tau_thermals>1.) then
591	lambda=exp(-ptimestep/tau_thermals)
592	fm0=(1.-lambda)fm+lambdafm0
593	entr0=(1.-lambda)entr+lambdaentr0
594	detr0=(1.-lambda)detr+lambdadetr0
595	else
596	fm0=fm
597	entr0=entr
598	detr0=detr
599	endif
600
601	!c------------------------------------------------------------------
602	! calcul du transport vertical
603	!------------------------------------------------------------------
604
605	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
606	& zthl,zdthladj,zta,lev_out)
607	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
608	& po,pdoadj,zoa,lev_out)
609
610	!------------------------------------------------------------------
611	! Calcul de la fraction de l'ascendance
612	!------------------------------------------------------------------
613	do ig=1,klon
614	fraca(ig,1)=0.
615	fraca(ig,nlay+1)=0.
616	enddo
617	do l=2,nlay
618	do ig=1,klon
619	if (zw2(ig,l).gt.1.e-10) then
620	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
621	else
622	fraca(ig,l)=0.
623	endif
624	enddo
625	enddo
626
627	!------------------------------------------------------------------
628	! calcul du transport vertical du moment horizontal
629	!------------------------------------------------------------------
630
631	!IM 090508
632	if (dvdq == 0 ) then
633
634	! Calcul du transport de V tenant compte d'echange par gradient
635	! de pression horizontal avec l'environnement
636
637	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
638	! & ,fraca*dvdq,zmax &
639	& ,fraca,zmax &
640	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
641
642	else
643
644	! calcul purement conservatif pour le transport de V
645	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
646	& ,zu,pduadj,zua,lev_out)
647	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
648	& ,zv,pdvadj,zva,lev_out)
649
650	endif
651
652	! print*,'13 OK convect8'
653	do l=1,nlay
654	do ig=1,ngrid
655	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
656	enddo
657	enddo
658
659	if (prt_level.ge.1) print*,'14 OK convect8'
660	!------------------------------------------------------------------
661	! Calculs de diagnostiques pour les sorties
662	!------------------------------------------------------------------
663	!calcul de fraca pour les sorties
664
665	if (sorties) then
666	if (prt_level.ge.1) print*,'14a OK convect8'
667	! calcul du niveau de condensation
668	! initialisation
669	do ig=1,ngrid
670	nivcon(ig)=0
671	zcon(ig)=0.
672	enddo
673	!nouveau calcul
674	do ig=1,ngrid
675	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
676	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
677	enddo
678	!IM do k=1,nlay
679	do k=1,nlay-1
680	do ig=1,ngrid
681	if ((pcon(ig).le.pplay(ig,k)) &
682	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
683	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
684	endif
685	enddo
686	enddo
687	!IM
688	ierr=0
689	do ig=1,ngrid
690	if (pcon(ig).le.pplay(ig,nlay)) then
691	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
692	ierr=1
693	endif
694	enddo
695	if (ierr==1) then
696	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
697	CALL abort_gcm (modname,abort_message,1)
698	endif
699
700	if (prt_level.ge.1) print*,'14b OK convect8'
701	do k=nlay,1,-1
702	do ig=1,ngrid
703	if (zqla(ig,k).gt.1e-10) then
704	nivcon(ig)=k
705	zcon(ig)=zlev(ig,k)
706	endif
707	enddo
708	enddo
709	if (prt_level.ge.1) print*,'14c OK convect8'
710	!calcul des moments
711	!initialisation
712	do l=1,nlay
713	do ig=1,ngrid
714	q2(ig,l)=0.
715	wth2(ig,l)=0.
716	wth3(ig,l)=0.
717	ratqscth(ig,l)=0.
718	ratqsdiff(ig,l)=0.
719	enddo
720	enddo
721	if (prt_level.ge.1) print*,'14d OK convect8'
722	if (prt_level.ge.10)write(lunout,*) &
723	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
724	do l=1,nlay
725	do ig=1,ngrid
726	zf=fraca(ig,l)
727	zf2=zf/(1.-zf)
728	!
729	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
730	if(zw2(ig,l).gt.1.e-10) then
731	wth2(ig,l)=zf2(zw2(ig,l))*2
732	else
733	wth2(ig,l)=0.
734	endif
735	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
736	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
737	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
738	!test: on calcul q2/po=ratqsc
739	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
740	enddo
741	enddo
742	!calcul des flux: q, thetal et thetav
743	do l=1,nlay
744	do ig=1,ngrid
745	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
746	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
747	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
748	enddo
749	enddo
750	!
751
752	!!! nrlmd le 10/04/2012
753
754	!------------Test sur le LCL des thermiques
755	do ig=1,ngrid
756	ok_lcl(ig)=.false.
757	if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
758	enddo
759
760	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
761	do l=1,nlay-1
762	do ig=1,ngrid
763	if (ok_lcl(ig)) then
764	!ATTENTION,zw2 calcule en pplev
765	! if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
766	! klcl(ig)=l
767	! interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
768	! endif
769	if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
770	klcl(ig)=l
771	interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
772	endif
773	endif
774	enddo
775	enddo
776
777	!------------Hauteur des thermiques
778	!!jyg le 27/04/2012
779	!! do ig =1,ngrid
780	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
781	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
782	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
783	!! zmax(ig)=pphi(ig,lmax(ig))/rg
784	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
785	!! enddo
786	do ig =1,ngrid
787	!CR:REHABILITATION ZMAX CONTINU
788	! zmax(ig)=pphi(ig,lmax(ig))/rg
789	if (ok_lcl(ig)) then
790	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
791	& -rhobarz(ig,klcl(ig)))*interp(ig)
792	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
793	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
794	else
795	rhobarz0(ig)=0.
796	zlcl(ig)=zmax(ig)
797	endif
798	enddo
799	!!jyg fin
800
801	!------------Calcul des propriétés du thermique au LCL
802	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
803
804	!-----Initialisation de la TKE moyenne
805	do l=1,nlay
806	do ig=1,ngrid
807	pbl_tke_max(ig,l)=0.
808	enddo
809	enddo
810
811	!-----Calcul de la TKE moyenne
812	do nsrf=1,nbsrf
813	do l=1,nlay
814	do ig=1,ngrid
815	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
816	enddo
817	enddo
818	enddo
819
820	!-----Initialisations des TKE dans et hors des thermiques
821	do l=1,nlay
822	do ig=1,ngrid
823	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
824	env_tke_max(ig,l)=pbl_tke_max(ig,l)
825	enddo
826	enddo
827
828	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
829	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
830	& rg,pplev,therm_tke_max)
831	! print *,' thermcell_tke_transport -> ' !!jyg
832
833	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
834	do l=1,nlay
835	do ig=1,ngrid
836	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
837	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
838	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
839	enddo
840	enddo
841	! print *,' apres w_ls = ' !!jyg
842
843	do ig=1,ngrid
844	if (ok_lcl(ig)) then
845	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
846	& -fraca(ig,klcl(ig)))*interp(ig)
847	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
848	& -zw2(ig,klcl(ig)))*interp(ig)
849	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
850	& -w_ls(ig,klcl(ig)))*interp(ig)
851	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
852	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
853	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
854	& -env_tke_max(ig,klcl(ig)))*interp(ig)
855	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
856	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
857	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
858	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
859	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
860	else
861	fraca0(ig)=0.
862	w0(ig)=0.
863	!!jyg le 27/04/2012
864	!! zlcl(ig)=0.
865	!!
866	endif
867	enddo
868
869	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
870	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
871
872	!------------Triggering------------------
873	IF (iflag_trig_bl.ge.1) THEN
874
875	!-----Initialisations
876	depth(:)=0.
877	n2(:)=0.
878	s2(:)=100. ! some low value, arbitrary
879	s_max(:)=0.
880
881	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
882	do ig=1,ngrid
883	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
884	depth(ig)=zmax_moy(ig)-zlcl(ig)
885	hmin(ig)=hmincoef*zlcl(ig)
886	if (depth(ig).ge.10.) then
887	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
888	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
889	!!
890	!!jyg le 27/04/2012
891	!! s_max(ig)=s2(ig)*log(n2(ig))
892	!! if (n2(ig) .lt. 1) s_max(ig)=0.
893	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
894	!!fin jyg
895	else
896	n2(ig)=0.
897	s_max(ig)=0.
898	endif
899	enddo
900	! print *,'avant Calcul de Wmax ' !!jyg
901
902	!-----Calcul de Wmax et ALE_BL_STAT associée
903	!!jyg le 30/04/2012
904	!! do ig=1,ngrid
905	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
906	!! 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))))
907	!! ale_bl_stat(ig)=0.5w_max(ig)*2
908	!! else
909	!! w_max(ig)=0.
910	!! ale_bl_stat(ig)=0.
911	!! endif
912	!! enddo
913	susqr2pi=susqrt(2.Rpi)
914	Reuler=exp(1.)
915	do ig=1,ngrid
916	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
917	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
918	ale_bl_stat(ig)=0.5w_max(ig)*2
919	else
920	w_max(ig)=0.
921	ale_bl_stat(ig)=0.
922	endif
923	enddo
924
925	ENDIF ! iflag_trig_bl
926	! print *,'ENDIF iflag_trig_bl' !!jyg
927
928	!------------Closure------------------
929
930	IF (iflag_clos_bl.ge.2) THEN
931
932	!-----Calcul de ALP_BL_STAT
933	do ig=1,ngrid
934	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
935	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
936	& (w0(ig)**2)
937	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
938	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
939	if (iflag_clos_bl.ge.2) then
940	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
941	& (w0(ig)**2)
942	else
943	alp_bl_conv(ig)=0.
944	endif
945	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
946	enddo
947
948	!-----Sécurité ALP infinie
949	do ig=1,ngrid
950	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
951	enddo
952
953	ENDIF ! (iflag_clos_bl.ge.2)
954
955	!!! fin nrlmd le 10/04/2012
956
957	if (prt_level.ge.10) then
958	ig=igout
959	do l=1,nlay
960	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
961	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
962	enddo
963	endif
964
965	! print*,'avant calcul ale et alp'
966	!calcul de ALE et ALP pour la convection
967	Alp_bl(:)=0.
968	Ale_bl(:)=0.
969	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
970	do l=1,nlay
971	do ig=1,ngrid
972	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
973	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
974	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
975	enddo
976	enddo
977
978	! Ale sec (max de wmax/2 sous la zone d'inhibition) dans
979	! le cas iflag_trig_bl=3
980	IF (iflag_trig_bl==3) Ale_bl(:)=0.5wmax_sec(:)*2
981
982	!test:calcul de la ponderation des couches pour KE
983	!initialisations
984
985	fm_tot(:)=0.
986	wght_th(:,:)=1.
987	lalim_conv(:)=lalim(:)
988
989	do k=1,klev
990	do ig=1,ngrid
991	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
992	enddo
993	enddo
994
995	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
996	! plus petit que 1.e-10, on prend wght_th=1.
997	do k=1,klev
998	do ig=1,ngrid
999	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1000	wght_th(ig,k)=alim_star(ig,k)
1001	endif
1002	enddo
1003	enddo
1004
1005	! print*,'apres wght_th'
1006	!test pour prolonger la convection
1007	do ig=1,ngrid
1008	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1009	if ((alim_star(ig,1).lt.1.e-10)) then
1010	lalim_conv(ig)=1
1011	wght_th(ig,1)=1.
1012	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1013	endif
1014	enddo
1015
1016	!------------------------------------------------------------------------
1017	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1018	! Integrale verticale de ALP.
1019	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1020	! couches
1021	!------------------------------------------------------------------------
1022
1023	alp_int(:)=0.
1024	dp_int(:)=0.
1025	do l=2,nlay
1026	do ig=1,ngrid
1027	if(l.LE.lmax(ig)) THEN
1028	zdp=pplay(ig,l-1)-pplay(ig,l)
1029	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1030	dp_int(ig)=dp_int(ig)+zdp
1031	endif
1032	enddo
1033	enddo
1034
1035	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1036	do ig=1,ngrid
1037	!valeur integree de alp_bl * 0.5:
1038	if (dp_int(ig)>0.) then
1039	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1040	endif
1041	enddo!
1042	endif
1043
1044
1045	! Facteur multiplicatif sur Alp_bl
1046	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1047
1048	!------------------------------------------------------------------------
1049
1050
1051	!calcul du ratqscdiff
1052	if (prt_level.ge.1) print*,'14e OK convect8'
1053	var=0.
1054	vardiff=0.
1055	ratqsdiff(:,:)=0.
1056
1057	do l=1,klev
1058	do ig=1,ngrid
1059	if (l<=lalim(ig)) then
1060	var=var+alim_star(ig,l)zqta(ig,l)1000.
1061	endif
1062	enddo
1063	enddo
1064
1065	if (prt_level.ge.1) print*,'14f OK convect8'
1066
1067	do l=1,klev
1068	do ig=1,ngrid
1069	if (l<=lalim(ig)) then
1070	zf=fraca(ig,l)
1071	zf2=zf/(1.-zf)
1072	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1073	endif
1074	enddo
1075	enddo
1076
1077	if (prt_level.ge.1) print*,'14g OK convect8'
1078	do l=1,nlay
1079	do ig=1,ngrid
1080	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1081	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1082	enddo
1083	enddo
1084	!--------------------------------------------------------------------
1085	!
1086	!ecriture des fichiers sortie
1087	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1088
1089	endif
1090
1091	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1092
1093	return
1094	end
1095
1096	!-----------------------------------------------------------------------------
1097
1098	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1099	IMPLICIT NONE
1100	#include "iniprint.h"
1101
1102	integer i, k, klon,klev
1103	real pplev(klon,klev+1),pplay(klon,klev)
1104	real ztv(klon,klev)
1105	real po(klon,klev)
1106	real ztva(klon,klev)
1107	real zqla(klon,klev)
1108	real f_star(klon,klev)
1109	real zw2(klon,klev)
1110	integer long(klon)
1111	real seuil
1112	character*21 comment
1113
1114	if (prt_level.ge.1) THEN
1115	print*,'WARNING !!! TEST ',comment
1116	endif
1117	return
1118
1119	! test sur la hauteur des thermiques ...
1120	do i=1,klon
1121	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1122	if (prt_level.ge.10) then
1123	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1124	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1125	do k=1,klev
1126	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)
1127	enddo
1128	endif
1129	enddo
1130
1131
1132	return
1133	end
1134
1135	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1136	! On transporte pbl_tke pour donner therm_tke
1137	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1138	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1139	& rg,pplev,therm_tke_max)
1140	implicit none
1141
1142	#include "iniprint.h"
1143	!=======================================================================
1144	!
1145	! Calcul du transport verticale dans la couche limite en presence
1146	! de "thermiques" explicitement representes
1147	! calcul du dq/dt une fois qu'on connait les ascendances
1148	!
1149	!=======================================================================
1150
1151	integer ngrid,nlay,nsrf
1152
1153	real ptimestep
1154	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1155	real entr0(ngrid,nlay),rg
1156	real therm_tke_max(ngrid,nlay)
1157	real detr0(ngrid,nlay)
1158
1159
1160	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1161	real entr(ngrid,nlay)
1162	real q(ngrid,nlay)
1163	integer lev_out ! niveau pour les print
1164
1165	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1166
1167	real zzm
1168
1169	integer ig,k
1170	integer isrf
1171
1172
1173	lev_out=0
1174
1175
1176	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1177
1178	! calcul du detrainement
1179	do k=1,nlay
1180	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1181	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1182	enddo
1183
1184
1185	! Decalage vertical des entrainements et detrainements.
1186	masse(:,1)=0.5*masse0(:,1)
1187	entr(:,1)=0.5*entr0(:,1)
1188	detr(:,1)=0.5*detr0(:,1)
1189	fm(:,1)=0.
1190	do k=1,nlay-1
1191	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1192	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1193	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1194	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1195	enddo
1196	fm(:,nlay+1)=0.
1197
1198	!!! nrlmd le 16/09/2010
1199	! calcul de la valeur dans les ascendances
1200	! do ig=1,ngrid
1201	! qa(ig,1)=q(ig,1)
1202	! enddo
1203	!!!
1204
1205	!do isrf=1,nsrf
1206
1207	! q(:,:)=therm_tke(:,:,isrf)
1208	q(:,:)=therm_tke_max(:,:)
1209	!!! nrlmd le 16/09/2010
1210	do ig=1,ngrid
1211	qa(ig,1)=q(ig,1)
1212	enddo
1213	!!!
1214
1215	if (1==1) then
1216	do k=2,nlay
1217	do ig=1,ngrid
1218	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1219	& 1.e-5*masse(ig,k)) then
1220	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1221	& /(fm(ig,k+1)+detr(ig,k))
1222	else
1223	qa(ig,k)=q(ig,k)
1224	endif
1225	if (qa(ig,k).lt.0.) then
1226	! print*,'qa<0!!!'
1227	endif
1228	if (q(ig,k).lt.0.) then
1229	! print*,'q<0!!!'
1230	endif
1231	enddo
1232	enddo
1233
1234	! Calcul du flux subsident
1235
1236	do k=2,nlay
1237	do ig=1,ngrid
1238	wqd(ig,k)=fm(ig,k)*q(ig,k)
1239	if (wqd(ig,k).lt.0.) then
1240	! print*,'wqd<0!!!'
1241	endif
1242	enddo
1243	enddo
1244	do ig=1,ngrid
1245	wqd(ig,1)=0.
1246	wqd(ig,nlay+1)=0.
1247	enddo
1248
1249	! Calcul des tendances
1250	do k=1,nlay
1251	do ig=1,ngrid
1252	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1253	& -wqd(ig,k)+wqd(ig,k+1)) &
1254	& *ptimestep/masse(ig,k)
1255	enddo
1256	enddo
1257
1258	endif
1259
1260	therm_tke_max(:,:)=q(:,:)
1261
1262	return
1263	!!! fin nrlmd le 10/04/2012
1264	end
1265

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

Download in other formats:

Original Format