Context Navigation

thermcell_main.F90 @ 1738

Last change on this file since 1738 was 1738, checked in by idelkadi, 11 years ago

Modifications for numerical stability of the boundary layer parameterizations.
Concerning yamada4.F :

option with asymptotic mixing length l0 imposed (iflag_pbl=8 and 9)
option with a new temporal scheme (iflag_pbl=10 and 11)
Correction for very stable PBLs (iflag_pbl=10 and 11) iflag_pbl=8 converges numerically with NPv3.1 iflag_pbl=11 -> now starts with NP from start files created by ce0l

-> the model can run with longer time-steps.

Concerning thermals :

Introduction of an implicit computation of vertical advection in
the environment of thermal plumes in thermcell_dq
impl = 0 : explicit, 1 : implicit, -1 : old version
controled by iflag_thermals =

15, 16 run with impl=-1 : numerical convergence with NPv3
17, 18 run with impl=1 : more stable

15 and 17 correspond to the activation of the stratocumulus "bidouille"

Modified routines (phylmd/):
calltherm.F90 : for managing the various options of thermcell_dq
coef_diff_turb_mod.F90 : yamada4 called for iflag_pbl<= 18 instead of 11
physiq.F : desactivation of the vertical diffusion of TKE for iflag_pbl=10
thermcellV0_main.F90 : calling thermcell_dq with implicit scheme
thermcell_dq.F90 : thermcell_dq with optional explicit/implicit scheme
thermcell_main.F90 : various option for vertical transport by thermals
yamada4.F : Yamada scheme with new options

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 40.9 KB

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

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

Download in other formats:

Original Format