Context Navigation

thermcell_main.F90 @ 132

Last change on this file since 132 was 132, checked in by lfita, 10 years ago
Adding checks for rhobarz0
File size: 54.7 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	! Lluis
254	INTEGER :: llp
255	CHARACTER(LEN=50) :: lvarname, lfname
256	REAL :: largest
257
258	llp = 734
259	lfname = 'physiq'
260	largest = 10.e5
261
262	! L. Fita, LMD July 2014. Initializing variables.
263	! Some not initializated according to values: iflag_trig_bl, iflag_clos_bl
264	zdthladj = 0.
265	pbl_tke_max0 = 0.
266	fraca0 = 0.
267	w_conv = 0.
268	w0 = 0.
269	therm_tke_max0 = 0.
270	env_tke_max0 = 0.
271	alp_bl_det = 0.
272	alp_bl_fluct_m = 0.
273	alp_bl_fluct_tke = 0.
274	alp_bl_conv = 0.
275	alp_bl_stat = 0.
276	interp = 0.
277	klcl = 0
278
279	!-----------------------------------------------------------------------
280	! initialisation:
281	! ---------------
282	!
283
284	seuil=0.25
285
286	if (debut) then
287	! call getin('dvdq',dvdq)
288	! call getin('dqimpl',dqimpl)
289
290	if (iflag_thermals==15.or.iflag_thermals==16) then
291	dvdq=0
292	dqimpl=-1
293	else
294	dvdq=1
295	dqimpl=1
296	endif
297
298	fm0=0.
299	entr0=0.
300	detr0=0.
301
302
303	#undef wrgrads_thermcell
304	#ifdef wrgrads_thermcell
305	! Initialisation des sorties grads pour les thermiques.
306	! Pour l'instant en 1D sur le point igout.
307	! Utilise par thermcell_out3d.h
308	str10='therm'
309	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
310	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
311	& ptimestep,str10,'therm ')
312	#endif
313
314
315
316	endif
317
318	fm=0. ; entr=0. ; detr=0.
319
320
321	icount=icount+1
322
323	!IM 090508 beg
324	!print*,'====================================================================='
325	!print*,'====================================================================='
326	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
327	!print*,'====================================================================='
328	!print*,'====================================================================='
329	!IM 090508 end
330
331	if (prt_level.ge.1) print*,'thermcell_main V4'
332
333	sorties=.true.
334	IF(ngrid.NE.ngrid) THEN
335	PRINT*
336	PRINT*,'STOP dans convadj'
337	PRINT*,'ngrid =',ngrid
338	PRINT*,'ngrid =',ngrid
339	ENDIF
340	!
341	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
342	do ig=1,ngrid
343	f0(ig)=max(f0(ig),1.e-2)
344	zmax0(ig)=max(zmax0(ig),40.)
345	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
346	enddo
347
348	if (prt_level.ge.20) then
349	do ig=1,ngrid
350	print*,'th_main ig f0',ig,f0(ig)
351	enddo
352	endif
353	!-----------------------------------------------------------------------
354	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
355	! --------------------------------------------------------------------
356	!
357	lfname='thermcell_main before thermcell_env'
358	lvarname = 'pt'
359	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
360	lvarname = 'pdtadj'
361	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
362	lvarname = 'pplev'
363	CALL check_var3D(lfname, lvarname, pplev, ngrid, nlay, largest, .FALSE.)
364
365	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
366	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
367
368	lfname='thermcell_main after thermcell_env'
369	lvarname = 'pt'
370	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
371	lvarname = 'pdtadj'
372	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
373	lvarname = 'pplev'
374	CALL check_var3D(lfname, lvarname, pplev, ngrid, nlay, largest, .FALSE.)
375
376	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
377
378	!------------------------------------------------------------------------
379	! --------------------
380	!
381	!
382	! + + + + + + + + + + +
383	!
384	!
385	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
386	! wh,wt,wo ...
387	!
388	! + + + + + + + + + + + zh,zu,zv,zo,rho
389	!
390	!
391	! -------------------- zlev(1)
392	! \\\\\\\\\\\\\\\\\\\\
393	!
394	!
395
396	!-----------------------------------------------------------------------
397	! Calcul des altitudes des couches
398	!-----------------------------------------------------------------------
399
400	do l=2,nlay
401	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
402	enddo
403	zlev(:,1)=0.
404	zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
405	do l=1,nlay
406	zlay(:,l)=pphi(:,l)/RG
407	enddo
408	!calcul de l epaisseur des couches
409	do l=1,nlay
410	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
411	enddo
412
413	! print*,'2 OK convect8'
414	!-----------------------------------------------------------------------
415	! Calcul des densites
416	!-----------------------------------------------------------------------
417
418	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
419
420	if (prt_level.ge.10)write(lunout,*) &
421	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
422	rhobarz(:,1)=rho(:,1)
423
424	do l=2,nlay
425	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
426	enddo
427
428	!calcul de la masse
429	do l=1,nlay
430	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
431	enddo
432
433	lfname='thermcell_main after initialization'
434	lvarname = 'rhobarz'
435	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
436	lvarname = 'rho'
437	CALL check_var3D(lfname, lvarname, rho, ngrid, nlay, largest, .FALSE.)
438	lvarname = 'zpspsk'
439	CALL check_var3D(lfname, lvarname, zpspsk, ngrid, nlay, largest, .FALSE.)
440	lvarname = 'pplay'
441	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest, .FALSE.)
442	lvarname = 'zw2'
443	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
444
445	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
446
447	!------------------------------------------------------------------
448	!
449	! /\|\
450	! -------- \| F_k+1 -------
451	! ----> D_k
452	! /\|\ <---- E_k , A_k
453	! -------- \| F_k ---------
454	! ----> D_k-1
455	! <---- E_k-1 , A_k-1
456	!
457	!
458	!
459	!
460	!
461	! ---------------------------
462	!
463	! ----- F_lmax+1=0 ---------- \
464	! lmax (zmax) \|
465	! --------------------------- \|
466	! \|
467	! --------------------------- \|
468	! \|
469	! --------------------------- \|
470	! \|
471	! --------------------------- \|
472	! \|
473	! --------------------------- \|
474	! \| E
475	! --------------------------- \| D
476	! \|
477	! --------------------------- \|
478	! \|
479	! --------------------------- \ \|
480	! lalim \| \|
481	! --------------------------- \| \|
482	! \| \|
483	! --------------------------- \| \|
484	! \| A \|
485	! --------------------------- \| \|
486	! \| \|
487	! --------------------------- \| \|
488	! lmin (=1 pour le moment) \| \|
489	! ----- F_lmin=0 ------------ / /
490	!
491	! ---------------------------
492	! //////////////////////////
493	!
494	!
495	!=============================================================================
496	! Calculs initiaux ne faisant pas intervenir les changements de phase
497	!=============================================================================
498
499	!------------------------------------------------------------------
500	! 1. alim_star est le profil vertical de l'alimentation a la base du
501	! panache thermique, calcule a partir de la flotabilite de l'air sec
502	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
503	!------------------------------------------------------------------
504	!
505	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
506	lmin=1
507
508	!-----------------------------------------------------------------------------
509	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
510	! panache sec conservatif (e=d=0) alimente selon alim_star
511	! Il s'agit d'un calcul de type CAPE
512	! zmax_sec est utilise pour determiner la geometrie du thermique.
513	!------------------------------------------------------------------------------
514	!---------------------------------------------------------------------------------
515	!calcul du melange et des variables dans le thermique
516	!--------------------------------------------------------------------------------
517	!
518	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
519	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
520
521	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
522	! de la variable iflag_thermals
523	lfname='thermcell_main before plume'
524	lvarname = 'alim_star'
525	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
526
527	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
528	if (iflag_thermals_ed<=9) then
529	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
530	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
531	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
532	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
533	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
534	& ,lev_out,lunout1,igout)
535
536	elseif (iflag_thermals_ed>9) then
537	! print*,'THERM RIO et al 2010, version d Arnaud'
538	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
539	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
540	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
541	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
542	& ,lev_out,lunout1,igout)
543
544	endif
545	lfname='thermcell_main after thermcell_plume'
546	lvarname = 'pt'
547	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
548	lvarname = 'pdtadj'
549	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
550	lvarname = 'zw2'
551	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
552	lvarname = 'alim_star'
553	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
554
555	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
556
557	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
558	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
559
560	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
561	if (prt_level.ge.10) then
562	write(lunout1,*) 'Dans thermcell_main 2'
563	write(lunout1,*) 'lmin ',lmin(igout)
564	write(lunout1,*) 'lalim ',lalim(igout)
565	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
566	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
567	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
568	endif
569
570	!-------------------------------------------------------------------------------
571	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
572	!-------------------------------------------------------------------------------
573	!
574	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
575	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
576	! Attention, w2 est transforme en sa racine carree dans cette routine
577	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
578	wmax_tmp=0.
579	do l=1,nlay
580	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
581	enddo
582	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
583
584	lfname='thermcell_main after thermcell_height'
585	lvarname = 'zw2'
586	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
587
588	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
589	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
590	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
591	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
592
593	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
594
595	!-------------------------------------------------------------------------------
596	! Fermeture,determination de f
597	!-------------------------------------------------------------------------------
598	!
599	!
600	!! write(lunout,*)'THERM NOUVEAU XXXXX'
601	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
602	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
603
604	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
605	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
606
607	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
608	if (prt_level.ge.10) then
609	write(lunout1,*) 'Dans thermcell_main 1b'
610	write(lunout1,*) 'lmin ',lmin(igout)
611	write(lunout1,*) 'lalim ',lalim(igout)
612	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
613	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
614	& ,l=1,lalim(igout)+4)
615	endif
616
617	lfname='thermcell_main after thermcell_dry'
618	lvarname = 'alim_star'
619	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
620
621
622
623	! Choix de la fonction d'alimentation utilisee pour la fermeture.
624	! Apparemment sans importance
625	alim_star_clos(:,:)=alim_star(:,:)
626	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
627
628	! Appel avec la version seche
629	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
630	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
631
632	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
633	! Appel avec les zmax et wmax tenant compte de la condensation
634	! Semble moins bien marcher
635	! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
636	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
637	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
638
639	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
640
641	if (tau_thermals>1.) then
642	lambda=exp(-ptimestep/tau_thermals)
643	f0=(1.-lambda)f+lambdaf0
644	else
645	f0=f
646	endif
647
648	! Test valable seulement en 1D mais pas genant
649	if (.not. (f0(1).ge.0.) ) then
650	abort_message = '.not. (f0(1).ge.0.)'
651	CALL abort_gcm (modname,abort_message,1)
652	endif
653
654	!-------------------------------------------------------------------------------
655	!deduction des flux
656	!-------------------------------------------------------------------------------
657	lfname='thermcell_main before flux'
658	lvarname = 'fm'
659	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
660	lvarname = 'entr'
661	CALL check_var3D(lfname, lvarname, entr, ngrid, nlay, largest, .FALSE.)
662	lvarname = 'detr'
663	CALL check_var3D(lfname, lvarname, detr, ngrid, nlay, largest, .FALSE.)
664	lvarname = 'zqla'
665	CALL check_var3D(lfname, lvarname, zqla, ngrid, nlay, largest, .FALSE.)
666	lvarname = 'zw2'
667	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
668	lvarname = 'alim_star'
669	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
670
671	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
672	& lalim,lmax,alim_star, &
673	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
674	& detr,zqla,lev_out,lunout1,igout)
675	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
676	lfname='thermcell_main after flux'
677	lvarname = 'fm'
678	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
679	lvarname = 'zw2'
680	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
681	lvarname = 'alim_star'
682	CALL check_var3D(lfname, lvarname, alim_star, ngrid, nlay, largest, .FALSE.)
683
684	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
685	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
686	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
687
688	!------------------------------------------------------------------
689	! On ne prend pas directement les profils issus des calculs precedents
690	! mais on s'autorise genereusement une relaxation vers ceci avec
691	! une constante de temps tau_thermals (typiquement 1800s).
692	!------------------------------------------------------------------
693
694	if (tau_thermals>1.) then
695	lambda=exp(-ptimestep/tau_thermals)
696	fm0=(1.-lambda)fm+lambdafm0
697	entr0=(1.-lambda)entr+lambdaentr0
698	detr0=(1.-lambda)detr+lambdadetr0
699	else
700	fm0=fm
701	entr0=entr
702	detr0=detr
703	endif
704
705	!c------------------------------------------------------------------
706	! calcul du transport vertical
707	!------------------------------------------------------------------
708
709	lfname='thermcell_main before transport_vertical'
710	lvarname = 'zdthladj'
711	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
712
713	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
714	& zthl,zdthladj,zta,lev_out)
715	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
716	& po,pdoadj,zoa,lev_out)
717	lfname='thermcell_main after transport_vertical'
718	lvarname = 'zdthladj'
719	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
720	lvarname = 'masse'
721	CALL check_var3D(lfname, lvarname, masse, ngrid, nlay, largest, .FALSE.)
722
723	lfname='thermcell_main before fraction ascendance'
724	lvarname = 'fm'
725	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
726	lvarname = 'zw2'
727	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
728
729	!------------------------------------------------------------------
730	! Calcul de la fraction de l'ascendance
731	!------------------------------------------------------------------
732	do ig=1,ngrid
733	fraca(ig,1)=0.
734	fraca(ig,nlay+1)=0.
735	enddo
736	do l=2,nlay
737	do ig=1,ngrid
738	if (zw2(ig,l).gt.1.e-10) then
739	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
740	else
741	fraca(ig,l)=0.
742	endif
743	enddo
744	enddo
745
746	!------------------------------------------------------------------
747	! calcul du transport vertical du moment horizontal
748	!------------------------------------------------------------------
749	lfname='before thermcell_dv2'
750	lvarname = 'pt'
751	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
752	lvarname = 'pdtadj'
753	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
754	lvarname = 'zw2'
755	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay+1, largest, .FALSE.)
756	lvarname = 'rhobarz'
757	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
758	lvarname = 'fraca'
759	CALL check_var3D(lfname, lvarname, fraca, ngrid, nlay+1, largest, .FALSE.)
760
761	!IM 090508
762	if (dvdq == 0 ) then
763
764	! Calcul du transport de V tenant compte d'echange par gradient
765	! de pression horizontal avec l'environnement
766
767	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
768	! & ,fraca*dvdq,zmax &
769	& ,fraca,zmax &
770	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
771
772	else
773
774	! calcul purement conservatif pour le transport de V
775	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
776	& ,zu,pduadj,zua,lev_out)
777	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
778	& ,zv,pdvadj,zva,lev_out)
779
780	endif
781
782	! print*,'13 OK convect8'
783	do l=1,nlay
784	do ig=1,ngrid
785	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
786	enddo
787	enddo
788	lfname='after thermcell_dv2'
789	lvarname = 'pt'
790	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
791	lvarname = 'pdtadj'
792	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
793	lvarname = 'zdthladj'
794	CALL check_var3D(lfname, lvarname, zdthladj, ngrid, nlay, largest, .FALSE.)
795	lvarname = 'zpspsk'
796	CALL check_var3D(lfname, lvarname, zpspsk, ngrid, nlay, largest, .FALSE.)
797	lvarname = 'pduadj'
798	CALL check_var3D(lfname, lvarname, pduadj, ngrid, nlay, largest, .FALSE.)
799	lvarname = 'pdvadj'
800	CALL check_var3D(lfname, lvarname, pdvadj, ngrid, nlay, largest, .FALSE.)
801
802	if (prt_level.ge.1) print*,'14 OK convect8'
803	!------------------------------------------------------------------
804	! Calculs de diagnostiques pour les sorties
805	!------------------------------------------------------------------
806	!calcul de fraca pour les sorties
807
808	if (sorties) then
809	if (prt_level.ge.1) print*,'14a OK convect8'
810	! calcul du niveau de condensation
811	! initialisation
812	do ig=1,ngrid
813	nivcon(ig)=0
814	zcon(ig)=0.
815	enddo
816	!nouveau calcul
817	do ig=1,ngrid
818	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
819	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
820	enddo
821	!IM do k=1,nlay
822	do k=1,nlay-1
823	do ig=1,ngrid
824	if ((pcon(ig).le.pplay(ig,k)) &
825	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
826	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
827	endif
828	enddo
829	enddo
830	!IM
831	ierr=0
832	do ig=1,ngrid
833	if (pcon(ig).le.pplay(ig,nlay)) then
834	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
835	ierr=1
836	endif
837	enddo
838	if (ierr==1) then
839	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
840	CALL abort_gcm (modname,abort_message,1)
841	endif
842
843	if (prt_level.ge.1) print*,'14b OK convect8'
844	do k=nlay,1,-1
845	do ig=1,ngrid
846	if (zqla(ig,k).gt.1e-10) then
847	nivcon(ig)=k
848	zcon(ig)=zlev(ig,k)
849	endif
850	enddo
851	enddo
852	if (prt_level.ge.1) print*,'14c OK convect8'
853	!calcul des moments
854	!initialisation
855	do l=1,nlay
856	do ig=1,ngrid
857	q2(ig,l)=0.
858	wth2(ig,l)=0.
859	wth3(ig,l)=0.
860	ratqscth(ig,l)=0.
861	ratqsdiff(ig,l)=0.
862	enddo
863	enddo
864	if (prt_level.ge.1) print*,'14d OK convect8'
865	if (prt_level.ge.10)write(lunout,*) &
866	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
867	do l=1,nlay
868	do ig=1,ngrid
869	zf=fraca(ig,l)
870	zf2=zf/(1.-zf)
871	!
872	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
873	if(zw2(ig,l).gt.1.e-10) then
874	wth2(ig,l)=zf2(zw2(ig,l))*2
875	else
876	wth2(ig,l)=0.
877	endif
878	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
879	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
880	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
881	!test: on calcul q2/po=ratqsc
882	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
883	enddo
884	enddo
885	lfname='thermcell_main calculation of wth3'
886	lvarname = 'wth3'
887	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
888	lvarname = 'fraca'
889	CALL check_var3D(lfname, lvarname, fraca, ngrid, nlay, largest, .FALSE.)
890	lvarname = 'zw2'
891	CALL check_var3D(lfname, lvarname, zw2, ngrid, nlay, largest, .FALSE.)
892	lvarname = '1-fraca'
893	CALL check_var3D(lfname, lvarname, 1./(1.-fraca), ngrid, nlay, largest, .FALSE.)
894	!calcul des flux: q, thetal et thetav
895	do l=1,nlay
896	do ig=1,ngrid
897	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
898	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
899	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
900	enddo
901	enddo
902	!
903
904	!!! nrlmd le 10/04/2012
905
906	!------------Test sur le LCL des thermiques
907	do ig=1,ngrid
908	ok_lcl(ig)=.false.
909	if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
910	enddo
911
912	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
913	do l=1,nlay-1
914	do ig=1,ngrid
915	if (ok_lcl(ig)) then
916	! if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
917	! L. Fita, LMD July 2014. Adding avoiding divisons by zero...
918	if ((pplay(ig,l).ge.pcon(ig)) .and. (pplay(ig,l+1).le.pcon(ig)) .and. &
919	(klcl(ig).gt.0).and.(klcl(ig)+1.le.nlay-1) .and. (klcl(ig)+1.gt.0) ) then
920	! (pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)).ne.0.) ) then
921	klcl(ig)=l
922	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
923	IF (interp(ig) /= interp(ig)) THEN
924	PRINT *,' Lluis wrong interp= ',interp(ig),' at ', ig,' klcl(ig): ', &
925	klcl(ig),' klcl(ig)+1: ', klcl(ig)+1
926	END IF
927	endif
928	endif
929	enddo
930	enddo
931	lfname='thermcell_main calculation of LCL'
932	lvarname = 'interp'
933	CALL check_var3D(lfname, lvarname, interp, ngrid, nlay, largest, .FALSE.)
934	lvarname = 'klcl'
935	CALL check_var(lfname, lvarname, REAL(klcl), ngrid, largest, .FALSE.)
936	lvarname = 'pcon'
937	CALL check_var(lfname, lvarname, pcon, ngrid, largest, .FALSE.)
938	lvarname = 'pplay'
939	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest, .FALSE.)
940	lvarname = '1/pplay'
941	CALL check_var3D(lfname, lvarname, 1./pplay, ngrid, nlay, largest, .FALSE.)
942
943
944	!------------Hauteur des thermiques
945	!!jyg le 27/04/2012
946	!! do ig =1,ngrid
947	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
948	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
949	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
950	!! zmax(ig)=pphi(ig,lmax(ig))/rg
951	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
952	!! enddo
953	do ig =1,ngrid
954	zmax(ig)=pphi(ig,lmax(ig))/rg
955	! if (ok_lcl(ig)) then
956	! L. Fita, LMD July 2014. Adding avoiding divisons by zero...
957	if ( ok_lcl(ig) .and. (klcl(ig).gt.0) .and. (klcl(ig)+1.le.nlay-1) .and. &
958	(klcl(ig)+1.gt.0) ) then
959	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
960	& -rhobarz(ig,klcl(ig)))*interp(ig)
961	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
962	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
963	else
964	rhobarz0(ig)=0.
965	zlcl(ig)=zmax(ig)
966	endif
967	enddo
968	!!jyg fin
969	lfname='thermcell_main before LCL'
970	lvarname = 'rhobarz0'
971	CALL check_var(lfname, lvarname, rhobarz0, ngrid, largest, .FALSE.)
972	lvarname = 'pcon'
973	CALL check_var(lfname, lvarname, pcon, ngrid, largest, .FALSE.)
974	lvarname = 'fraca'
975	CALL check_var3D(lfname, lvarname, fraca, ngrid, nlay+1, largest, .FALSE.)
976	lvarname = 'fraca0'
977	CALL check_var(lfname, lvarname, fraca0, ngrid, largest, .FALSE.)
978	lvarname = 'w_conv'
979	CALL check_var(lfname, lvarname, w_conv, ngrid, largest, .FALSE.)
980	lvarname = 'interp'
981	CALL check_var(lfname, lvarname, interp, ngrid, largest, .FALSE.)
982	lvarname = 'zlcl'
983	CALL check_var(lfname, lvarname, zlcl, ngrid, largest, .FALSE.)
984
985	!------------Calcul des propriétés du thermique au LCL
986	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
987
988	!-----Initialisation de la TKE moyenne
989	do l=1,nlay
990	do ig=1,ngrid
991	pbl_tke_max(ig,l)=0.
992	enddo
993	enddo
994
995	!-----Calcul de la TKE moyenne
996	do nsrf=1,nbsrf
997	do l=1,nlay
998	do ig=1,ngrid
999	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
1000	enddo
1001	enddo
1002	enddo
1003
1004	!-----Initialisations des TKE dans et hors des thermiques
1005	do l=1,nlay
1006	do ig=1,ngrid
1007	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
1008	env_tke_max(ig,l)=pbl_tke_max(ig,l)
1009	enddo
1010	enddo
1011
1012	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
1013	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1014	& rg,pplev,therm_tke_max)
1015	! print *,' thermcell_tke_transport -> ' !!jyg
1016
1017	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
1018	do l=1,nlay
1019	do ig=1,ngrid
1020	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
1021	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
1022	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
1023	enddo
1024	enddo
1025	! print *,' apres w_ls = ' !!jyg
1026
1027	do ig=1,ngrid
1028	if (ok_lcl(ig)) then
1029	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
1030	& -fraca(ig,klcl(ig)))*interp(ig)
1031	IF (fraca0(ig) /= fraca0(ig) .OR. ABS(fraca0(ig)) > largest*10.e5) THEN
1032	PRINT *,' Lluis wrong fraca0(ig): ',fraca0(ig),' at : ',ig
1033	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1034	' fraca(ig,klcl(ig)): ',fraca(ig,klcl(ig)),' fraca(ig,klcl(ig)+1): ', &
1035	fraca(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1036	END IF
1037
1038	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
1039	& -zw2(ig,klcl(ig)))*interp(ig)
1040	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
1041	& -w_ls(ig,klcl(ig)))*interp(ig)
1042	IF (w_conv(ig) /= w_conv(ig) .OR. ABS(w_conv(ig)) > largest*10.e5) THEN
1043	PRINT *,' Lluis wrong w_conv(ig): ',w_conv(ig),' at : ',ig
1044	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1045	' w_ls(ig,klcl(ig)): ',w_ls(ig,klcl(ig)),' w_ls(ig,klcl(ig)+1): ', &
1046	w_ls(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1047	END IF
1048	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
1049	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
1050	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
1051	& -env_tke_max(ig,klcl(ig)))*interp(ig)
1052	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
1053	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
1054	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
1055	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
1056	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
1057	else
1058	IF (fraca0(ig) /= fraca0(ig) .OR. ABS(fraca0(ig)) > largest*10.e5) THEN
1059	PRINT *,' Lluis wrong fraca0(ig): ',fraca0(ig),' at : ',ig
1060	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1061	' fraca(ig,klcl(ig)): ',fraca(ig,klcl(ig)),' fraca(ig,klcl(ig)+1): ', &
1062	fraca(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1063	END IF
1064	IF (w_conv(ig) /= w_conv(ig) .OR. ABS(w_conv(ig)) > largest*10.e5) THEN
1065	PRINT *,' Lluis wrong w_conv(ig): ',w_conv(ig),' at : ',ig
1066	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1067	' w_ls(ig,klcl(ig)): ',w_ls(ig,klcl(ig)),' w_ls(ig,klcl(ig)+1): ', &
1068	w_ls(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1069	END IF
1070	fraca0(ig)=0.
1071	w0(ig)=0.
1072	! L. Fita, LMD July 2014. Adding zero value for stability issues
1073	w_conv(ig) = 0.
1074	!!jyg le 27/04/2012
1075	!! zlcl(ig)=0.
1076	!!
1077	endif
1078	enddo
1079
1080	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
1081	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
1082
1083	!------------Triggering------------------
1084	IF (iflag_trig_bl.ge.1) THEN
1085
1086	!-----Initialisations
1087	depth(:)=0.
1088	n2(:)=0.
1089	s2(:)=0.
1090	s_max(:)=0.
1091
1092	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
1093	do ig=1,ngrid
1094	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
1095	depth(ig)=zmax_moy(ig)-zlcl(ig)
1096	hmin(ig)=hmincoef*zlcl(ig)
1097	if (depth(ig).ge.10.) then
1098	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
1099	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
1100	!!
1101	!!jyg le 27/04/2012
1102	!! s_max(ig)=s2(ig)*log(n2(ig))
1103	!! if (n2(ig) .lt. 1) s_max(ig)=0.
1104	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
1105	!!fin jyg
1106	else
1107	s2(ig)=0.
1108	n2(ig)=0.
1109	s_max(ig)=0.
1110	endif
1111	enddo
1112	! print *,'avant Calcul de Wmax ' !!jyg
1113
1114	!-----Calcul de Wmax et ALE_BL_STAT associée
1115	!!jyg le 30/04/2012
1116	!! do ig=1,ngrid
1117	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
1118	!! 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))))
1119	!! ale_bl_stat(ig)=0.5w_max(ig)*2
1120	!! else
1121	!! w_max(ig)=0.
1122	!! ale_bl_stat(ig)=0.
1123	!! endif
1124	!! enddo
1125	susqr2pi=susqrt(2.Rpi)
1126	Reuler=exp(1.)
1127	do ig=1,ngrid
1128	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
1129	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
1130	ale_bl_stat(ig)=0.5w_max(ig)*2
1131	else
1132	w_max(ig)=0.
1133	ale_bl_stat(ig)=0.
1134	endif
1135	enddo
1136
1137	ENDIF ! iflag_trig_bl
1138	! print *,'ENDIF iflag_trig_bl' !!jyg
1139	lfname='thermcell_main before closure'
1140	lvarname = 'pt'
1141	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1142	lvarname = 'pdtadj'
1143	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1144	lvarname = 'rhobarz0'
1145	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1146	lvarname = 'fraca0'
1147	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1148	lvarname = 'w_conv'
1149	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1150	lvarname = 'interp'
1151	CALL check_var(lfname, lvarname, interp,ngrid, largest, .FALSE.)
1152	lvarname = 'w0'
1153	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1154	lvarname = 'therm_tke_max0'
1155	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1156	lvarname = 'env_tke_max0'
1157	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1158	lvarname = 'pbl_tke_max0'
1159	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1160	!------------Closure------------------
1161
1162	IF (iflag_clos_bl.ge.1) THEN
1163
1164	!-----Calcul de ALP_BL_STAT
1165	do ig=1,ngrid
1166	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
1167	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
1168	& (w0(ig)**2)
1169	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
1170	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
1171	if (iflag_clos_bl.ge.2) then
1172	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
1173	& (w0(ig)**2)
1174	else
1175	alp_bl_conv(ig)=0.
1176	endif
1177	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
1178	enddo
1179
1180	!-----Sécurité ALP infinie
1181	do ig=1,ngrid
1182	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
1183	enddo
1184
1185	ENDIF ! (iflag_clos_bl.ge.1)
1186	lfname='thermcell main after closure'
1187	lvarname = 'pt'
1188	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1189	lvarname = 'pdtadj'
1190	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1191	lvarname = 'alp_bl_det'
1192	CALL check_var(lfname, lvarname, alp_bl_det,ngrid, largest, .FALSE.)
1193	lvarname = 'rhobarz0'
1194	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1195	lvarname = 'fraca0'
1196	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1197	lvarname = 'w_conv'
1198	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1199	lvarname = 'w0'
1200	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1201	lvarname = 'alp_bl_fluct_m'
1202	CALL check_var(lfname, lvarname, alp_bl_fluct_m,ngrid, largest, .FALSE.)
1203	lvarname = 'alp_bl_fluct_tke'
1204	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,ngrid, largest, .FALSE.)
1205	lvarname = 'therm_tke_max0'
1206	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1207	lvarname = 'env_tke_max0'
1208	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1209	lvarname = 'pbl_tke_max0'
1210	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1211	lvarname = 'alp_bl_conv'
1212	CALL check_var(lfname, lvarname, alp_bl_conv,ngrid, largest, .FALSE.)
1213	lvarname = 'alp_bl_stat'
1214	CALL check_var(lfname, lvarname, alp_bl_stat,ngrid, largest, .FALSE.)
1215
1216	!!! fin nrlmd le 10/04/2012
1217
1218	if (prt_level.ge.10) then
1219	ig=igout
1220	do l=1,nlay
1221	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
1222	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
1223	enddo
1224	endif
1225
1226	! print*,'avant calcul ale et alp'
1227	!calcul de ALE et ALP pour la convection
1228	Alp_bl(:)=0.
1229	Ale_bl(:)=0.
1230	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
1231	do l=1,nlay
1232	do ig=1,ngrid
1233	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
1234	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
1235	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
1236	enddo
1237	enddo
1238
1239	!test:calcul de la ponderation des couches pour KE
1240	!initialisations
1241
1242	fm_tot(:)=0.
1243	wght_th(:,:)=1.
1244	lalim_conv(:)=lalim(:)
1245
1246	do k=1,nlay
1247	do ig=1,ngrid
1248	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
1249	enddo
1250	enddo
1251	lfname='after calculation of Al[p/e]_bl'
1252	lvarname = 'fm'
1253	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
1254	lvarname = 'rhobarz'
1255	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
1256	lvarname = 'wth3'
1257	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
1258	lvarname = 'Alp_bl'
1259	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1260	lvarname = 'Ale_bl'
1261	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1262
1263	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
1264	! plus petit que 1.e-10, on prend wght_th=1.
1265	do k=1,nlay
1266	do ig=1,ngrid
1267	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1268	wght_th(ig,k)=alim_star(ig,k)
1269	endif
1270	enddo
1271	enddo
1272
1273	! print*,'apres wght_th'
1274	!test pour prolonger la convection
1275	do ig=1,ngrid
1276	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1277	if ((alim_star(ig,1).lt.1.e-10)) then
1278	lalim_conv(ig)=1
1279	wght_th(ig,1)=1.
1280	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1281	endif
1282	enddo
1283
1284	!------------------------------------------------------------------------
1285	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1286	! Integrale verticale de ALP.
1287	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1288	! couches
1289	!------------------------------------------------------------------------
1290
1291	alp_int(:)=0.
1292	dp_int(:)=0.
1293	do l=2,nlay
1294	do ig=1,ngrid
1295	if(l.LE.lmax(ig)) THEN
1296	zdp=pplay(ig,l-1)-pplay(ig,l)
1297	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1298	dp_int(ig)=dp_int(ig)+zdp
1299	endif
1300	enddo
1301	enddo
1302
1303	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1304	do ig=1,ngrid
1305	!valeur integree de alp_bl * 0.5:
1306	if (dp_int(ig)>0.) then
1307	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1308	endif
1309	enddo!
1310	endif
1311
1312
1313	! Facteur multiplicatif sur Alp_bl
1314	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1315
1316	lfname='thermcell_main last computations on Alp_bl'
1317	lvarname = 'pplay'
1318	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest*10.e4, .FALSE.)
1319	lvarname = 'alp_int'
1320	CALL check_var(lfname, lvarname, alp_int, ngrid, largest*10.e4, .FALSE.)
1321	lvarname = 'dp_int'
1322	CALL check_var(lfname, lvarname, dp_int, ngrid, largest, .FALSE.)
1323	lvarname = 'Alp_bl'
1324	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1325
1326	!------------------------------------------------------------------------
1327
1328
1329	!calcul du ratqscdiff
1330	if (prt_level.ge.1) print*,'14e OK convect8'
1331	var=0.
1332	vardiff=0.
1333	ratqsdiff(:,:)=0.
1334
1335	do l=1,nlay
1336	do ig=1,ngrid
1337	if (l<=lalim(ig)) then
1338	var=var+alim_star(ig,l)zqta(ig,l)1000.
1339	endif
1340	enddo
1341	enddo
1342
1343	if (prt_level.ge.1) print*,'14f OK convect8'
1344
1345	do l=1,nlay
1346	do ig=1,ngrid
1347	if (l<=lalim(ig)) then
1348	zf=fraca(ig,l)
1349	zf2=zf/(1.-zf)
1350	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1351	endif
1352	enddo
1353	enddo
1354
1355	if (prt_level.ge.1) print*,'14g OK convect8'
1356	do l=1,nlay
1357	do ig=1,ngrid
1358	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1359	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1360	enddo
1361	enddo
1362	!--------------------------------------------------------------------
1363	!
1364	!ecriture des fichiers sortie
1365	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1366
1367	#ifdef wrgrads_thermcell
1368	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1369	#include "thermcell_out3d.h"
1370	#endif
1371
1372	endif
1373
1374	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1375
1376	lfname='before leaving thermcell_main'
1377	lvarname = 'pt'
1378	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1379	lvarname = 'pdtadj'
1380	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1381	lvarname = 'Alp_bl'
1382	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1383	lvarname = 'Ale_bl'
1384	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1385
1386
1387	return
1388	end
1389
1390	!-----------------------------------------------------------------------------
1391
1392	subroutine test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1393	IMPLICIT NONE
1394	#include "iniprint.h"
1395
1396	integer i, k, ngrid,nlay
1397	real pplev(ngrid,nlay+1),pplay(ngrid,nlay)
1398	real ztv(ngrid,nlay)
1399	real po(ngrid,nlay)
1400	real ztva(ngrid,nlay)
1401	real zqla(ngrid,nlay)
1402	real f_star(ngrid,nlay)
1403	real zw2(ngrid,nlay)
1404	integer long(ngrid)
1405	real seuil
1406	character*21 comment
1407
1408	if (prt_level.ge.1) THEN
1409	print*,'WARNING !!! TEST ',comment
1410	endif
1411	return
1412
1413	! test sur la hauteur des thermiques ...
1414	do i=1,ngrid
1415	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1416	if (prt_level.ge.10) then
1417	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1418	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1419	do k=1,nlay
1420	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)
1421	enddo
1422	endif
1423	enddo
1424
1425
1426	return
1427	end
1428
1429	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1430	! On transporte pbl_tke pour donner therm_tke
1431	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1432	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1433	& rg,pplev,therm_tke_max)
1434	implicit none
1435
1436	#include "iniprint.h"
1437	!=======================================================================
1438	!
1439	! Calcul du transport verticale dans la couche limite en presence
1440	! de "thermiques" explicitement representes
1441	! calcul du dq/dt une fois qu'on connait les ascendances
1442	!
1443	!=======================================================================
1444
1445	integer ngrid,nlay,nsrf
1446
1447	real ptimestep
1448	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1449	real entr0(ngrid,nlay),rg
1450	real therm_tke_max(ngrid,nlay)
1451	real detr0(ngrid,nlay)
1452
1453
1454	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1455	real entr(ngrid,nlay)
1456	real q(ngrid,nlay)
1457	integer lev_out ! niveau pour les print
1458
1459	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1460
1461	real zzm
1462
1463	integer ig,k
1464	integer isrf
1465
1466
1467	lev_out=0
1468
1469
1470	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1471
1472	! calcul du detrainement
1473	do k=1,nlay
1474	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1475	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1476	enddo
1477
1478
1479	! Decalage vertical des entrainements et detrainements.
1480	masse(:,1)=0.5*masse0(:,1)
1481	entr(:,1)=0.5*entr0(:,1)
1482	detr(:,1)=0.5*detr0(:,1)
1483	fm(:,1)=0.
1484	do k=1,nlay-1
1485	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1486	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1487	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1488	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1489	enddo
1490	fm(:,nlay+1)=0.
1491
1492	!!! nrlmd le 16/09/2010
1493	! calcul de la valeur dans les ascendances
1494	! do ig=1,ngrid
1495	! qa(ig,1)=q(ig,1)
1496	! enddo
1497	!!!
1498
1499	!do isrf=1,nsrf
1500
1501	! q(:,:)=therm_tke(:,:,isrf)
1502	q(:,:)=therm_tke_max(:,:)
1503	!!! nrlmd le 16/09/2010
1504	do ig=1,ngrid
1505	qa(ig,1)=q(ig,1)
1506	enddo
1507	!!!
1508
1509	if (1==1) then
1510	do k=2,nlay
1511	do ig=1,ngrid
1512	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1513	& 1.e-5*masse(ig,k)) then
1514	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1515	& /(fm(ig,k+1)+detr(ig,k))
1516	else
1517	qa(ig,k)=q(ig,k)
1518	endif
1519	if (qa(ig,k).lt.0.) then
1520	! print*,'qa<0!!!'
1521	endif
1522	if (q(ig,k).lt.0.) then
1523	! print*,'q<0!!!'
1524	endif
1525	enddo
1526	enddo
1527
1528	! Calcul du flux subsident
1529
1530	do k=2,nlay
1531	do ig=1,ngrid
1532	wqd(ig,k)=fm(ig,k)*q(ig,k)
1533	if (wqd(ig,k).lt.0.) then
1534	! print*,'wqd<0!!!'
1535	endif
1536	enddo
1537	enddo
1538	do ig=1,ngrid
1539	wqd(ig,1)=0.
1540	wqd(ig,nlay+1)=0.
1541	enddo
1542
1543	! Calcul des tendances
1544	do k=1,nlay
1545	do ig=1,ngrid
1546	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1547	& -wqd(ig,k)+wqd(ig,k+1)) &
1548	& *ptimestep/masse(ig,k)
1549	enddo
1550	enddo
1551
1552	endif
1553
1554	therm_tke_max(:,:)=q(:,:)
1555
1556	return
1557	!!! fin nrlmd le 10/04/2012
1558	end
1559

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

Download in other formats:

Original Format