Context Navigation

thermcell_main.F90 @ 2196

Last change on this file since 2196 was 1998, checked in by fhourdin, 11 years ago

Reactivation du calcul d'un zmax continu pour les thermiques

(thermcell_height.F90, thermcell_plume.F90, thermcell_dry.F90)
ouvre la voie à la réactivation d'une fermeture humide des thermiques
iflag_thermals_closure=2
(conf_phys_m.F90, thermcell.h, thermcell_main.F90)

Modification liée à la conservation de l'eau

(add_phys_tend.F90, add_pbl_tend.F90, physiq.F90)

Modifications liées au déclenchement stochastique

possibilité de revenir à la Ale déterministe pour le criter ALE>|CIN| iflag_trig_bl=2, 1 par défaut)
possibilité d'activer une fermeture statistique où ALP est divisé par la probabilité de déclenchement iflag_clos_bl=1 (0 par défaut, ancienne option 1 passée en =2)

Modification de l'entrainemement dans la version "stratocumulus" du

modèle du thermique (quand iflag_thermals_ed=8).
(modifie thermcell_plume.F90)

Catherine, Jean-Yves 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.6 KB

Line
1	!
2	! $Id: thermcell_main.F90 1998 2014-03-19 14:14:30Z jescribano $
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(:)=0.
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	s2(ig)=0.
897	n2(ig)=0.
898	s_max(ig)=0.
899	endif
900	enddo
901	! print *,'avant Calcul de Wmax ' !!jyg
902
903	!-----Calcul de Wmax et ALE_BL_STAT associée
904	!!jyg le 30/04/2012
905	!! do ig=1,ngrid
906	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
907	!! 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))))
908	!! ale_bl_stat(ig)=0.5w_max(ig)*2
909	!! else
910	!! w_max(ig)=0.
911	!! ale_bl_stat(ig)=0.
912	!! endif
913	!! enddo
914	susqr2pi=susqrt(2.Rpi)
915	Reuler=exp(1.)
916	do ig=1,ngrid
917	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
918	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
919	ale_bl_stat(ig)=0.5w_max(ig)*2
920	else
921	w_max(ig)=0.
922	ale_bl_stat(ig)=0.
923	endif
924	enddo
925
926	ENDIF ! iflag_trig_bl
927	! print *,'ENDIF iflag_trig_bl' !!jyg
928
929	!------------Closure------------------
930
931	IF (iflag_clos_bl.ge.2) THEN
932
933	!-----Calcul de ALP_BL_STAT
934	do ig=1,ngrid
935	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
936	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
937	& (w0(ig)**2)
938	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
939	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
940	if (iflag_clos_bl.ge.2) then
941	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
942	& (w0(ig)**2)
943	else
944	alp_bl_conv(ig)=0.
945	endif
946	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
947	enddo
948
949	!-----Sécurité ALP infinie
950	do ig=1,ngrid
951	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
952	enddo
953
954	ENDIF ! (iflag_clos_bl.ge.2)
955
956	!!! fin nrlmd le 10/04/2012
957
958	if (prt_level.ge.10) then
959	ig=igout
960	do l=1,nlay
961	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
962	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
963	enddo
964	endif
965
966	! print*,'avant calcul ale et alp'
967	!calcul de ALE et ALP pour la convection
968	Alp_bl(:)=0.
969	Ale_bl(:)=0.
970	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
971	do l=1,nlay
972	do ig=1,ngrid
973	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
974	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
975	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
976	enddo
977	enddo
978
979	!test:calcul de la ponderation des couches pour KE
980	!initialisations
981
982	fm_tot(:)=0.
983	wght_th(:,:)=1.
984	lalim_conv(:)=lalim(:)
985
986	do k=1,klev
987	do ig=1,ngrid
988	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
989	enddo
990	enddo
991
992	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
993	! plus petit que 1.e-10, on prend wght_th=1.
994	do k=1,klev
995	do ig=1,ngrid
996	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
997	wght_th(ig,k)=alim_star(ig,k)
998	endif
999	enddo
1000	enddo
1001
1002	! print*,'apres wght_th'
1003	!test pour prolonger la convection
1004	do ig=1,ngrid
1005	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1006	if ((alim_star(ig,1).lt.1.e-10)) then
1007	lalim_conv(ig)=1
1008	wght_th(ig,1)=1.
1009	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1010	endif
1011	enddo
1012
1013	!------------------------------------------------------------------------
1014	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1015	! Integrale verticale de ALP.
1016	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1017	! couches
1018	!------------------------------------------------------------------------
1019
1020	alp_int(:)=0.
1021	dp_int(:)=0.
1022	do l=2,nlay
1023	do ig=1,ngrid
1024	if(l.LE.lmax(ig)) THEN
1025	zdp=pplay(ig,l-1)-pplay(ig,l)
1026	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1027	dp_int(ig)=dp_int(ig)+zdp
1028	endif
1029	enddo
1030	enddo
1031
1032	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1033	do ig=1,ngrid
1034	!valeur integree de alp_bl * 0.5:
1035	if (dp_int(ig)>0.) then
1036	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1037	endif
1038	enddo!
1039	endif
1040
1041
1042	! Facteur multiplicatif sur Alp_bl
1043	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1044
1045	!------------------------------------------------------------------------
1046
1047
1048	!calcul du ratqscdiff
1049	if (prt_level.ge.1) print*,'14e OK convect8'
1050	var=0.
1051	vardiff=0.
1052	ratqsdiff(:,:)=0.
1053
1054	do l=1,klev
1055	do ig=1,ngrid
1056	if (l<=lalim(ig)) then
1057	var=var+alim_star(ig,l)zqta(ig,l)1000.
1058	endif
1059	enddo
1060	enddo
1061
1062	if (prt_level.ge.1) print*,'14f OK convect8'
1063
1064	do l=1,klev
1065	do ig=1,ngrid
1066	if (l<=lalim(ig)) then
1067	zf=fraca(ig,l)
1068	zf2=zf/(1.-zf)
1069	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1070	endif
1071	enddo
1072	enddo
1073
1074	if (prt_level.ge.1) print*,'14g OK convect8'
1075	do l=1,nlay
1076	do ig=1,ngrid
1077	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1078	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1079	enddo
1080	enddo
1081	!--------------------------------------------------------------------
1082	!
1083	!ecriture des fichiers sortie
1084	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1085
1086	endif
1087
1088	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1089
1090	return
1091	end
1092
1093	!-----------------------------------------------------------------------------
1094
1095	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1096	IMPLICIT NONE
1097	#include "iniprint.h"
1098
1099	integer i, k, klon,klev
1100	real pplev(klon,klev+1),pplay(klon,klev)
1101	real ztv(klon,klev)
1102	real po(klon,klev)
1103	real ztva(klon,klev)
1104	real zqla(klon,klev)
1105	real f_star(klon,klev)
1106	real zw2(klon,klev)
1107	integer long(klon)
1108	real seuil
1109	character*21 comment
1110
1111	if (prt_level.ge.1) THEN
1112	print*,'WARNING !!! TEST ',comment
1113	endif
1114	return
1115
1116	! test sur la hauteur des thermiques ...
1117	do i=1,klon
1118	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1119	if (prt_level.ge.10) then
1120	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1121	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1122	do k=1,klev
1123	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)
1124	enddo
1125	endif
1126	enddo
1127
1128
1129	return
1130	end
1131
1132	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1133	! On transporte pbl_tke pour donner therm_tke
1134	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1135	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1136	& rg,pplev,therm_tke_max)
1137	implicit none
1138
1139	#include "iniprint.h"
1140	!=======================================================================
1141	!
1142	! Calcul du transport verticale dans la couche limite en presence
1143	! de "thermiques" explicitement representes
1144	! calcul du dq/dt une fois qu'on connait les ascendances
1145	!
1146	!=======================================================================
1147
1148	integer ngrid,nlay,nsrf
1149
1150	real ptimestep
1151	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1152	real entr0(ngrid,nlay),rg
1153	real therm_tke_max(ngrid,nlay)
1154	real detr0(ngrid,nlay)
1155
1156
1157	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1158	real entr(ngrid,nlay)
1159	real q(ngrid,nlay)
1160	integer lev_out ! niveau pour les print
1161
1162	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1163
1164	real zzm
1165
1166	integer ig,k
1167	integer isrf
1168
1169
1170	lev_out=0
1171
1172
1173	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1174
1175	! calcul du detrainement
1176	do k=1,nlay
1177	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1178	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1179	enddo
1180
1181
1182	! Decalage vertical des entrainements et detrainements.
1183	masse(:,1)=0.5*masse0(:,1)
1184	entr(:,1)=0.5*entr0(:,1)
1185	detr(:,1)=0.5*detr0(:,1)
1186	fm(:,1)=0.
1187	do k=1,nlay-1
1188	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1189	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1190	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1191	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1192	enddo
1193	fm(:,nlay+1)=0.
1194
1195	!!! nrlmd le 16/09/2010
1196	! calcul de la valeur dans les ascendances
1197	! do ig=1,ngrid
1198	! qa(ig,1)=q(ig,1)
1199	! enddo
1200	!!!
1201
1202	!do isrf=1,nsrf
1203
1204	! q(:,:)=therm_tke(:,:,isrf)
1205	q(:,:)=therm_tke_max(:,:)
1206	!!! nrlmd le 16/09/2010
1207	do ig=1,ngrid
1208	qa(ig,1)=q(ig,1)
1209	enddo
1210	!!!
1211
1212	if (1==1) then
1213	do k=2,nlay
1214	do ig=1,ngrid
1215	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1216	& 1.e-5*masse(ig,k)) then
1217	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1218	& /(fm(ig,k+1)+detr(ig,k))
1219	else
1220	qa(ig,k)=q(ig,k)
1221	endif
1222	if (qa(ig,k).lt.0.) then
1223	! print*,'qa<0!!!'
1224	endif
1225	if (q(ig,k).lt.0.) then
1226	! print*,'q<0!!!'
1227	endif
1228	enddo
1229	enddo
1230
1231	! Calcul du flux subsident
1232
1233	do k=2,nlay
1234	do ig=1,ngrid
1235	wqd(ig,k)=fm(ig,k)*q(ig,k)
1236	if (wqd(ig,k).lt.0.) then
1237	! print*,'wqd<0!!!'
1238	endif
1239	enddo
1240	enddo
1241	do ig=1,ngrid
1242	wqd(ig,1)=0.
1243	wqd(ig,nlay+1)=0.
1244	enddo
1245
1246	! Calcul des tendances
1247	do k=1,nlay
1248	do ig=1,ngrid
1249	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1250	& -wqd(ig,k)+wqd(ig,k+1)) &
1251	& *ptimestep/masse(ig,k)
1252	enddo
1253	enddo
1254
1255	endif
1256
1257	therm_tke_max(:,:)=q(:,:)
1258
1259	return
1260	!!! fin nrlmd le 10/04/2012
1261	end
1262

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

Download in other formats:

Original Format