Context Navigation

thermcell_main.F90 @ 5427

Last change on this file since 5427 was 1790, checked in by Ehouarn Millour, 11 years ago

Ré-implémentation des commissions écrasées par la rev 1785. (Notamment lmdz1d.F)
UG
...................................................

Reimplementation of the commitions undonned by rev 1785. (Especially lmdy1d.F)
UG

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

File size: 40.8 KB

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

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

Download in other formats:

Original Format