Context Navigation

thermcell_main.F90 @ 2153

Last change on this file since 2153 was 2079, checked in by lguez, 10 years ago
Protect against division by 0 in computation of proba_notrig in procedure physiq: initialize s2 at some low value instead of 0.
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.6 KB

Line
1	!
2	! $Id: thermcell_main.F90 2079 2014-07-07 13:38:36Z lguez $
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	!test:calcul de la ponderation des couches pour KE
979	!initialisations
980
981	fm_tot(:)=0.
982	wght_th(:,:)=1.
983	lalim_conv(:)=lalim(:)
984
985	do k=1,klev
986	do ig=1,ngrid
987	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
988	enddo
989	enddo
990
991	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
992	! plus petit que 1.e-10, on prend wght_th=1.
993	do k=1,klev
994	do ig=1,ngrid
995	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
996	wght_th(ig,k)=alim_star(ig,k)
997	endif
998	enddo
999	enddo
1000
1001	! print*,'apres wght_th'
1002	!test pour prolonger la convection
1003	do ig=1,ngrid
1004	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1005	if ((alim_star(ig,1).lt.1.e-10)) then
1006	lalim_conv(ig)=1
1007	wght_th(ig,1)=1.
1008	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1009	endif
1010	enddo
1011
1012	!------------------------------------------------------------------------
1013	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1014	! Integrale verticale de ALP.
1015	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1016	! couches
1017	!------------------------------------------------------------------------
1018
1019	alp_int(:)=0.
1020	dp_int(:)=0.
1021	do l=2,nlay
1022	do ig=1,ngrid
1023	if(l.LE.lmax(ig)) THEN
1024	zdp=pplay(ig,l-1)-pplay(ig,l)
1025	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1026	dp_int(ig)=dp_int(ig)+zdp
1027	endif
1028	enddo
1029	enddo
1030
1031	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1032	do ig=1,ngrid
1033	!valeur integree de alp_bl * 0.5:
1034	if (dp_int(ig)>0.) then
1035	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1036	endif
1037	enddo!
1038	endif
1039
1040
1041	! Facteur multiplicatif sur Alp_bl
1042	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1043
1044	!------------------------------------------------------------------------
1045
1046
1047	!calcul du ratqscdiff
1048	if (prt_level.ge.1) print*,'14e OK convect8'
1049	var=0.
1050	vardiff=0.
1051	ratqsdiff(:,:)=0.
1052
1053	do l=1,klev
1054	do ig=1,ngrid
1055	if (l<=lalim(ig)) then
1056	var=var+alim_star(ig,l)zqta(ig,l)1000.
1057	endif
1058	enddo
1059	enddo
1060
1061	if (prt_level.ge.1) print*,'14f OK convect8'
1062
1063	do l=1,klev
1064	do ig=1,ngrid
1065	if (l<=lalim(ig)) then
1066	zf=fraca(ig,l)
1067	zf2=zf/(1.-zf)
1068	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1069	endif
1070	enddo
1071	enddo
1072
1073	if (prt_level.ge.1) print*,'14g OK convect8'
1074	do l=1,nlay
1075	do ig=1,ngrid
1076	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1077	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1078	enddo
1079	enddo
1080	!--------------------------------------------------------------------
1081	!
1082	!ecriture des fichiers sortie
1083	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1084
1085	endif
1086
1087	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1088
1089	return
1090	end
1091
1092	!-----------------------------------------------------------------------------
1093
1094	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1095	IMPLICIT NONE
1096	#include "iniprint.h"
1097
1098	integer i, k, klon,klev
1099	real pplev(klon,klev+1),pplay(klon,klev)
1100	real ztv(klon,klev)
1101	real po(klon,klev)
1102	real ztva(klon,klev)
1103	real zqla(klon,klev)
1104	real f_star(klon,klev)
1105	real zw2(klon,klev)
1106	integer long(klon)
1107	real seuil
1108	character*21 comment
1109
1110	if (prt_level.ge.1) THEN
1111	print*,'WARNING !!! TEST ',comment
1112	endif
1113	return
1114
1115	! test sur la hauteur des thermiques ...
1116	do i=1,klon
1117	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1118	if (prt_level.ge.10) then
1119	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1120	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1121	do k=1,klev
1122	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)
1123	enddo
1124	endif
1125	enddo
1126
1127
1128	return
1129	end
1130
1131	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1132	! On transporte pbl_tke pour donner therm_tke
1133	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1134	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1135	& rg,pplev,therm_tke_max)
1136	implicit none
1137
1138	#include "iniprint.h"
1139	!=======================================================================
1140	!
1141	! Calcul du transport verticale dans la couche limite en presence
1142	! de "thermiques" explicitement representes
1143	! calcul du dq/dt une fois qu'on connait les ascendances
1144	!
1145	!=======================================================================
1146
1147	integer ngrid,nlay,nsrf
1148
1149	real ptimestep
1150	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1151	real entr0(ngrid,nlay),rg
1152	real therm_tke_max(ngrid,nlay)
1153	real detr0(ngrid,nlay)
1154
1155
1156	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1157	real entr(ngrid,nlay)
1158	real q(ngrid,nlay)
1159	integer lev_out ! niveau pour les print
1160
1161	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1162
1163	real zzm
1164
1165	integer ig,k
1166	integer isrf
1167
1168
1169	lev_out=0
1170
1171
1172	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1173
1174	! calcul du detrainement
1175	do k=1,nlay
1176	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1177	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1178	enddo
1179
1180
1181	! Decalage vertical des entrainements et detrainements.
1182	masse(:,1)=0.5*masse0(:,1)
1183	entr(:,1)=0.5*entr0(:,1)
1184	detr(:,1)=0.5*detr0(:,1)
1185	fm(:,1)=0.
1186	do k=1,nlay-1
1187	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1188	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1189	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1190	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1191	enddo
1192	fm(:,nlay+1)=0.
1193
1194	!!! nrlmd le 16/09/2010
1195	! calcul de la valeur dans les ascendances
1196	! do ig=1,ngrid
1197	! qa(ig,1)=q(ig,1)
1198	! enddo
1199	!!!
1200
1201	!do isrf=1,nsrf
1202
1203	! q(:,:)=therm_tke(:,:,isrf)
1204	q(:,:)=therm_tke_max(:,:)
1205	!!! nrlmd le 16/09/2010
1206	do ig=1,ngrid
1207	qa(ig,1)=q(ig,1)
1208	enddo
1209	!!!
1210
1211	if (1==1) then
1212	do k=2,nlay
1213	do ig=1,ngrid
1214	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1215	& 1.e-5*masse(ig,k)) then
1216	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1217	& /(fm(ig,k+1)+detr(ig,k))
1218	else
1219	qa(ig,k)=q(ig,k)
1220	endif
1221	if (qa(ig,k).lt.0.) then
1222	! print*,'qa<0!!!'
1223	endif
1224	if (q(ig,k).lt.0.) then
1225	! print*,'q<0!!!'
1226	endif
1227	enddo
1228	enddo
1229
1230	! Calcul du flux subsident
1231
1232	do k=2,nlay
1233	do ig=1,ngrid
1234	wqd(ig,k)=fm(ig,k)*q(ig,k)
1235	if (wqd(ig,k).lt.0.) then
1236	! print*,'wqd<0!!!'
1237	endif
1238	enddo
1239	enddo
1240	do ig=1,ngrid
1241	wqd(ig,1)=0.
1242	wqd(ig,nlay+1)=0.
1243	enddo
1244
1245	! Calcul des tendances
1246	do k=1,nlay
1247	do ig=1,ngrid
1248	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1249	& -wqd(ig,k)+wqd(ig,k+1)) &
1250	& *ptimestep/masse(ig,k)
1251	enddo
1252	enddo
1253
1254	endif
1255
1256	therm_tke_max(:,:)=q(:,:)
1257
1258	return
1259	!!! fin nrlmd le 10/04/2012
1260	end
1261

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

Download in other formats:

Original Format