Context Navigation

thermcell_main.F90 @ 2348

Last change on this file since 2348 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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.8 KB

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

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

Download in other formats:

Original Format