Context Navigation

thermcell_main.F90 @ 125

Last change on this file since 125 was 125, checked in by lfita, 10 years ago
Adding the last initializations?
File size: 53.9 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	inerp = 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	klcl(ig)=l
918	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
919	endif
920	endif
921	enddo
922	enddo
923	lfname='thermcell_main calculation of LCL'
924	lvarname = 'inerp'
925	CALL check_var3D(lfname, lvarname, interp, ngrid, nlay, largest, .FALSE.)
926	lvarname = 'klcl'
927	CALL check_var(lfname, lvarname, REAL(klcl), ngrid, largest, .FALSE.)
928	lvarname = 'pcon'
929	CALL check_var(lfname, lvarname, pcon, ngrid, largest, .FALSE.)
930	lvarname = 'pplay'
931	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest, .FALSE.)
932
933	!------------Hauteur des thermiques
934	!!jyg le 27/04/2012
935	!! do ig =1,ngrid
936	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
937	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
938	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
939	!! zmax(ig)=pphi(ig,lmax(ig))/rg
940	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
941	!! enddo
942	do ig =1,ngrid
943	zmax(ig)=pphi(ig,lmax(ig))/rg
944	if (ok_lcl(ig)) then
945	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
946	& -rhobarz(ig,klcl(ig)))*interp(ig)
947	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
948	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
949	else
950	rhobarz0(ig)=0.
951	zlcl(ig)=zmax(ig)
952	endif
953	enddo
954	!!jyg fin
955	lfname='thermcell_main before LCL'
956	lvarname = 'rhobarz0'
957	CALL check_var(lfname, lvarname, rhobarz0, ngrid, largest, .FALSE.)
958	lvarname = 'pcon'
959	CALL check_var(lfname, lvarname, pcon, ngrid, largest, .FALSE.)
960	lvarname = 'fraca'
961	CALL check_var3D(lfname, lvarname, fraca, ngrid, nlay+1, largest, .FALSE.)
962	lvarname = 'fraca0'
963	CALL check_var(lfname, lvarname, fraca0, ngrid, largest, .FALSE.)
964	lvarname = 'w_conv'
965	CALL check_var(lfname, lvarname, w_conv, ngrid, largest, .FALSE.)
966	lvarname = 'interp'
967	CALL check_var(lfname, lvarname, interp, ngrid, largest, .FALSE.)
968	lvarname = 'zlcl'
969	CALL check_var(lfname, lvarname, zlcl, ngrid, largest, .FALSE.)
970
971	!------------Calcul des propriétés du thermique au LCL
972	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
973
974	!-----Initialisation de la TKE moyenne
975	do l=1,nlay
976	do ig=1,ngrid
977	pbl_tke_max(ig,l)=0.
978	enddo
979	enddo
980
981	!-----Calcul de la TKE moyenne
982	do nsrf=1,nbsrf
983	do l=1,nlay
984	do ig=1,ngrid
985	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
986	enddo
987	enddo
988	enddo
989
990	!-----Initialisations des TKE dans et hors des thermiques
991	do l=1,nlay
992	do ig=1,ngrid
993	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
994	env_tke_max(ig,l)=pbl_tke_max(ig,l)
995	enddo
996	enddo
997
998	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
999	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1000	& rg,pplev,therm_tke_max)
1001	! print *,' thermcell_tke_transport -> ' !!jyg
1002
1003	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
1004	do l=1,nlay
1005	do ig=1,ngrid
1006	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
1007	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
1008	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
1009	enddo
1010	enddo
1011	! print *,' apres w_ls = ' !!jyg
1012
1013	do ig=1,ngrid
1014	if (ok_lcl(ig)) then
1015	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
1016	& -fraca(ig,klcl(ig)))*interp(ig)
1017	IF (fraca0(ig) /= fraca0(ig) .OR. ABS(fraca0(ig)) > largest*10.e5) THEN
1018	PRINT *,' Lluis wrong fraca0(ig): ',fraca0(ig),' at : ',ig
1019	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1020	' fraca(ig,klcl(ig)): ',fraca(ig,klcl(ig)),' fraca(ig,klcl(ig)+1): ', &
1021	fraca(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1022	END IF
1023
1024	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
1025	& -zw2(ig,klcl(ig)))*interp(ig)
1026	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
1027	& -w_ls(ig,klcl(ig)))*interp(ig)
1028	IF (w_conv(ig) /= w_conv(ig) .OR. ABS(w_conv(ig)) > largest*10.e5) THEN
1029	PRINT *,' Lluis wrong w_conv(ig): ',w_conv(ig),' at : ',ig
1030	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1031	' w_ls(ig,klcl(ig)): ',w_ls(ig,klcl(ig)),' w_ls(ig,klcl(ig)+1): ', &
1032	w_ls(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1033	END IF
1034	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
1035	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
1036	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
1037	& -env_tke_max(ig,klcl(ig)))*interp(ig)
1038	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
1039	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
1040	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
1041	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
1042	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
1043	else
1044	IF (fraca0(ig) /= fraca0(ig) .OR. ABS(fraca0(ig)) > largest*10.e5) THEN
1045	PRINT *,' Lluis wrong fraca0(ig): ',fraca0(ig),' at : ',ig
1046	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1047	' fraca(ig,klcl(ig)): ',fraca(ig,klcl(ig)),' fraca(ig,klcl(ig)+1): ', &
1048	fraca(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1049	END IF
1050	IF (w_conv(ig) /= w_conv(ig) .OR. ABS(w_conv(ig)) > largest*10.e5) THEN
1051	PRINT *,' Lluis wrong w_conv(ig): ',w_conv(ig),' at : ',ig
1052	PRINT *,' klcl(ig): ', klcl(ig),' klcl(ig)+1: ',klcl(ig)+1, &
1053	' w_ls(ig,klcl(ig)): ',w_ls(ig,klcl(ig)),' w_ls(ig,klcl(ig)+1): ', &
1054	w_ls(ig,klcl(ig)+1), ' interp(ig): ',interp(ig)
1055	END IF
1056	fraca0(ig)=0.
1057	w0(ig)=0.
1058	! L. Fita, LMD July 2014. Adding zero value for stability issues
1059	w_conv(ig) = 0.
1060	!!jyg le 27/04/2012
1061	!! zlcl(ig)=0.
1062	!!
1063	endif
1064	enddo
1065
1066	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
1067	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
1068
1069	!------------Triggering------------------
1070	IF (iflag_trig_bl.ge.1) THEN
1071
1072	!-----Initialisations
1073	depth(:)=0.
1074	n2(:)=0.
1075	s2(:)=0.
1076	s_max(:)=0.
1077
1078	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
1079	do ig=1,ngrid
1080	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
1081	depth(ig)=zmax_moy(ig)-zlcl(ig)
1082	hmin(ig)=hmincoef*zlcl(ig)
1083	if (depth(ig).ge.10.) then
1084	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
1085	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
1086	!!
1087	!!jyg le 27/04/2012
1088	!! s_max(ig)=s2(ig)*log(n2(ig))
1089	!! if (n2(ig) .lt. 1) s_max(ig)=0.
1090	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
1091	!!fin jyg
1092	else
1093	s2(ig)=0.
1094	n2(ig)=0.
1095	s_max(ig)=0.
1096	endif
1097	enddo
1098	! print *,'avant Calcul de Wmax ' !!jyg
1099
1100	!-----Calcul de Wmax et ALE_BL_STAT associée
1101	!!jyg le 30/04/2012
1102	!! do ig=1,ngrid
1103	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
1104	!! 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))))
1105	!! ale_bl_stat(ig)=0.5w_max(ig)*2
1106	!! else
1107	!! w_max(ig)=0.
1108	!! ale_bl_stat(ig)=0.
1109	!! endif
1110	!! enddo
1111	susqr2pi=susqrt(2.Rpi)
1112	Reuler=exp(1.)
1113	do ig=1,ngrid
1114	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
1115	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
1116	ale_bl_stat(ig)=0.5w_max(ig)*2
1117	else
1118	w_max(ig)=0.
1119	ale_bl_stat(ig)=0.
1120	endif
1121	enddo
1122
1123	ENDIF ! iflag_trig_bl
1124	! print *,'ENDIF iflag_trig_bl' !!jyg
1125	lfname='thermcell_main before closure'
1126	lvarname = 'pt'
1127	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1128	lvarname = 'pdtadj'
1129	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1130	lvarname = 'rhobarz0'
1131	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1132	lvarname = 'fraca0'
1133	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1134	lvarname = 'w_conv'
1135	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1136	lvarname = 'interp'
1137	CALL check_var(lfname, lvarname, interp,ngrid, largest, .FALSE.)
1138	lvarname = 'w0'
1139	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1140	lvarname = 'therm_tke_max0'
1141	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1142	lvarname = 'env_tke_max0'
1143	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1144	lvarname = 'pbl_tke_max0'
1145	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1146	!------------Closure------------------
1147
1148	IF (iflag_clos_bl.ge.1) THEN
1149
1150	!-----Calcul de ALP_BL_STAT
1151	do ig=1,ngrid
1152	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
1153	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
1154	& (w0(ig)**2)
1155	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
1156	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
1157	if (iflag_clos_bl.ge.2) then
1158	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
1159	& (w0(ig)**2)
1160	else
1161	alp_bl_conv(ig)=0.
1162	endif
1163	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
1164	enddo
1165
1166	!-----Sécurité ALP infinie
1167	do ig=1,ngrid
1168	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
1169	enddo
1170
1171	ENDIF ! (iflag_clos_bl.ge.1)
1172	lfname='thermcell main after closure'
1173	lvarname = 'pt'
1174	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1175	lvarname = 'pdtadj'
1176	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1177	lvarname = 'alp_bl_det'
1178	CALL check_var(lfname, lvarname, alp_bl_det,ngrid, largest, .FALSE.)
1179	lvarname = 'rhobarz0'
1180	CALL check_var(lfname, lvarname, rhobarz0,ngrid, largest, .FALSE.)
1181	lvarname = 'fraca0'
1182	CALL check_var(lfname, lvarname, fraca0,ngrid, largest, .FALSE.)
1183	lvarname = 'w_conv'
1184	CALL check_var(lfname, lvarname, w_conv,ngrid, largest, .FALSE.)
1185	lvarname = 'w0'
1186	CALL check_var(lfname, lvarname, w0,ngrid, largest, .FALSE.)
1187	lvarname = 'alp_bl_fluct_m'
1188	CALL check_var(lfname, lvarname, alp_bl_fluct_m,ngrid, largest, .FALSE.)
1189	lvarname = 'alp_bl_fluct_tke'
1190	CALL check_var(lfname, lvarname, alp_bl_fluct_tke,ngrid, largest, .FALSE.)
1191	lvarname = 'therm_tke_max0'
1192	CALL check_var(lfname, lvarname, therm_tke_max0,ngrid, largest, .FALSE.)
1193	lvarname = 'env_tke_max0'
1194	CALL check_var(lfname, lvarname, env_tke_max0,ngrid, largest, .FALSE.)
1195	lvarname = 'pbl_tke_max0'
1196	CALL check_var(lfname, lvarname, pbl_tke_max0,ngrid, largest, .FALSE.)
1197	lvarname = 'alp_bl_conv'
1198	CALL check_var(lfname, lvarname, alp_bl_conv,ngrid, largest, .FALSE.)
1199	lvarname = 'alp_bl_stat'
1200	CALL check_var(lfname, lvarname, alp_bl_stat,ngrid, largest, .FALSE.)
1201
1202	!!! fin nrlmd le 10/04/2012
1203
1204	if (prt_level.ge.10) then
1205	ig=igout
1206	do l=1,nlay
1207	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
1208	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
1209	enddo
1210	endif
1211
1212	! print*,'avant calcul ale et alp'
1213	!calcul de ALE et ALP pour la convection
1214	Alp_bl(:)=0.
1215	Ale_bl(:)=0.
1216	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
1217	do l=1,nlay
1218	do ig=1,ngrid
1219	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
1220	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
1221	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
1222	enddo
1223	enddo
1224
1225	!test:calcul de la ponderation des couches pour KE
1226	!initialisations
1227
1228	fm_tot(:)=0.
1229	wght_th(:,:)=1.
1230	lalim_conv(:)=lalim(:)
1231
1232	do k=1,nlay
1233	do ig=1,ngrid
1234	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
1235	enddo
1236	enddo
1237	lfname='after calculation of Al[p/e]_bl'
1238	lvarname = 'fm'
1239	CALL check_var3D(lfname, lvarname, fm, ngrid, nlay, largest, .FALSE.)
1240	lvarname = 'rhobarz'
1241	CALL check_var3D(lfname, lvarname, rhobarz, ngrid, nlay, largest, .FALSE.)
1242	lvarname = 'wth3'
1243	CALL check_var3D(lfname, lvarname, wth3, ngrid, nlay, largest, .FALSE.)
1244	lvarname = 'Alp_bl'
1245	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1246	lvarname = 'Ale_bl'
1247	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1248
1249	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
1250	! plus petit que 1.e-10, on prend wght_th=1.
1251	do k=1,nlay
1252	do ig=1,ngrid
1253	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
1254	wght_th(ig,k)=alim_star(ig,k)
1255	endif
1256	enddo
1257	enddo
1258
1259	! print*,'apres wght_th'
1260	!test pour prolonger la convection
1261	do ig=1,ngrid
1262	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
1263	if ((alim_star(ig,1).lt.1.e-10)) then
1264	lalim_conv(ig)=1
1265	wght_th(ig,1)=1.
1266	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
1267	endif
1268	enddo
1269
1270	!------------------------------------------------------------------------
1271	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1272	! Integrale verticale de ALP.
1273	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1274	! couches
1275	!------------------------------------------------------------------------
1276
1277	alp_int(:)=0.
1278	dp_int(:)=0.
1279	do l=2,nlay
1280	do ig=1,ngrid
1281	if(l.LE.lmax(ig)) THEN
1282	zdp=pplay(ig,l-1)-pplay(ig,l)
1283	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1284	dp_int(ig)=dp_int(ig)+zdp
1285	endif
1286	enddo
1287	enddo
1288
1289	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1290	do ig=1,ngrid
1291	!valeur integree de alp_bl * 0.5:
1292	if (dp_int(ig)>0.) then
1293	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1294	endif
1295	enddo!
1296	endif
1297
1298
1299	! Facteur multiplicatif sur Alp_bl
1300	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1301
1302	lfname='thermcell_main last computations on Alp_bl'
1303	lvarname = 'pplay'
1304	CALL check_var3D(lfname, lvarname, pplay, ngrid, nlay, largest*10.e4, .FALSE.)
1305	lvarname = 'alp_int'
1306	CALL check_var(lfname, lvarname, alp_int, ngrid, largest*10.e4, .FALSE.)
1307	lvarname = 'dp_int'
1308	CALL check_var(lfname, lvarname, dp_int, ngrid, largest, .FALSE.)
1309	lvarname = 'Alp_bl'
1310	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1311
1312	!------------------------------------------------------------------------
1313
1314
1315	!calcul du ratqscdiff
1316	if (prt_level.ge.1) print*,'14e OK convect8'
1317	var=0.
1318	vardiff=0.
1319	ratqsdiff(:,:)=0.
1320
1321	do l=1,nlay
1322	do ig=1,ngrid
1323	if (l<=lalim(ig)) then
1324	var=var+alim_star(ig,l)zqta(ig,l)1000.
1325	endif
1326	enddo
1327	enddo
1328
1329	if (prt_level.ge.1) print*,'14f OK convect8'
1330
1331	do l=1,nlay
1332	do ig=1,ngrid
1333	if (l<=lalim(ig)) then
1334	zf=fraca(ig,l)
1335	zf2=zf/(1.-zf)
1336	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1337	endif
1338	enddo
1339	enddo
1340
1341	if (prt_level.ge.1) print*,'14g OK convect8'
1342	do l=1,nlay
1343	do ig=1,ngrid
1344	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1345	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1346	enddo
1347	enddo
1348	!--------------------------------------------------------------------
1349	!
1350	!ecriture des fichiers sortie
1351	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1352
1353	#ifdef wrgrads_thermcell
1354	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1355	#include "thermcell_out3d.h"
1356	#endif
1357
1358	endif
1359
1360	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1361
1362	lfname='before leaving thermcell_main'
1363	lvarname = 'pt'
1364	CALL check_var3D(lfname, lvarname, pt, ngrid, nlay, largest, .FALSE.)
1365	lvarname = 'pdtadj'
1366	CALL check_var3D(lfname, lvarname, pdtadj, ngrid, nlay, largest, .FALSE.)
1367	lvarname = 'Alp_bl'
1368	CALL check_var(lfname, lvarname, Alp_bl, ngrid, largest, .FALSE.)
1369	lvarname = 'Ale_bl'
1370	CALL check_var(lfname, lvarname, Ale_bl, ngrid, largest, .FALSE.)
1371
1372
1373	return
1374	end
1375
1376	!-----------------------------------------------------------------------------
1377
1378	subroutine test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1379	IMPLICIT NONE
1380	#include "iniprint.h"
1381
1382	integer i, k, ngrid,nlay
1383	real pplev(ngrid,nlay+1),pplay(ngrid,nlay)
1384	real ztv(ngrid,nlay)
1385	real po(ngrid,nlay)
1386	real ztva(ngrid,nlay)
1387	real zqla(ngrid,nlay)
1388	real f_star(ngrid,nlay)
1389	real zw2(ngrid,nlay)
1390	integer long(ngrid)
1391	real seuil
1392	character*21 comment
1393
1394	if (prt_level.ge.1) THEN
1395	print*,'WARNING !!! TEST ',comment
1396	endif
1397	return
1398
1399	! test sur la hauteur des thermiques ...
1400	do i=1,ngrid
1401	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1402	if (prt_level.ge.10) then
1403	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1404	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1405	do k=1,nlay
1406	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)
1407	enddo
1408	endif
1409	enddo
1410
1411
1412	return
1413	end
1414
1415	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1416	! On transporte pbl_tke pour donner therm_tke
1417	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1418	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1419	& rg,pplev,therm_tke_max)
1420	implicit none
1421
1422	#include "iniprint.h"
1423	!=======================================================================
1424	!
1425	! Calcul du transport verticale dans la couche limite en presence
1426	! de "thermiques" explicitement representes
1427	! calcul du dq/dt une fois qu'on connait les ascendances
1428	!
1429	!=======================================================================
1430
1431	integer ngrid,nlay,nsrf
1432
1433	real ptimestep
1434	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1435	real entr0(ngrid,nlay),rg
1436	real therm_tke_max(ngrid,nlay)
1437	real detr0(ngrid,nlay)
1438
1439
1440	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1441	real entr(ngrid,nlay)
1442	real q(ngrid,nlay)
1443	integer lev_out ! niveau pour les print
1444
1445	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1446
1447	real zzm
1448
1449	integer ig,k
1450	integer isrf
1451
1452
1453	lev_out=0
1454
1455
1456	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1457
1458	! calcul du detrainement
1459	do k=1,nlay
1460	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1461	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1462	enddo
1463
1464
1465	! Decalage vertical des entrainements et detrainements.
1466	masse(:,1)=0.5*masse0(:,1)
1467	entr(:,1)=0.5*entr0(:,1)
1468	detr(:,1)=0.5*detr0(:,1)
1469	fm(:,1)=0.
1470	do k=1,nlay-1
1471	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1472	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1473	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1474	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1475	enddo
1476	fm(:,nlay+1)=0.
1477
1478	!!! nrlmd le 16/09/2010
1479	! calcul de la valeur dans les ascendances
1480	! do ig=1,ngrid
1481	! qa(ig,1)=q(ig,1)
1482	! enddo
1483	!!!
1484
1485	!do isrf=1,nsrf
1486
1487	! q(:,:)=therm_tke(:,:,isrf)
1488	q(:,:)=therm_tke_max(:,:)
1489	!!! nrlmd le 16/09/2010
1490	do ig=1,ngrid
1491	qa(ig,1)=q(ig,1)
1492	enddo
1493	!!!
1494
1495	if (1==1) then
1496	do k=2,nlay
1497	do ig=1,ngrid
1498	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1499	& 1.e-5*masse(ig,k)) then
1500	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1501	& /(fm(ig,k+1)+detr(ig,k))
1502	else
1503	qa(ig,k)=q(ig,k)
1504	endif
1505	if (qa(ig,k).lt.0.) then
1506	! print*,'qa<0!!!'
1507	endif
1508	if (q(ig,k).lt.0.) then
1509	! print*,'q<0!!!'
1510	endif
1511	enddo
1512	enddo
1513
1514	! Calcul du flux subsident
1515
1516	do k=2,nlay
1517	do ig=1,ngrid
1518	wqd(ig,k)=fm(ig,k)*q(ig,k)
1519	if (wqd(ig,k).lt.0.) then
1520	! print*,'wqd<0!!!'
1521	endif
1522	enddo
1523	enddo
1524	do ig=1,ngrid
1525	wqd(ig,1)=0.
1526	wqd(ig,nlay+1)=0.
1527	enddo
1528
1529	! Calcul des tendances
1530	do k=1,nlay
1531	do ig=1,ngrid
1532	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1533	& -wqd(ig,k)+wqd(ig,k+1)) &
1534	& *ptimestep/masse(ig,k)
1535	enddo
1536	enddo
1537
1538	endif
1539
1540	therm_tke_max(:,:)=q(:,:)
1541
1542	return
1543	!!! fin nrlmd le 10/04/2012
1544	end
1545

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

Download in other formats:

Original Format