Context Navigation

thermcell_main.F90 @ 2351

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

More on physics/dynamics separation and cleanup:

Set things up so that all physics-related initializations are done via iniphysiq.
Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
removed obsolete comgeomphy.h and comgeomphy.F90

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

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

Download in other formats:

Original Format