Context Navigation

thermcell_main.F90 @ 1939

Last change on this file since 1939 was 179, checked in by lfita, 10 years ago
No Ale, Alp, ... values because pphi is at the intermediate level!!! Removing raw prints of checkig for Ale, Alp
File size: 42.0 KB

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

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

Download in other formats:

Original Format