Context Navigation

thermcell_main.F90 @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1920:1997 into testing branch
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.2 KB

Line
1	!
2	! $Id: thermcell_main.F90 1999 2014-03-20 09:57:19Z fairhead $
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	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
516	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
517
518	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
519	if (prt_level.ge.10) then
520	write(lunout1,*) 'Dans thermcell_main 1b'
521	write(lunout1,*) 'lmin ',lmin(igout)
522	write(lunout1,*) 'lalim ',lalim(igout)
523	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
524	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
525	& ,l=1,lalim(igout)+4)
526	endif
527
528
529
530
531	! Choix de la fonction d'alimentation utilisee pour la fermeture.
532	! Apparemment sans importance
533	alim_star_clos(:,:)=alim_star(:,:)
534	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
535
536	! Appel avec la version seche
537	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
538	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
539
540	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
541	! Appel avec les zmax et wmax tenant compte de la condensation
542	! Semble moins bien marcher
543	! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
544	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
545	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
546
547	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
548
549	if (tau_thermals>1.) then
550	lambda=exp(-ptimestep/tau_thermals)
551	f0=(1.-lambda)f+lambdaf0
552	else
553	f0=f
554	endif
555
556	! Test valable seulement en 1D mais pas genant
557	if (.not. (f0(1).ge.0.) ) then
558	abort_message = '.not. (f0(1).ge.0.)'
559	CALL abort_gcm (modname,abort_message,1)
560	endif
561
562	!-------------------------------------------------------------------------------
563	!deduction des flux
564	!-------------------------------------------------------------------------------
565
566	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
567	& lalim,lmax,alim_star, &
568	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
569	& detr,zqla,lev_out,lunout1,igout)
570	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
571
572	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
573	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
574	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
575
576	!------------------------------------------------------------------
577	! On ne prend pas directement les profils issus des calculs precedents
578	! mais on s'autorise genereusement une relaxation vers ceci avec
579	! une constante de temps tau_thermals (typiquement 1800s).
580	!------------------------------------------------------------------
581
582	if (tau_thermals>1.) then
583	lambda=exp(-ptimestep/tau_thermals)
584	fm0=(1.-lambda)fm+lambdafm0
585	entr0=(1.-lambda)entr+lambdaentr0
586	detr0=(1.-lambda)detr+lambdadetr0
587	else
588	fm0=fm
589	entr0=entr
590	detr0=detr
591	endif
592
593	!c------------------------------------------------------------------
594	! calcul du transport vertical
595	!------------------------------------------------------------------
596
597	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
598	& zthl,zdthladj,zta,lev_out)
599	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
600	& po,pdoadj,zoa,lev_out)
601
602	!------------------------------------------------------------------
603	! Calcul de la fraction de l'ascendance
604	!------------------------------------------------------------------
605	do ig=1,klon
606	fraca(ig,1)=0.
607	fraca(ig,nlay+1)=0.
608	enddo
609	do l=2,nlay
610	do ig=1,klon
611	if (zw2(ig,l).gt.1.e-10) then
612	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
613	else
614	fraca(ig,l)=0.
615	endif
616	enddo
617	enddo
618
619	!------------------------------------------------------------------
620	! calcul du transport vertical du moment horizontal
621	!------------------------------------------------------------------
622
623	!IM 090508
624	if (dvdq == 0 ) then
625
626	! Calcul du transport de V tenant compte d'echange par gradient
627	! de pression horizontal avec l'environnement
628
629	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
630	! & ,fraca*dvdq,zmax &
631	& ,fraca,zmax &
632	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
633
634	else
635
636	! calcul purement conservatif pour le transport de V
637	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
638	& ,zu,pduadj,zua,lev_out)
639	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
640	& ,zv,pdvadj,zva,lev_out)
641
642	endif
643
644	! print*,'13 OK convect8'
645	do l=1,nlay
646	do ig=1,ngrid
647	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
648	enddo
649	enddo
650
651	if (prt_level.ge.1) print*,'14 OK convect8'
652	!------------------------------------------------------------------
653	! Calculs de diagnostiques pour les sorties
654	!------------------------------------------------------------------
655	!calcul de fraca pour les sorties
656
657	if (sorties) then
658	if (prt_level.ge.1) print*,'14a OK convect8'
659	! calcul du niveau de condensation
660	! initialisation
661	do ig=1,ngrid
662	nivcon(ig)=0
663	zcon(ig)=0.
664	enddo
665	!nouveau calcul
666	do ig=1,ngrid
667	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
668	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
669	enddo
670	!IM do k=1,nlay
671	do k=1,nlay-1
672	do ig=1,ngrid
673	if ((pcon(ig).le.pplay(ig,k)) &
674	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
675	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
676	endif
677	enddo
678	enddo
679	!IM
680	ierr=0
681	do ig=1,ngrid
682	if (pcon(ig).le.pplay(ig,nlay)) then
683	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
684	ierr=1
685	endif
686	enddo
687	if (ierr==1) then
688	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
689	CALL abort_gcm (modname,abort_message,1)
690	endif
691
692	if (prt_level.ge.1) print*,'14b OK convect8'
693	do k=nlay,1,-1
694	do ig=1,ngrid
695	if (zqla(ig,k).gt.1e-10) then
696	nivcon(ig)=k
697	zcon(ig)=zlev(ig,k)
698	endif
699	enddo
700	enddo
701	if (prt_level.ge.1) print*,'14c OK convect8'
702	!calcul des moments
703	!initialisation
704	do l=1,nlay
705	do ig=1,ngrid
706	q2(ig,l)=0.
707	wth2(ig,l)=0.
708	wth3(ig,l)=0.
709	ratqscth(ig,l)=0.
710	ratqsdiff(ig,l)=0.
711	enddo
712	enddo
713	if (prt_level.ge.1) print*,'14d OK convect8'
714	if (prt_level.ge.10)write(lunout,*) &
715	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
716	do l=1,nlay
717	do ig=1,ngrid
718	zf=fraca(ig,l)
719	zf2=zf/(1.-zf)
720	!
721	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
722	if(zw2(ig,l).gt.1.e-10) then
723	wth2(ig,l)=zf2(zw2(ig,l))*2
724	else
725	wth2(ig,l)=0.
726	endif
727	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
728	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
729	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
730	!test: on calcul q2/po=ratqsc
731	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
732	enddo
733	enddo
734	!calcul des flux: q, thetal et thetav
735	do l=1,nlay
736	do ig=1,ngrid
737	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
738	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
739	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
740	enddo
741	enddo
742	!
743
744	!!! nrlmd le 10/04/2012
745
746	!------------Test sur le LCL des thermiques
747	do ig=1,ngrid
748	ok_lcl(ig)=.false.
749	if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
750	enddo
751
752	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
753	do l=1,nlay-1
754	do ig=1,ngrid
755	if (ok_lcl(ig)) then
756	if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
757	klcl(ig)=l
758	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
759	endif
760	endif
761	enddo
762	enddo
763
764	!------------Hauteur des thermiques
765	!!jyg le 27/04/2012
766	!! do ig =1,ngrid
767	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
768	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
769	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
770	!! zmax(ig)=pphi(ig,lmax(ig))/rg
771	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
772	!! enddo
773	do ig =1,ngrid
774	zmax(ig)=pphi(ig,lmax(ig))/rg
775	if (ok_lcl(ig)) then
776	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
777	& -rhobarz(ig,klcl(ig)))*interp(ig)
778	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
779	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
780	else
781	rhobarz0(ig)=0.
782	zlcl(ig)=zmax(ig)
783	endif
784	enddo
785	!!jyg fin
786
787	!------------Calcul des propriétés du thermique au LCL
788	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
789
790	!-----Initialisation de la TKE moyenne
791	do l=1,nlay
792	do ig=1,ngrid
793	pbl_tke_max(ig,l)=0.
794	enddo
795	enddo
796
797	!-----Calcul de la TKE moyenne
798	do nsrf=1,nbsrf
799	do l=1,nlay
800	do ig=1,ngrid
801	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
802	enddo
803	enddo
804	enddo
805
806	!-----Initialisations des TKE dans et hors des thermiques
807	do l=1,nlay
808	do ig=1,ngrid
809	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
810	env_tke_max(ig,l)=pbl_tke_max(ig,l)
811	enddo
812	enddo
813
814	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
815	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
816	& rg,pplev,therm_tke_max)
817	! print *,' thermcell_tke_transport -> ' !!jyg
818
819	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
820	do l=1,nlay
821	do ig=1,ngrid
822	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
823	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
824	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
825	enddo
826	enddo
827	! print *,' apres w_ls = ' !!jyg
828
829	do ig=1,ngrid
830	if (ok_lcl(ig)) then
831	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
832	& -fraca(ig,klcl(ig)))*interp(ig)
833	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
834	& -zw2(ig,klcl(ig)))*interp(ig)
835	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
836	& -w_ls(ig,klcl(ig)))*interp(ig)
837	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
838	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
839	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
840	& -env_tke_max(ig,klcl(ig)))*interp(ig)
841	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
842	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
843	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
844	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
845	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
846	else
847	fraca0(ig)=0.
848	w0(ig)=0.
849	!!jyg le 27/04/2012
850	!! zlcl(ig)=0.
851	!!
852	endif
853	enddo
854
855	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
856	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
857
858	!------------Triggering------------------
859	IF (iflag_trig_bl.ge.1) THEN
860
861	!-----Initialisations
862	depth(:)=0.
863	n2(:)=0.
864	s2(:)=0.
865	s_max(:)=0.
866
867	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
868	do ig=1,ngrid
869	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
870	depth(ig)=zmax_moy(ig)-zlcl(ig)
871	hmin(ig)=hmincoef*zlcl(ig)
872	if (depth(ig).ge.10.) then
873	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
874	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
875	!!
876	!!jyg le 27/04/2012
877	!! s_max(ig)=s2(ig)*log(n2(ig))
878	!! if (n2(ig) .lt. 1) s_max(ig)=0.
879	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
880	!!fin jyg
881	else
882	s2(ig)=0.
883	n2(ig)=0.
884	s_max(ig)=0.
885	endif
886	enddo
887	! print *,'avant Calcul de Wmax ' !!jyg
888
889	!-----Calcul de Wmax et ALE_BL_STAT associée
890	!!jyg le 30/04/2012
891	!! do ig=1,ngrid
892	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
893	!! 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))))
894	!! ale_bl_stat(ig)=0.5w_max(ig)*2
895	!! else
896	!! w_max(ig)=0.
897	!! ale_bl_stat(ig)=0.
898	!! endif
899	!! enddo
900	susqr2pi=susqrt(2.Rpi)
901	Reuler=exp(1.)
902	do ig=1,ngrid
903	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
904	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
905	ale_bl_stat(ig)=0.5w_max(ig)*2
906	else
907	w_max(ig)=0.
908	ale_bl_stat(ig)=0.
909	endif
910	enddo
911
912	ENDIF ! iflag_trig_bl
913	! print *,'ENDIF iflag_trig_bl' !!jyg
914
915	!------------Closure------------------
916
917	IF (iflag_clos_bl.ge.1) THEN
918
919	!-----Calcul de ALP_BL_STAT
920	do ig=1,ngrid
921	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
922	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
923	& (w0(ig)**2)
924	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
925	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
926	if (iflag_clos_bl.ge.2) then
927	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
928	& (w0(ig)**2)
929	else
930	alp_bl_conv(ig)=0.
931	endif
932	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
933	enddo
934
935	!-----Sécurité ALP infinie
936	do ig=1,ngrid
937	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
938	enddo
939
940	ENDIF ! (iflag_clos_bl.ge.1)
941
942	!!! fin nrlmd le 10/04/2012
943
944	if (prt_level.ge.10) then
945	ig=igout
946	do l=1,nlay
947	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
948	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
949	enddo
950	endif
951
952	! print*,'avant calcul ale et alp'
953	!calcul de ALE et ALP pour la convection
954	Alp_bl(:)=0.
955	Ale_bl(:)=0.
956	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
957	do l=1,nlay
958	do ig=1,ngrid
959	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
960	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
961	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
962	enddo
963	enddo
964
965	!test:calcul de la ponderation des couches pour KE
966	!initialisations
967
968	fm_tot(:)=0.
969	wght_th(:,:)=1.
970	lalim_conv(:)=lalim(:)
971
972	do k=1,klev
973	do ig=1,ngrid
974	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
975	enddo
976	enddo
977
978	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
979	! plus petit que 1.e-10, on prend wght_th=1.
980	do k=1,klev
981	do ig=1,ngrid
982	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
983	wght_th(ig,k)=alim_star(ig,k)
984	endif
985	enddo
986	enddo
987
988	! print*,'apres wght_th'
989	!test pour prolonger la convection
990	do ig=1,ngrid
991	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
992	if ((alim_star(ig,1).lt.1.e-10)) then
993	lalim_conv(ig)=1
994	wght_th(ig,1)=1.
995	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
996	endif
997	enddo
998
999	!------------------------------------------------------------------------
1000	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1001	! Integrale verticale de ALP.
1002	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1003	! couches
1004	!------------------------------------------------------------------------
1005
1006	alp_int(:)=0.
1007	dp_int(:)=0.
1008	do l=2,nlay
1009	do ig=1,ngrid
1010	if(l.LE.lmax(ig)) THEN
1011	zdp=pplay(ig,l-1)-pplay(ig,l)
1012	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1013	dp_int(ig)=dp_int(ig)+zdp
1014	endif
1015	enddo
1016	enddo
1017
1018	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1019	do ig=1,ngrid
1020	!valeur integree de alp_bl * 0.5:
1021	if (dp_int(ig)>0.) then
1022	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1023	endif
1024	enddo!
1025	endif
1026
1027
1028	! Facteur multiplicatif sur Alp_bl
1029	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1030
1031	!------------------------------------------------------------------------
1032
1033
1034	!calcul du ratqscdiff
1035	if (prt_level.ge.1) print*,'14e OK convect8'
1036	var=0.
1037	vardiff=0.
1038	ratqsdiff(:,:)=0.
1039
1040	do l=1,klev
1041	do ig=1,ngrid
1042	if (l<=lalim(ig)) then
1043	var=var+alim_star(ig,l)zqta(ig,l)1000.
1044	endif
1045	enddo
1046	enddo
1047
1048	if (prt_level.ge.1) print*,'14f OK convect8'
1049
1050	do l=1,klev
1051	do ig=1,ngrid
1052	if (l<=lalim(ig)) then
1053	zf=fraca(ig,l)
1054	zf2=zf/(1.-zf)
1055	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1056	endif
1057	enddo
1058	enddo
1059
1060	if (prt_level.ge.1) print*,'14g OK convect8'
1061	do l=1,nlay
1062	do ig=1,ngrid
1063	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1064	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1065	enddo
1066	enddo
1067	!--------------------------------------------------------------------
1068	!
1069	!ecriture des fichiers sortie
1070	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1071
1072	endif
1073
1074	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1075
1076	return
1077	end
1078
1079	!-----------------------------------------------------------------------------
1080
1081	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1082	IMPLICIT NONE
1083	#include "iniprint.h"
1084
1085	integer i, k, klon,klev
1086	real pplev(klon,klev+1),pplay(klon,klev)
1087	real ztv(klon,klev)
1088	real po(klon,klev)
1089	real ztva(klon,klev)
1090	real zqla(klon,klev)
1091	real f_star(klon,klev)
1092	real zw2(klon,klev)
1093	integer long(klon)
1094	real seuil
1095	character*21 comment
1096
1097	if (prt_level.ge.1) THEN
1098	print*,'WARNING !!! TEST ',comment
1099	endif
1100	return
1101
1102	! test sur la hauteur des thermiques ...
1103	do i=1,klon
1104	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1105	if (prt_level.ge.10) then
1106	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1107	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1108	do k=1,klev
1109	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)
1110	enddo
1111	endif
1112	enddo
1113
1114
1115	return
1116	end
1117
1118	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1119	! On transporte pbl_tke pour donner therm_tke
1120	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1121	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1122	& rg,pplev,therm_tke_max)
1123	implicit none
1124
1125	#include "iniprint.h"
1126	!=======================================================================
1127	!
1128	! Calcul du transport verticale dans la couche limite en presence
1129	! de "thermiques" explicitement representes
1130	! calcul du dq/dt une fois qu'on connait les ascendances
1131	!
1132	!=======================================================================
1133
1134	integer ngrid,nlay,nsrf
1135
1136	real ptimestep
1137	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1138	real entr0(ngrid,nlay),rg
1139	real therm_tke_max(ngrid,nlay)
1140	real detr0(ngrid,nlay)
1141
1142
1143	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1144	real entr(ngrid,nlay)
1145	real q(ngrid,nlay)
1146	integer lev_out ! niveau pour les print
1147
1148	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1149
1150	real zzm
1151
1152	integer ig,k
1153	integer isrf
1154
1155
1156	lev_out=0
1157
1158
1159	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1160
1161	! calcul du detrainement
1162	do k=1,nlay
1163	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1164	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1165	enddo
1166
1167
1168	! Decalage vertical des entrainements et detrainements.
1169	masse(:,1)=0.5*masse0(:,1)
1170	entr(:,1)=0.5*entr0(:,1)
1171	detr(:,1)=0.5*detr0(:,1)
1172	fm(:,1)=0.
1173	do k=1,nlay-1
1174	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1175	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1176	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1177	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1178	enddo
1179	fm(:,nlay+1)=0.
1180
1181	!!! nrlmd le 16/09/2010
1182	! calcul de la valeur dans les ascendances
1183	! do ig=1,ngrid
1184	! qa(ig,1)=q(ig,1)
1185	! enddo
1186	!!!
1187
1188	!do isrf=1,nsrf
1189
1190	! q(:,:)=therm_tke(:,:,isrf)
1191	q(:,:)=therm_tke_max(:,:)
1192	!!! nrlmd le 16/09/2010
1193	do ig=1,ngrid
1194	qa(ig,1)=q(ig,1)
1195	enddo
1196	!!!
1197
1198	if (1==1) then
1199	do k=2,nlay
1200	do ig=1,ngrid
1201	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1202	& 1.e-5*masse(ig,k)) then
1203	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1204	& /(fm(ig,k+1)+detr(ig,k))
1205	else
1206	qa(ig,k)=q(ig,k)
1207	endif
1208	if (qa(ig,k).lt.0.) then
1209	! print*,'qa<0!!!'
1210	endif
1211	if (q(ig,k).lt.0.) then
1212	! print*,'q<0!!!'
1213	endif
1214	enddo
1215	enddo
1216
1217	! Calcul du flux subsident
1218
1219	do k=2,nlay
1220	do ig=1,ngrid
1221	wqd(ig,k)=fm(ig,k)*q(ig,k)
1222	if (wqd(ig,k).lt.0.) then
1223	! print*,'wqd<0!!!'
1224	endif
1225	enddo
1226	enddo
1227	do ig=1,ngrid
1228	wqd(ig,1)=0.
1229	wqd(ig,nlay+1)=0.
1230	enddo
1231
1232	! Calcul des tendances
1233	do k=1,nlay
1234	do ig=1,ngrid
1235	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1236	& -wqd(ig,k)+wqd(ig,k+1)) &
1237	& *ptimestep/masse(ig,k)
1238	enddo
1239	enddo
1240
1241	endif
1242
1243	therm_tke_max(:,:)=q(:,:)
1244
1245	return
1246	!!! fin nrlmd le 10/04/2012
1247	end
1248

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

Download in other formats:

Original Format