Context Navigation

thermcell_old.F @ 979

Last change on this file since 979 was 940, checked in by Laurent Fairhead, 17 years ago
On remplace le fichier include dimphy.h par le module dimphy.F90i pour etre coherent avec le partout LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 188.8 KB

Line
1	SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep
2	s ,pplay,pplev,pphi
3	s ,pu,pv,pt,po
4	s ,pduadj,pdvadj,pdtadj,pdoadj
5	s ,fm0,entr0
6	c s ,pu_therm,pv_therm
7	s ,r_aspect,l_mix,w2di,tho)
8
9	USE dimphy
10	IMPLICIT NONE
11
12	c=======================================================================
13	c
14	c Calcul du transport verticale dans la couche limite en presence
15	c de "thermiques" explicitement representes
16	c
17	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
18	c
19	c le thermique est supposé homogène et dissipé par mélange avec
20	c son environnement. la longueur l_mix contrôle l'efficacité du
21	c mélange
22	c
23	c Le calcul du transport des différentes espèces se fait en prenant
24	c en compte:
25	c 1. un flux de masse montant
26	c 2. un flux de masse descendant
27	c 3. un entrainement
28	c 4. un detrainement
29	c
30	c=======================================================================
31
32	c-----------------------------------------------------------------------
33	c declarations:
34	c -------------
35
36	#include "dimensions.h"
37	cccc#include "dimphy.h"
38	#include "YOMCST.h"
39
40	c arguments:
41	c ----------
42
43	INTEGER ngrid,nlay,w2di,tho
44	real ptimestep,l_mix,r_aspect
45	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
46	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
47	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
48	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
49	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
50	real pphi(ngrid,nlay)
51
52	integer idetr
53	save idetr
54	data idetr/3/
55
56	c local:
57	c ------
58
59	INTEGER ig,k,l,lmax(klon,klev+1),lmaxa(klon),lmix(klon)
60	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
61
62	real zlev(klon,klev+1),zlay(klon,klev)
63	REAL zh(klon,klev),zdhadj(klon,klev)
64	REAL ztv(klon,klev)
65	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
66	REAL wh(klon,klev+1)
67	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
68	real zla(klon,klev+1)
69	real zwa(klon,klev+1)
70	real zld(klon,klev+1)
71	real zwd(klon,klev+1)
72	real zsortie(klon,klev)
73	real zva(klon,klev)
74	real zua(klon,klev)
75	real zoa(klon,klev)
76
77	real zha(klon,klev)
78	real wa_moy(klon,klev+1)
79	real fraca(klon,klev+1)
80	real fracc(klon,klev+1)
81	real zf,zf2
82	real thetath2(klon,klev),wth2(klon,klev)
83	! common/comtherm/thetath2,wth2
84
85	real count_time
86	integer isplit,nsplit,ialt
87	parameter (nsplit=10)
88	data isplit/0/
89	save isplit
90
91	logical sorties
92	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
93	real zpspsk(klon,klev)
94
95	real wmax(klon,klev),wmaxa(klon)
96
97	real wa(klon,klev,klev+1)
98	real wd(klon,klev+1)
99	real larg_part(klon,klev,klev+1)
100	real fracd(klon,klev+1)
101	real xxx(klon,klev+1)
102	real larg_cons(klon,klev+1)
103	real larg_detr(klon,klev+1)
104	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
105	real pu_therm(klon,klev),pv_therm(klon,klev)
106	real fm(klon,klev+1),entr(klon,klev)
107	real fmc(klon,klev+1)
108
109	character (len=2) :: str2
110	character (len=10) :: str10
111
112	LOGICAL vtest(klon),down
113
114	EXTERNAL SCOPY
115
116	integer ncorrec,ll
117	save ncorrec
118	data ncorrec/0/
119	c
120	c-----------------------------------------------------------------------
121	c initialisation:
122	c ---------------
123	c
124	sorties=.true.
125	IF(ngrid.NE.klon) THEN
126	PRINT*
127	PRINT*,'STOP dans convadj'
128	PRINT*,'ngrid =',ngrid
129	PRINT*,'klon =',klon
130	ENDIF
131	c
132	c-----------------------------------------------------------------------
133	c incrementation eventuelle de tendances precedentes:
134	c ---------------------------------------------------
135
136	print*,'0 OK convect8'
137
138	DO 1010 l=1,nlay
139	DO 1015 ig=1,ngrid
140	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
141	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
142	zu(ig,l)=pu(ig,l)
143	zv(ig,l)=pv(ig,l)
144	zo(ig,l)=po(ig,l)
145	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
146	1015 CONTINUE
147	1010 CONTINUE
148
149	c print*,'1 OK convect8'
150	c --------------------
151	c
152	c
153	c + + + + + + + + + + +
154	c
155	c
156	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
157	c wh,wt,wo ...
158	c
159	c + + + + + + + + + + + zh,zu,zv,zo,rho
160	c
161	c
162	c -------------------- zlev(1)
163	c \\\\\\\\\\\\\\\\\\\\
164	c
165	c
166
167	c-----------------------------------------------------------------------
168	c Calcul des altitudes des couches
169	c-----------------------------------------------------------------------
170
171	do l=2,nlay
172	do ig=1,ngrid
173	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
174	enddo
175	enddo
176	do ig=1,ngrid
177	zlev(ig,1)=0.
178	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
179	enddo
180	do l=1,nlay
181	do ig=1,ngrid
182	zlay(ig,l)=pphi(ig,l)/RG
183	enddo
184	enddo
185
186	c print*,'2 OK convect8'
187	c-----------------------------------------------------------------------
188	c Calcul des densites
189	c-----------------------------------------------------------------------
190
191	do l=1,nlay
192	do ig=1,ngrid
193	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
194	enddo
195	enddo
196
197	do l=2,nlay
198	do ig=1,ngrid
199	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
200	enddo
201	enddo
202
203	do k=1,nlay
204	do l=1,nlay+1
205	do ig=1,ngrid
206	wa(ig,k,l)=0.
207	enddo
208	enddo
209	enddo
210
211	c print*,'3 OK convect8'
212	c------------------------------------------------------------------
213	c Calcul de w2, quarre de w a partir de la cape
214	c a partir de w2, on calcule wa, vitesse de l'ascendance
215	c
216	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
217	c w2 est stoke dans wa
218	c
219	c ATTENTION: dans convect8, on n'utilise le calcule des wa
220	c independants par couches que pour calculer l'entrainement
221	c a la base et la hauteur max de l'ascendance.
222	c
223	c Indicages:
224	c l'ascendance provenant du niveau k traverse l'interface l avec
225	c une vitesse wa(k,l).
226	c
227	c --------------------
228	c
229	c + + + + + + + + + +
230	c
231	c wa(k,l) ---- -------------------- l
232	c /\
233	c /\|\|\ + + + + + + + + + +
234	c \|\|
235	c \|\| --------------------
236	c \|\|
237	c \|\| + + + + + + + + + +
238	c \|\|
239	c \|\| --------------------
240	c \|\|__
241	c \|___ + + + + + + + + + + k
242	c
243	c --------------------
244	c
245	c
246	c
247	c------------------------------------------------------------------
248
249
250	do k=1,nlay-1
251	do ig=1,ngrid
252	wa(ig,k,k)=0.
253	wa(ig,k,k+1)=2.RG(ztv(ig,k)-ztv(ig,k+1))/ztv(ig,k+1)
254	s *(zlev(ig,k+1)-zlev(ig,k))
255	enddo
256	do l=k+1,nlay-1
257	do ig=1,ngrid
258	wa(ig,k,l+1)=wa(ig,k,l)+
259	s 2.RG(ztv(ig,k)-ztv(ig,l))/ztv(ig,l)
260	s *(zlev(ig,l+1)-zlev(ig,l))
261	enddo
262	enddo
263	do ig=1,ngrid
264	wa(ig,k,nlay+1)=0.
265	enddo
266	enddo
267
268	c print*,'4 OK convect8'
269	c Calcul de la couche correspondant a la hauteur du thermique
270	do k=1,nlay-1
271	do ig=1,ngrid
272	lmax(ig,k)=k
273	enddo
274	do l=nlay,k+1,-1
275	do ig=1,ngrid
276	if(wa(ig,k,l).le.1.e-10) lmax(ig,k)=l-1
277	enddo
278	enddo
279	enddo
280
281	c print*,'5 OK convect8'
282	c Calcule du w max du thermique
283	do k=1,nlay
284	do ig=1,ngrid
285	wmax(ig,k)=0.
286	enddo
287	enddo
288
289	do k=1,nlay-1
290	do l=k,nlay
291	do ig=1,ngrid
292	if (l.le.lmax(ig,k)) then
293	wa(ig,k,l)=sqrt(wa(ig,k,l))
294	wmax(ig,k)=max(wmax(ig,k),wa(ig,k,l))
295	else
296	wa(ig,k,l)=0.
297	endif
298	enddo
299	enddo
300	enddo
301
302	do k=1,nlay-1
303	do ig=1,ngrid
304	pu_therm(ig,k)=sqrt(wmax(ig,k))
305	pv_therm(ig,k)=sqrt(wmax(ig,k))
306	enddo
307	enddo
308
309	c print*,'6 OK convect8'
310	c Longueur caracteristique correspondant a la hauteur des thermiques.
311	do ig=1,ngrid
312	zmax(ig)=500.
313	enddo
314	c print*,'LMAX LMAX LMAX '
315	do k=1,nlay-1
316	do ig=1,ngrid
317	zmax(ig)=max(zmax(ig),zlev(ig,lmax(ig,k))-zlev(ig,k))
318	enddo
319	c print*,k,lmax(1,k)
320	enddo
321	c print*,'ZMAX ZMAX ZMAX ',zmax
322	c call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX ')
323
324	c Calcul de l'entrainement.
325	c Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
326	c de la couche d'alimentation en partant du principe que la vitesse
327	c maximum dans l'ascendance est la vitesse d'entrainement horizontale.
328	do k=1,nlay
329	do ig=1,ngrid
330	zzz=rho(ig,k)wmax(ig,k)(zlev(ig,k+1)-zlev(ig,k))
331	s /(zmax(ig)*r_aspect)
332	if(w2di.eq.2) then
333	entr(ig,k)=entr(ig,k)+
334	s ptimestep*(zzz-entr(ig,k))/float(tho)
335	else
336	entr(ig,k)=zzz
337	endif
338	ztva(ig,k)=ztv(ig,k)
339	enddo
340	enddo
341
342	c print*,'7 OK convect8'
343	do k=1,klev+1
344	do ig=1,ngrid
345	zw2(ig,k)=0.
346	fmc(ig,k)=0.
347	larg_cons(ig,k)=0.
348	larg_detr(ig,k)=0.
349	wa_moy(ig,k)=0.
350	enddo
351	enddo
352
353	c print*,'8 OK convect8'
354	do ig=1,ngrid
355	lmaxa(ig)=1
356	lmix(ig)=1
357	wmaxa(ig)=0.
358	enddo
359
360
361	do l=1,nlay-2
362	do ig=1,ngrid
363	c if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)) then
364	c print*,'COUCOU ',l,zw2(ig,l),ztv(ig,l),ztv(ig,l+1)
365	if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)
366	s .and.entr(ig,l).gt.1.e-10) then
367	c print*,'COUCOU cas 1'
368	c Initialisation de l'ascendance
369	c lmix(ig)=1
370	ztva(ig,l)=ztv(ig,l)
371	fmc(ig,l)=0.
372	fmc(ig,l+1)=entr(ig,l)
373	zw2(ig,l)=0.
374	c if (.not.ztv(ig,l+1).gt.150.) then
375	c print*,'ig,l+1,ztv(ig,l+1)'
376	c print*, ig,l+1,ztv(ig,l+1)
377	c stop'dans thermiques'
378	c endif
379	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
380	s *(zlev(ig,l+1)-zlev(ig,l))
381	larg_detr(ig,l)=0.
382	else if (zw2(ig,l).ge.1.e-10.and.
383	. fmc(ig,l)+entr(ig,l).gt.1.e-10) then
384	c Incrementation...
385	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
386	c if (.not.fmc(ig,l+1).gt.1.e-15) then
387	c print*,'ig,l+1,fmc(ig,l+1)'
388	c print*, ig,l+1,fmc(ig,l+1)
389	c print*,'Fmc ',(fmc(ig,ll),ll=1,klev+1)
390	c print*,'W2 ',(zw2(ig,ll),ll=1,klev+1)
391	c print*,'Tv ',(ztv(ig,ll),ll=1,klev)
392	c print*,'Entr ',(entr(ig,ll),ll=1,klev)
393	c stop'dans thermiques'
394	c endif
395	ztva(ig,l)=(fmc(ig,l)ztva(ig,l-1)+entr(ig,l)ztv(ig,l))
396	s /fmc(ig,l+1)
397	c mise a jour de la vitesse ascendante (l'air entraine de la couche
398	c consideree commence avec une vitesse nulle).
399	zw2(ig,l+1)=zw2(ig,l)(fmc(ig,l)/fmc(ig,l+1))*2+
400	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
401	s *(zlev(ig,l+1)-zlev(ig,l))
402	endif
403	if (zw2(ig,l+1).lt.0.) then
404	zw2(ig,l+1)=0.
405	lmaxa(ig)=l
406	else
407	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
408	endif
409	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
410	c lmix est le niveau de la couche ou w (wa_moy) est maximum
411	lmix(ig)=l+1
412	wmaxa(ig)=wa_moy(ig,l+1)
413	endif
414	c print*,'COUCOU cas 2 LMIX=',lmix(ig),wa_moy(ig,l+1),wmaxa(ig)
415	enddo
416	enddo
417
418	c print*,'9 OK convect8'
419	c print*,'WA1 ',wa_moy
420
421	c determination de l'indice du debut de la mixed layer ou w decroit
422
423	c calcul de la largeur de chaque ascendance dans le cas conservatif.
424	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
425	c d'une couche est égale à la hauteur de la couche alimentante.
426	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
427	c de la vitesse d'entrainement horizontal dans la couche alimentante.
428
429	do l=2,nlay
430	do ig=1,ngrid
431	if (l.le.lmaxa(ig)) then
432	zw=max(wa_moy(ig,l),1.e-10)
433	larg_cons(ig,l)=zmax(ig)*r_aspect
434	s fmc(ig,l)/(rhobarz(ig,l)zw)
435	endif
436	enddo
437	enddo
438
439	do l=2,nlay
440	do ig=1,ngrid
441	if (l.le.lmaxa(ig)) then
442	c if (idetr.eq.0) then
443	c cette option est finalement en dur.
444	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
445	c else if (idetr.eq.1) then
446	c larg_detr(ig,l)=larg_cons(ig,l)
447	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
448	c else if (idetr.eq.2) then
449	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
450	c s *sqrt(wa_moy(ig,l))
451	c else if (idetr.eq.4) then
452	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
453	c s *wa_moy(ig,l)
454	c endif
455	endif
456	enddo
457	enddo
458
459	c print*,'10 OK convect8'
460	c print*,'WA2 ',wa_moy
461	c calcul de la fraction de la maille concernée par l'ascendance en tenant
462	c compte de l'epluchage du thermique.
463
464	do l=2,nlay
465	do ig=1,ngrid
466	if(larg_cons(ig,l).gt.1.) then
467	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
468	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
469	s /(r_aspect*zmax(ig))
470	if(l.gt.lmix(ig)) then
471	xxx(ig,l)=(lmaxa(ig)+1.-l) / (lmaxa(ig)+1.-lmix(ig))
472	if (idetr.eq.0) then
473	fraca(ig,l)=fraca(ig,lmix(ig))
474	else if (idetr.eq.1) then
475	fraca(ig,l)=fraca(ig,lmix(ig))*xxx(ig,l)
476	else if (idetr.eq.2) then
477	fraca(ig,l)=fraca(ig,lmix(ig))(1.-(1.-xxx(ig,l))*2)
478	else
479	fraca(ig,l)=fraca(ig,lmix(ig))xxx(ig,l)*2
480	endif
481	endif
482	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
483	fraca(ig,l)=max(fraca(ig,l),0.)
484	fraca(ig,l)=min(fraca(ig,l),0.5)
485	fracd(ig,l)=1.-fraca(ig,l)
486	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
487	else
488	c wa_moy(ig,l)=0.
489	fraca(ig,l)=0.
490	fracc(ig,l)=0.
491	fracd(ig,l)=1.
492	endif
493	enddo
494	enddo
495
496	c print*,'11 OK convect8'
497	c print*,'Ea3 ',wa_moy
498	c------------------------------------------------------------------
499	c Calcul de fracd, wd
500	c somme wa - wd = 0
501	c------------------------------------------------------------------
502
503
504	do ig=1,ngrid
505	fm(ig,1)=0.
506	fm(ig,nlay+1)=0.
507	enddo
508
509	do l=2,nlay
510	do ig=1,ngrid
511	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
512	enddo
513	do ig=1,ngrid
514	if(fracd(ig,l).lt.0.1) then
515	stop'fracd trop petit'
516	else
517	c vitesse descendante "diagnostique"
518	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
519	endif
520	enddo
521	enddo
522
523	do l=1,nlay
524	do ig=1,ngrid
525	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
526	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
527	enddo
528	enddo
529
530	c print*,'12 OK convect8'
531	c print*,'WA4 ',wa_moy
532	cc------------------------------------------------------------------
533	c calcul du transport vertical
534	c------------------------------------------------------------------
535
536	go to 4444
537	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
538	do l=2,nlay-1
539	do ig=1,ngrid
540	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
541	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
542	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
543	c s ,fm(ig,l+1)*ptimestep
544	c s ,' M=',masse(ig,l),masse(ig,l+1)
545	endif
546	enddo
547	enddo
548
549	do l=1,nlay
550	do ig=1,ngrid
551	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
552	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
553	c s ,entr(ig,l)*ptimestep
554	c s ,' M=',masse(ig,l)
555	endif
556	enddo
557	enddo
558
559	do l=1,nlay
560	do ig=1,ngrid
561	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
562	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
563	c s ,' FM=',fm(ig,l)
564	endif
565	if(.not.masse(ig,l).ge.1.e-10
566	s .or..not.masse(ig,l).le.1.e4) then
567	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
568	c s ,' M=',masse(ig,l)
569	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
570	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
571	c print*,'zlev(ig,l+1),zlev(ig,l)'
572	c s ,zlev(ig,l+1),zlev(ig,l)
573	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
574	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
575	endif
576	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
577	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
578	c s ,' E=',entr(ig,l)
579	endif
580	enddo
581	enddo
582
583	4444 continue
584
585	if (w2di.eq.1) then
586	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
587	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
588	else
589	fm0=fm
590	entr0=entr
591	endif
592
593	if (1.eq.1) then
594	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
595	. ,zh,zdhadj,zha)
596	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
597	. ,zo,pdoadj,zoa)
598	else
599	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
600	. ,zh,zdhadj,zha)
601	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
602	. ,zo,pdoadj,zoa)
603	endif
604
605	if (1.eq.0) then
606	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
607	. ,fraca,zmax
608	. ,zu,zv,pduadj,pdvadj,zua,zva)
609	else
610	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
611	. ,zu,pduadj,zua)
612	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
613	. ,zv,pdvadj,zva)
614	endif
615
616	do l=1,nlay
617	do ig=1,ngrid
618	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
619	zf2=zf/(1.-zf)
620	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
621	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
622	enddo
623	enddo
624
625
626
627	c print*,'13 OK convect8'
628	c print*,'WA5 ',wa_moy
629	do l=1,nlay
630	do ig=1,ngrid
631	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
632	enddo
633	enddo
634
635
636	c do l=1,nlay
637	c do ig=1,ngrid
638	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
639	c print*,'WARN!!! ig=',ig,' l=',l
640	c s ,' pdtadj=',pdtadj(ig,l)
641	c endif
642	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
643	c print*,'WARN!!! ig=',ig,' l=',l
644	c s ,' pdoadj=',pdoadj(ig,l)
645	c endif
646	c enddo
647	c enddo
648
649	c print*,'14 OK convect8'
650	c------------------------------------------------------------------
651	c Calculs pour les sorties
652	c------------------------------------------------------------------
653
654	if(sorties) then
655	do l=1,nlay
656	do ig=1,ngrid
657	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
658	zld(ig,l)=fracd(ig,l)*zmax(ig)
659	if(1.-fracd(ig,l).gt.1.e-10)
660	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
661	enddo
662	enddo
663
664	do l=1,nlay
665	do ig=1,ngrid
666	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
667	if (detr(ig,l).lt.0.) then
668	entr(ig,l)=entr(ig,l)-detr(ig,l)
669	detr(ig,l)=0.
670	c print*,'WARNING !!! detrainement negatif ',ig,l
671	endif
672	enddo
673	enddo
674
675	c print*,'15 OK convect8'
676
677	isplit=isplit+1
678
679
680	c #define und
681	goto 123
682	#ifdef und
683	CALL writeg1d(1,nlay,wd,'wd ','wd ')
684	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
685	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
686	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
687	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
688	CALL writeg1d(1,nlay,zla,'la ','la ')
689	CALL writeg1d(1,nlay,zld,'ld ','ld ')
690	CALL writeg1d(1,nlay,pt,'pt ','pt ')
691	CALL writeg1d(1,nlay,zh,'zh ','zh ')
692	CALL writeg1d(1,nlay,zha,'zha ','zha ')
693	CALL writeg1d(1,nlay,zu,'zu ','zu ')
694	CALL writeg1d(1,nlay,zv,'zv ','zv ')
695	CALL writeg1d(1,nlay,zo,'zo ','zo ')
696	CALL writeg1d(1,nlay,wh,'wh ','wh ')
697	CALL writeg1d(1,nlay,wu,'wu ','wu ')
698	CALL writeg1d(1,nlay,wv,'wv ','wv ')
699	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
700	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
701	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
702	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
703	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
704	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
705	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
706	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
707
708	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
709	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
710	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
711	c recalcul des flux en diagnostique...
712	c print*,'PAS DE TEMPS ',ptimestep
713	call dt2F(pplev,pplay,pt,pdtadj,wh)
714	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
715	#endif
716	123 continue
717	! #define troisD
718	#ifdef troisD
719	c if (sorties) then
720	print*,'Debut des wrgradsfi'
721
722	c print*,'16 OK convect8'
723	call wrgradsfi(1,nlay,wd,'wd ','wd ')
724	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
725	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
726	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
727	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
728	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
729	c print*,'WA6 ',wa_moy
730	call wrgradsfi(1,nlay,zla,'la ','la ')
731	call wrgradsfi(1,nlay,zld,'ld ','ld ')
732	call wrgradsfi(1,nlay,pt,'pt ','pt ')
733	call wrgradsfi(1,nlay,zh,'zh ','zh ')
734	call wrgradsfi(1,nlay,zha,'zha ','zha ')
735	call wrgradsfi(1,nlay,zua,'zua ','zua ')
736	call wrgradsfi(1,nlay,zva,'zva ','zva ')
737	call wrgradsfi(1,nlay,zu,'zu ','zu ')
738	call wrgradsfi(1,nlay,zv,'zv ','zv ')
739	call wrgradsfi(1,nlay,zo,'zo ','zo ')
740	call wrgradsfi(1,nlay,wh,'wh ','wh ')
741	call wrgradsfi(1,nlay,wu,'wu ','wu ')
742	call wrgradsfi(1,nlay,wv,'wv ','wv ')
743	call wrgradsfi(1,nlay,wo,'wo ','wo ')
744	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
745	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
746	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
747	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
748	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
749	call wrgradsfi(1,nlay,entr,'entr ','entr ')
750	call wrgradsfi(1,nlay,detr,'detr ','detr ')
751	call wrgradsfi(1,nlay,fm,'fm ','fm ')
752	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
753	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
754	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
755	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
756
757	call wrgradsfi(1,nlay,zo,'zo ','zo ')
758	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
759	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
760
761
762	c print*,'17 OK convect8'
763
764	do k=1,klev/10
765	write(str2,'(i2.2)') k
766	str10='wa'//str2
767	do l=1,nlay
768	do ig=1,ngrid
769	zsortie(ig,l)=wa(ig,k,l)
770	enddo
771	enddo
772	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
773	do l=1,nlay
774	do ig=1,ngrid
775	zsortie(ig,l)=larg_part(ig,k,l)
776	enddo
777	enddo
778	str10='la'//str2
779	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
780	enddo
781
782
783	c print*,'18 OK convect8'
784	c endif
785	print*,'Fin des wrgradsfi'
786	#endif
787
788	endif
789
790	c if(wa_moy(1,4).gt.1.e-10) stop
791
792	c print*,'19 OK convect8'
793	return
794	end
795
796	SUBROUTINE thermcell_cld(ngrid,nlay,ptimestep
797	s ,pplay,pplev,pphi,zlev,debut
798	s ,pu,pv,pt,po
799	s ,pduadj,pdvadj,pdtadj,pdoadj
800	s ,fm0,entr0,zqla,lmax
801	s ,zmax_sec,wmax_sec,zw_sec,lmix_sec
802	s ,ratqscth,ratqsdiff
803	c s ,pu_therm,pv_therm
804	s ,r_aspect,l_mix,w2di,tho)
805
806	USE dimphy
807	IMPLICIT NONE
808
809	c=======================================================================
810	c
811	c Calcul du transport verticale dans la couche limite en presence
812	c de "thermiques" explicitement representes
813	c
814	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
815	c
816	c le thermique est supposé homogène et dissipé par mélange avec
817	c son environnement. la longueur l_mix contrôle l'efficacité du
818	c mélange
819	c
820	c Le calcul du transport des différentes espèces se fait en prenant
821	c en compte:
822	c 1. un flux de masse montant
823	c 2. un flux de masse descendant
824	c 3. un entrainement
825	c 4. un detrainement
826	c
827	c=======================================================================
828
829	c-----------------------------------------------------------------------
830	c declarations:
831	c -------------
832
833	#include "dimensions.h"
834	cccc#include "dimphy.h"
835	#include "YOMCST.h"
836	#include "YOETHF.h"
837	#include "FCTTRE.h"
838
839	c arguments:
840	c ----------
841
842	INTEGER ngrid,nlay,w2di,tho
843	real ptimestep,l_mix,r_aspect
844	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
845	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
846	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
847	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
848	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
849	real pphi(ngrid,nlay)
850
851	integer idetr
852	save idetr
853	data idetr/3/
854
855	c local:
856	c ------
857
858	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
859	real zsortie1d(klon)
860	c CR: on remplace lmax(klon,klev+1)
861	INTEGER lmax(klon),lmin(klon),lentr(klon)
862	real linter(klon)
863	real zmix(klon), fracazmix(klon)
864	real alpha
865	save alpha
866	data alpha/1./
867	c RC
868	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
869	real zmax_sec(klon)
870	real zmax_sec2(klon)
871	real zw_sec(klon,klev+1)
872	INTEGER lmix_sec(klon)
873	real w_est(klon,klev+1)
874	con garde le zmax du pas de temps precedent
875	c real zmax0(klon)
876	c save zmax0
877	c real zmix0(klon)
878	c save zmix0
879	REAL, SAVE, ALLOCATABLE :: zmax0(:), zmix0(:)
880	c$OMP THREADPRIVATE(zmax0, zmix0)
881
882	real zlev(klon,klev+1),zlay(klon,klev)
883	real deltaz(klon,klev)
884	REAL zh(klon,klev),zdhadj(klon,klev)
885	real zthl(klon,klev),zdthladj(klon,klev)
886	REAL ztv(klon,klev)
887	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
888	real zl(klon,klev)
889	REAL wh(klon,klev+1)
890	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
891	real zla(klon,klev+1)
892	real zwa(klon,klev+1)
893	real zld(klon,klev+1)
894	real zwd(klon,klev+1)
895	real zsortie(klon,klev)
896	real zva(klon,klev)
897	real zua(klon,klev)
898	real zoa(klon,klev)
899
900	real zta(klon,klev)
901	real zha(klon,klev)
902	real wa_moy(klon,klev+1)
903	real fraca(klon,klev+1)
904	real fracc(klon,klev+1)
905	real zf,zf2
906	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
907	real q2(klon,klev)
908	real dtheta(klon,klev)
909	! common/comtherm/thetath2,wth2
910
911	real ratqscth(klon,klev)
912	real sum
913	real sumdiff
914	real ratqsdiff(klon,klev)
915	real count_time
916	integer isplit,nsplit,ialt
917	parameter (nsplit=10)
918	data isplit/0/
919	save isplit
920
921	logical sorties
922	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
923	real zpspsk(klon,klev)
924
925	c real wmax(klon,klev),wmaxa(klon)
926	real wmax(klon),wmaxa(klon)
927	real wmax_sec(klon)
928	real wmax_sec2(klon)
929	real wa(klon,klev,klev+1)
930	real wd(klon,klev+1)
931	real larg_part(klon,klev,klev+1)
932	real fracd(klon,klev+1)
933	real xxx(klon,klev+1)
934	real larg_cons(klon,klev+1)
935	real larg_detr(klon,klev+1)
936	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
937	real massetot(klon,klev)
938	real detr0(klon,klev)
939	real alim0(klon,klev)
940	real pu_therm(klon,klev),pv_therm(klon,klev)
941	real fm(klon,klev+1),entr(klon,klev)
942	real fmc(klon,klev+1)
943
944	real zcor,zdelta,zcvm5,qlbef
945	real Tbef(klon),qsatbef(klon)
946	real dqsat_dT,DT,num,denom
947	REAL REPS,RLvCp,DDT0
948	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
949	cCR niveau de condensation
950	real nivcon(klon)
951	real zcon(klon)
952	real zqsat(klon,klev)
953	real zqsatth(klon,klev)
954	PARAMETER (DDT0=.01)
955
956
957	cCR:nouvelles variables
958	real f_star(klon,klev+1),entr_star(klon,klev)
959	real detr_star(klon,klev)
960	real alim_star_tot(klon),alim_star2(klon)
961	real entr_star_tot(klon)
962	real detr_star_tot(klon)
963	real alim_star(klon,klev)
964	real alim(klon,klev)
965	real nu(klon,klev)
966	real nu_e(klon,klev)
967	real nu_min
968	real nu_max
969	real nu_r
970	real f(klon)
971	c real f(klon), f0(klon)
972	c save f0
973	REAL,SAVE, ALLOCATABLE :: f0(:)
974	c$OMP THREADPRIVATE(f0)
975
976	real f_old
977	real zlevinter(klon)
978	logical, save :: first = .true.
979	c data first /.false./
980	c save first
981	logical nuage
982	c save nuage
983	logical boucle
984	logical therm
985	logical debut
986	logical rale
987	integer test(klon)
988	integer signe_zw2
989	cRC
990
991	character*2 str2
992	character*10 str10
993
994	LOGICAL vtest(klon),down
995	LOGICAL Zsat(klon)
996
997	EXTERNAL SCOPY
998
999	integer ncorrec,ll
1000	save ncorrec
1001	data ncorrec/0/
1002	c
1003
1004	c-----------------------------------------------------------------------
1005	c initialisation:
1006	c ---------------
1007	c
1008	if (first) then
1009	allocate(zmix0(klon))
1010	allocate(zmax0(klon))
1011	allocate(f0(klon))
1012	first=.false.
1013	endif
1014
1015	sorties=.false.
1016	c print*,'NOUVEAU DETR PLUIE '
1017	IF(ngrid.NE.klon) THEN
1018	PRINT*
1019	PRINT*,'STOP dans convadj'
1020	PRINT*,'ngrid =',ngrid
1021	PRINT*,'klon =',klon
1022	ENDIF
1023	c
1024	c Initialisation
1025	RLvCp = RLVTT/RCPD
1026	REPS = RD/RV
1027	cinitialisations de zqsat
1028	DO ll=1,nlay
1029	DO ig=1,ngrid
1030	zqsat(ig,ll)=0.
1031	zqsatth(ig,ll)=0.
1032	ENDDO
1033	ENDDO
1034	c
1035	con met le first a true pour le premier passage de la journée
1036	do ig=1,klon
1037	test(ig)=0
1038	enddo
1039	if (debut) then
1040	do ig=1,klon
1041	test(ig)=1
1042	f0(ig)=0.
1043	zmax0(ig)=0.
1044	enddo
1045	endif
1046	do ig=1,klon
1047	if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
1048	test(ig)=1
1049	endif
1050	enddo
1051	c do ig=1,klon
1052	c print*,'test(ig)',test(ig),zmax0(ig)
1053	c enddo
1054	nuage=.false.
1055	c-----------------------------------------------------------------------
1056	cAM Calcul de T,q,ql a partir de Tl et qT
1057	c ---------------------------------------------------
1058	c
1059	c Pr Tprec=Tl calcul de qsat
1060	c Si qsat>qT T=Tl, q=qT
1061	c Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
1062	c On cherche DDT < DDT0
1063	c
1064	c defaut
1065	DO ll=1,nlay
1066	DO ig=1,ngrid
1067	zo(ig,ll)=po(ig,ll)
1068	zl(ig,ll)=0.
1069	zh(ig,ll)=pt(ig,ll)
1070	EndDO
1071	EndDO
1072	do ig=1,ngrid
1073	Zsat(ig)=.false.
1074	enddo
1075	c
1076	c
1077	DO ll=1,nlay
1078	c les points insatures sont definitifs
1079	DO ig=1,ngrid
1080	Tbef(ig)=pt(ig,ll)
1081	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1082	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
1083	qsatbef(ig)=MIN(0.5,qsatbef(ig))
1084	zcor=1./(1.-retv*qsatbef(ig))
1085	qsatbef(ig)=qsatbef(ig)*zcor
1086	Zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig)) .gt. 1.e-10)
1087	EndDO
1088
1089	DO ig=1,ngrid
1090	if (Zsat(ig).and.(1.eq.1)) then
1091	qlbef=max(0.,po(ig,ll)-qsatbef(ig))
1092	c si sature: ql est surestime, d'ou la sous-relax
1093	DT = 0.5RLvCpqlbef
1094	c write(18,*),'DT0=',DT
1095	c on pourra enchainer 2 ou 3 calculs sans Do while
1096	do while (abs(DT).gt.DDT0)
1097	c il faut verifier si c,a conserve quand on repasse en insature ...
1098	Tbef(ig)=Tbef(ig)+DT
1099	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1100	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
1101	qsatbef(ig)=MIN(0.5,qsatbef(ig))
1102	zcor=1./(1.-retv*qsatbef(ig))
1103	qsatbef(ig)=qsatbef(ig)*zcor
1104	c on veut le signe de qlbef
1105	qlbef=po(ig,ll)-qsatbef(ig)
1106	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1107	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
1108	zcor=1./(1.-retv*qsatbef(ig))
1109	dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
1110	num=-Tbef(ig)+pt(ig,ll)+RLvCp*qlbef
1111	denom=1.+RLvCp*dqsat_dT
1112	if (denom.lt.1.e-10) then
1113	print*,'pb denom'
1114	endif
1115	DT=num/denom
1116	enddo
1117	c on ecrit de maniere conservative (sat ou non)
1118	zl(ig,ll) = max(0.,qlbef)
1119	c T = Tl +Lv/Cp ql
1120	zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
1121	zo(ig,ll) = po(ig,ll)-zl(ig,ll)
1122	endif
1123	con ecrit zqsat
1124	zqsat(ig,ll)=qsatbef(ig)
1125	EndDO
1126	EndDO
1127	cAM fin
1128	c
1129	c-----------------------------------------------------------------------
1130	c incrementation eventuelle de tendances precedentes:
1131	c ---------------------------------------------------
1132
1133	c print*,'0 OK convect8'
1134
1135	DO 1010 l=1,nlay
1136	DO 1015 ig=1,ngrid
1137	zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
1138	c zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
1139	c zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
1140	zu(ig,l)=pu(ig,l)
1141	zv(ig,l)=pv(ig,l)
1142	c zo(ig,l)=po(ig,l)
1143	c ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
1144	cAM attention zh est maintenant le profil de T et plus le profil de theta !
1145	c
1146	c T-> Theta
1147	ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
1148	cAM Theta_v
1149	ztv(ig,l)=ztv(ig,l)(1.+RETV(zo(ig,l))
1150	s -zl(ig,l))
1151	cAM Thetal
1152	zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
1153	c
1154	1015 CONTINUE
1155	1010 CONTINUE
1156
1157	c print*,'1 OK convect8'
1158	c --------------------
1159	c
1160	c
1161	c + + + + + + + + + + +
1162	c
1163	c
1164	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
1165	c wh,wt,wo ...
1166	c
1167	c + + + + + + + + + + + zh,zu,zv,zo,rho
1168	c
1169	c
1170	c -------------------- zlev(1)
1171	c \\\\\\\\\\\\\\\\\\\\
1172	c
1173	c
1174
1175	c-----------------------------------------------------------------------
1176	c Calcul des altitudes des couches
1177	c-----------------------------------------------------------------------
1178
1179	do l=2,nlay
1180	do ig=1,ngrid
1181	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
1182	enddo
1183	enddo
1184	do ig=1,ngrid
1185	zlev(ig,1)=0.
1186	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
1187	enddo
1188	do l=1,nlay
1189	do ig=1,ngrid
1190	zlay(ig,l)=pphi(ig,l)/RG
1191	enddo
1192	enddo
1193	ccalcul de deltaz
1194	do l=1,nlay
1195	do ig=1,ngrid
1196	deltaz(ig,l)=zlev(ig,l+1)-zlev(ig,l)
1197	enddo
1198	enddo
1199
1200	c print*,'2 OK convect8'
1201	c-----------------------------------------------------------------------
1202	c Calcul des densites
1203	c-----------------------------------------------------------------------
1204
1205	do l=1,nlay
1206	do ig=1,ngrid
1207	c rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
1208	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDztv(ig,l))
1209	enddo
1210	enddo
1211
1212	do l=2,nlay
1213	do ig=1,ngrid
1214	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
1215	enddo
1216	enddo
1217
1218	do k=1,nlay
1219	do l=1,nlay+1
1220	do ig=1,ngrid
1221	wa(ig,k,l)=0.
1222	enddo
1223	enddo
1224	enddo
1225	cCr:ajout:calcul de la masse
1226	do l=1,nlay
1227	do ig=1,ngrid
1228	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1229	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
1230	enddo
1231	enddo
1232	c print*,'3 OK convect8'
1233	c------------------------------------------------------------------
1234	c Calcul de w2, quarre de w a partir de la cape
1235	c a partir de w2, on calcule wa, vitesse de l'ascendance
1236	c
1237	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
1238	c w2 est stoke dans wa
1239	c
1240	c ATTENTION: dans convect8, on n'utilise le calcule des wa
1241	c independants par couches que pour calculer l'entrainement
1242	c a la base et la hauteur max de l'ascendance.
1243	c
1244	c Indicages:
1245	c l'ascendance provenant du niveau k traverse l'interface l avec
1246	c une vitesse wa(k,l).
1247	c
1248	c --------------------
1249	c
1250	c + + + + + + + + + +
1251	c
1252	c wa(k,l) ---- -------------------- l
1253	c /\
1254	c /\|\|\ + + + + + + + + + +
1255	c \|\|
1256	c \|\| --------------------
1257	c \|\|
1258	c \|\| + + + + + + + + + +
1259	c \|\|
1260	c \|\| --------------------
1261	c \|\|__
1262	c \|___ + + + + + + + + + + k
1263	c
1264	c --------------------
1265	c
1266	c
1267	c
1268	c------------------------------------------------------------------
1269
1270	cCR: ponderation entrainement des couches instables
1271	cdef des alim_star tels que alim=f*alim_star
1272	do l=1,klev
1273	do ig=1,ngrid
1274	alim_star(ig,l)=0.
1275	alim(ig,l)=0.
1276	enddo
1277	enddo
1278	c determination de la longueur de la couche d entrainement
1279	do ig=1,ngrid
1280	lentr(ig)=1
1281	enddo
1282
1283	con ne considere que les premieres couches instables
1284	therm=.false.
1285	do k=nlay-2,1,-1
1286	do ig=1,ngrid
1287	if (ztv(ig,k).gt.ztv(ig,k+1).and.
1288	s ztv(ig,k+1).le.ztv(ig,k+2)) then
1289	lentr(ig)=k+1
1290	therm=.true.
1291	endif
1292	enddo
1293	enddo
1294	c
1295	c determination du lmin: couche d ou provient le thermique
1296	do ig=1,ngrid
1297	lmin(ig)=1
1298	enddo
1299	do ig=1,ngrid
1300	do l=nlay,2,-1
1301	if (ztv(ig,l-1).gt.ztv(ig,l)) then
1302	lmin(ig)=l-1
1303	endif
1304	enddo
1305	enddo
1306	c
1307	c definition de l'entrainement des couches
1308	do l=1,klev-1
1309	do ig=1,ngrid
1310	if (ztv(ig,l).gt.ztv(ig,l+1).and.
1311	s l.ge.lmin(ig).and.l.lt.lentr(ig)) then
1312	cdef possibles pour alim_star: zdthetadz, dthetadz, zdtheta
1313	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)
1314	c s *(zlev(ig,l+1)-zlev(ig,l))
1315	s *sqrt(zlev(ig,l+1))
1316	c alim_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
1317	c s /zlev(ig,lentr(ig)+2)))**(3./2.)
1318	endif
1319	enddo
1320	enddo
1321
1322	c pas de thermique si couche 1 stable
1323	do ig=1,ngrid
1324	c if (lmin(ig).gt.1) then
1325	cCRnouveau test
1326	if (alim_star(ig,1).lt.1.e-10) then
1327	do l=1,klev
1328	alim_star(ig,l)=0.
1329	enddo
1330	endif
1331	enddo
1332	c calcul de l entrainement total
1333	do ig=1,ngrid
1334	alim_star_tot(ig)=0.
1335	entr_star_tot(ig)=0.
1336	detr_star_tot(ig)=0.
1337	enddo
1338	do ig=1,ngrid
1339	do k=1,klev
1340	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,k)
1341	enddo
1342	enddo
1343	c
1344	c Calcul entrainement normalise
1345	do ig=1,ngrid
1346	if (alim_star_tot(ig).gt.1.e-10) then
1347	c do l=1,lentr(ig)
1348	do l=1,klev
1349	cdef possibles pour entr_star: zdthetadz, dthetadz, zdtheta
1350	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
1351	enddo
1352	endif
1353	enddo
1354
1355	c print*,'fin calcul alim_star'
1356
1357	cAM:initialisations
1358	do k=1,nlay
1359	do ig=1,ngrid
1360	ztva(ig,k)=ztv(ig,k)
1361	ztla(ig,k)=zthl(ig,k)
1362	zqla(ig,k)=0.
1363	zqta(ig,k)=po(ig,k)
1364	Zsat(ig) =.false.
1365	enddo
1366	enddo
1367	do k=1,klev
1368	do ig=1,ngrid
1369	detr_star(ig,k)=0.
1370	entr_star(ig,k)=0.
1371	detr(ig,k)=0.
1372	entr(ig,k)=0.
1373	enddo
1374	enddo
1375	c print*,'7 OK convect8'
1376	do k=1,klev+1
1377	do ig=1,ngrid
1378	zw2(ig,k)=0.
1379	fmc(ig,k)=0.
1380	cCR
1381	f_star(ig,k)=0.
1382	cRC
1383	larg_cons(ig,k)=0.
1384	larg_detr(ig,k)=0.
1385	wa_moy(ig,k)=0.
1386	enddo
1387	enddo
1388
1389	cn print*,'8 OK convect8'
1390	do ig=1,ngrid
1391	linter(ig)=1.
1392	lmaxa(ig)=1
1393	lmix(ig)=1
1394	wmaxa(ig)=0.
1395	enddo
1396
1397	nu_min=l_mix
1398	nu_max=1000.
1399	c do ig=1,ngrid
1400	c nu_max=wmax_sec(ig)
1401	c enddo
1402	do ig=1,ngrid
1403	do k=1,klev
1404	nu(ig,k)=0.
1405	nu_e(ig,k)=0.
1406	enddo
1407	enddo
1408	cCalcul de l'excès de température du à la diffusion turbulente
1409	do ig=1,ngrid
1410	do l=1,klev
1411	dtheta(ig,l)=0.
1412	enddo
1413	enddo
1414	do ig=1,ngrid
1415	do l=1,lentr(ig)-1
1416	dtheta(ig,l)=sqrt(10.0.4zlev(ig,l+1)*21.
1417	s ((ztv(ig,l+1)-ztv(ig,l))/(zlev(ig,l+1)-zlev(ig,l)))*2)
1418	enddo
1419	enddo
1420	c do l=1,nlay-2
1421	do l=1,klev-1
1422	do ig=1,ngrid
1423	if (ztv(ig,l).gt.ztv(ig,l+1)
1424	s .and.alim_star(ig,l).gt.1.e-10
1425	s .and.zw2(ig,l).lt.1e-10) then
1426	cAM
1427	ctest:on rajoute un excès de T dans couche alim
1428	c ztla(ig,l)=zthl(ig,l)+dtheta(ig,l)
1429	ztla(ig,l)=zthl(ig,l)
1430	ctest: on rajoute un excès de q dans la couche alim
1431	c zqta(ig,l)=po(ig,l)+0.001
1432	zqta(ig,l)=po(ig,l)
1433	zqla(ig,l)=zl(ig,l)
1434	cAM
1435	f_star(ig,l+1)=alim_star(ig,l)
1436	ctest:calcul de dteta
1437	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
1438	s *(zlev(ig,l+1)-zlev(ig,l))
1439	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1440	w_est(ig,l+1)=zw2(ig,l+1)
1441	larg_detr(ig,l)=0.
1442	c print*,'coucou boucle 1'
1443	else if ((zw2(ig,l).ge.1e-10).and.
1444	s (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
1445	c print*,'coucou boucle 2'
1446	cestimation du detrainement a partir de la geometrie du pas precedent
1447	if ((test(ig).eq.1).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
1448	detr_star(ig,l)=0.
1449	entr_star(ig,l)=0.
1450	c print*,'coucou test(ig)',test(ig),f0(ig),zmax0(ig)
1451	else
1452	c print*,'coucou debut detr'
1453	ctests sur la definition du detr
1454	if (zqla(ig,l-1).gt.1.e-10) then
1455	nuage=.true.
1456	endif
1457
1458	w_est(ig,l+1)=zw2(ig,l)*
1459	s ((f_star(ig,l))**2)
1460	s /(f_star(ig,l)+alim_star(ig,l))**2+
1461	s 2.RG(ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l)
1462	s *(zlev(ig,l+1)-zlev(ig,l))
1463	if (w_est(ig,l+1).lt.0.) then
1464	w_est(ig,l+1)=zw2(ig,l)
1465	endif
1466	if (l.gt.2) then
1467	if ((w_est(ig,l+1).gt.w_est(ig,l)).and.
1468	s (zlev(ig,l+1).lt.zmax_sec(ig)).and.
1469	s (zqla(ig,l-1).lt.1.e-10)) then
1470	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)
1471	s sqrt(w_est(ig,l+1))sqrt(nu(ig,l)*zlev(ig,l+1))
1472	s -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(nu(ig,l)*zlev(ig,l)))
1473	s /(r_aspect*zmax_sec(ig)))
1474	else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.
1475	s (zqla(ig,l-1).lt.1.e-10)) then
1476	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))
1477	s /(rhobarz(ig,lmix(ig))wmaxa(ig))
1478	s (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))
1479	s *((zmax_sec(ig)-zlev(ig,l+1))/((zmax_sec(ig)-zlev(ig,lmix(ig)))))
1480	s **2.
1481	s -rhobarz(ig,l)*sqrt(w_est(ig,l))
1482	s *((zmax_sec(ig)-zlev(ig,l))/((zmax_sec(ig)-zlev(ig,lmix(ig)))))
1483	s **2.)
1484	else
1485	detr_star(ig,l)=0.002f0(ig)f_star(ig,l)
1486	s *(zlev(ig,l+1)-zlev(ig,l))
1487
1488	endif
1489	else
1490	detr_star(ig,l)=0.
1491	endif
1492
1493	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
1494	if (nuage) then
1495	entr_star(ig,l)=0.4*detr_star(ig,l)
1496	else
1497	entr_star(ig,l)=0.4*detr_star(ig,l)
1498	endif
1499
1500	if ((detr_star(ig,l)).gt.f_star(ig,l)) then
1501	detr_star(ig,l)=f_star(ig,l)
1502	c entr_star(ig,l)=0.
1503	endif
1504
1505	if ((l.lt.lentr(ig))) then
1506	entr_star(ig,l)=0.
1507	c detr_star(ig,l)=0.
1508	endif
1509
1510	c print*,'ok detr_star'
1511	endif
1512	cprise en compte du detrainement dans le calcul du flux
1513	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
1514	s -detr_star(ig,l)
1515	ctest
1516	c if (f_star(ig,l+1).lt.0.) then
1517	c f_star(ig,l+1)=0.
1518	c entr_star(ig,l)=0.
1519	c detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)
1520	c endif
1521	ctest sur le signe de f_star
1522	if (f_star(ig,l+1).gt.1.e-10) then
1523	c then
1524	ctest
1525	c if (((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10)) then
1526	cAM on melange Tl et qt du thermique
1527	con rajoute un excès de T dans la couche alim
1528	c if (l.lt.lentr(ig)) then
1529	c ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1530	c s (alim_star(ig,l)+entr_star(ig,l))*(zthl(ig,l)+dtheta(ig,l)))
1531	c s /(f_star(ig,l+1)+detr_star(ig,l))
1532	c else
1533	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1534	s (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))
1535	s /(f_star(ig,l+1)+detr_star(ig,l))
1536	c s /(f_star(ig,l+1))
1537	c endif
1538	con rajoute un excès de q dans la couche alim
1539	c if (l.lt.lentr(ig)) then
1540	c zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1541	c s (alim_star(ig,l)+entr_star(ig,l))*(po(ig,l)+0.001))
1542	c s /(f_star(ig,l+1)+detr_star(ig,l))
1543	c else
1544	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1545	s (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))
1546	s /(f_star(ig,l+1)+detr_star(ig,l))
1547	c s /(f_star(ig,l+1))
1548	c endif
1549	cAM on en deduit thetav et ql du thermique
1550	cCR test
1551	c Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1552	Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1553	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1554	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
1555	qsatbef(ig)=MIN(0.5,qsatbef(ig))
1556	zcor=1./(1.-retv*qsatbef(ig))
1557	qsatbef(ig)=qsatbef(ig)*zcor
1558	Zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig)) .gt. 1.e-10)
1559
1560	if (Zsat(ig).and.(1.eq.1)) then
1561	qlbef=max(0.,zqta(ig,l)-qsatbef(ig))
1562	DT = 0.5RLvCpqlbef
1563	c write(17,*)'DT0=',DT
1564	do while (abs(DT).gt.DDT0)
1565	c print*,'aie'
1566	Tbef(ig)=Tbef(ig)+DT
1567	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1568	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
1569	qsatbef(ig)=MIN(0.5,qsatbef(ig))
1570	zcor=1./(1.-retv*qsatbef(ig))
1571	qsatbef(ig)=qsatbef(ig)*zcor
1572	qlbef=zqta(ig,l)-qsatbef(ig)
1573
1574	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
1575	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
1576	zcor=1./(1.-retv*qsatbef(ig))
1577	dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
1578	num=-Tbef(ig)+ztla(ig,l)zpspsk(ig,l)+RLvCpqlbef
1579	denom=1.+RLvCp*dqsat_dT
1580	if (denom.lt.1.e-10) then
1581	print*,'pb denom'
1582	endif
1583	DT=num/denom
1584	c write(17,*)'DT=',DT
1585	enddo
1586	zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1587	zqla(ig,l) = max(0.,qlbef)
1588	c zqla(ig,l)=0.
1589	endif
1590	c zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1591	c
1592	c on ecrit de maniere conservative (sat ou non)
1593	c T = Tl +Lv/Cp ql
1594	cCR rq utilisation de humidite specifique ou rapport de melange?
1595	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
1596	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1597	con rajoute le calcul de zha pour diagnostiques (temp potentielle)
1598	zha(ig,l) = ztva(ig,l)
1599	c if (l.lt.lentr(ig)) then
1600	c ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)
1601	c s -zqla(ig,l))-zqla(ig,l)) + 0.1
1602	c else
1603	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)
1604	s -zqla(ig,l))-zqla(ig,l))
1605	c endif
1606	c ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
1607	c s /(1.-retv*zqla(ig,l))
1608	c ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1609	c ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)
1610	c s /(1.-retv*zqta(ig,l))
1611	c s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1612	c s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1613	c write(13,)zqla(ig,l),zqla(ig,l)/(1.-retvzqla(ig,l))
1614	con ecrit zqsat
1615	zqsatth(ig,l)=qsatbef(ig)
1616	c enddo
1617	c DO ig=1,ngrid
1618	c if (zw2(ig,l).ge.1.e-10.and.
1619	c s f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
1620	c mise a jour de la vitesse ascendante (l'air entraine de la couche
1621	c consideree commence avec une vitesse nulle).
1622	c
1623	c if (f_star(ig,l+1).gt.1.e-10) then
1624	zw2(ig,l+1)=zw2(ig,l)*
1625	c s ((f_star(ig,l)-detr_star(ig,l))**2)
1626	c s /f_star(ig,l+1)**2+
1627	s ((f_star(ig,l))**2)
1628	s /(f_star(ig,l+1)+detr_star(ig,l))**2+
1629	c s /(f_star(ig,l+1))**2+
1630	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1631	s *(zlev(ig,l+1)-zlev(ig,l))
1632	c s (f_star(ig,l)/f_star(ig,l+1))*2
1633
1634	endif
1635	endif
1636	c
1637	if (zw2(ig,l+1).lt.0.) then
1638	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
1639	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1640	zw2(ig,l+1)=0.
1641	c print*,'linter=',linter(ig)
1642	c else if ((zw2(ig,l+1).lt.1.e-10).and.(zw2(ig,l+1).ge.0.)) then
1643	c linter(ig)=l+1
1644	c print*,'linter=l',zw2(ig,l),zw2(ig,l+1)
1645	else
1646	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1647	c wa_moy(ig,l+1)=zw2(ig,l+1)
1648	endif
1649	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1650	c lmix est le niveau de la couche ou w (wa_moy) est maximum
1651	lmix(ig)=l+1
1652	wmaxa(ig)=wa_moy(ig,l+1)
1653	endif
1654	enddo
1655	enddo
1656	print*,'fin calcul zw2'
1657	c
1658	c Calcul de la couche correspondant a la hauteur du thermique
1659	do ig=1,ngrid
1660	lmax(ig)=lentr(ig)
1661	enddo
1662	do ig=1,ngrid
1663	do l=nlay,lentr(ig)+1,-1
1664	if (zw2(ig,l).le.1.e-10) then
1665	lmax(ig)=l-1
1666	endif
1667	enddo
1668	enddo
1669	c pas de thermique si couche 1 stable
1670	do ig=1,ngrid
1671	if (lmin(ig).gt.1) then
1672	lmax(ig)=1
1673	lmin(ig)=1
1674	lentr(ig)=1
1675	endif
1676	enddo
1677	c
1678	c Determination de zw2 max
1679	do ig=1,ngrid
1680	wmax(ig)=0.
1681	enddo
1682
1683	do l=1,nlay
1684	do ig=1,ngrid
1685	if (l.le.lmax(ig)) then
1686	if (zw2(ig,l).lt.0.)then
1687	print*,'pb2 zw2<0'
1688	endif
1689	zw2(ig,l)=sqrt(zw2(ig,l))
1690	wmax(ig)=max(wmax(ig),zw2(ig,l))
1691	else
1692	zw2(ig,l)=0.
1693	endif
1694	enddo
1695	enddo
1696
1697	c Longueur caracteristique correspondant a la hauteur des thermiques.
1698	do ig=1,ngrid
1699	zmax(ig)=0.
1700	zlevinter(ig)=zlev(ig,1)
1701	enddo
1702	do ig=1,ngrid
1703	c calcul de zlevinter
1704	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
1705	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
1706	s -zlev(ig,lmax(ig)))
1707	cpour le cas ou on prend tjs lmin=1
1708	c zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1709	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
1710	zmax0(ig)=zmax(ig)
1711	write(11,*)'ig,lmax,linter',ig,lmax(ig),linter(ig)
1712	write(12,*)'ig,zlevinter,zmax',ig,zmax(ig),zlevinter(ig)
1713	enddo
1714
1715	cCalcul de zmax_sec et wmax_sec
1716	call fermeture_seche(ngrid,nlay
1717	s ,pplay,pplev,pphi,zlev,rhobarz,f0,zpspsk
1718	s ,alim,zh,zo,lentr,lmin,nu_min,nu_max,r_aspect
1719	s ,zmax_sec2,wmax_sec2)
1720
1721	print*,'avant fermeture'
1722	c Fermeture,determination de f
1723	c en lmax f=d-e
1724	do ig=1,ngrid
1725	c entr_star(ig,lmax(ig))=0.
1726	c f_star(ig,lmax(ig)+1)=0.
1727	c detr_star(ig,lmax(ig))=f_star(ig,lmax(ig))+entr_star(ig,lmax(ig))
1728	c s +alim_star(ig,lmax(ig))
1729	enddo
1730	c
1731	do ig=1,ngrid
1732	alim_star2(ig)=0.
1733	enddo
1734	ccalcul de entr_star_tot
1735	do ig=1,ngrid
1736	do k=1,lmix(ig)
1737	entr_star_tot(ig)=entr_star_tot(ig)
1738	c s +entr_star(ig,k)
1739	s +alim_star(ig,k)
1740	c s -detr_star(ig,k)
1741	detr_star_tot(ig)=detr_star_tot(ig)
1742	c s +alim_star(ig,k)
1743	s -detr_star(ig,k)
1744	s +entr_star(ig,k)
1745	enddo
1746	enddo
1747
1748	do ig=1,ngrid
1749	if (alim_star_tot(ig).LT.1.e-10) then
1750	f(ig)=0.
1751	else
1752	c do k=lmin(ig),lentr(ig)
1753	do k=1,lentr(ig)
1754	alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2
1755	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
1756	enddo
1757	if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then
1758	f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect
1759	s *alim_star2(ig))
1760	f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/
1761	s zmax_sec(ig))*wmax_sec(ig))
1762	else
1763	f(ig)=wmax(ig)/(max(500.,zmax(ig))r_aspectalim_star2(ig))
1764	f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/
1765	s zmax(ig))*wmax(ig))
1766	endif
1767	endif
1768	f0(ig)=f(ig)
1769	enddo
1770	print*,'apres fermeture'
1771	c Calcul de l'entrainement
1772	do ig=1,ngrid
1773	do k=1,klev
1774	alim(ig,k)=f(ig)*alim_star(ig,k)
1775	enddo
1776	enddo
1777	cCR:test pour entrainer moins que la masse
1778	c do ig=1,ngrid
1779	c do l=1,lentr(ig)
1780	c if ((alim(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
1781	c alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
1782	c s -0.9*masse(ig,l)/ptimestep
1783	c alim(ig,l)=0.9*masse(ig,l)/ptimestep
1784	c endif
1785	c enddo
1786	c enddo
1787	c calcul du détrainement
1788	do ig=1,klon
1789	do k=1,klev
1790	detr(ig,k)=f(ig)*detr_star(ig,k)
1791	if (detr(ig,k).lt.0.) then
1792	c print*,'detr1<0!!!'
1793	endif
1794	enddo
1795	do k=1,klev
1796	entr(ig,k)=f(ig)*entr_star(ig,k)
1797	if (entr(ig,k).lt.0.) then
1798	c print*,'entr1<0!!!'
1799	endif
1800	enddo
1801	enddo
1802	c
1803	c do ig=1,ngrid
1804	c do l=1,klev
1805	c if (((detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep).gt.
1806	c s (masse(ig,l))) then
1807	c print*,'d2+e2+a2>m2','ig=',ig,'l=',l,'lmax(ig)=',lmax(ig),'d+e+a='
1808	c s,(detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep,'m=',masse(ig,l)
1809	c endif
1810	c enddo
1811	c enddo
1812	c Calcul des flux
1813
1814	do ig=1,ngrid
1815	do l=1,lmax(ig)
1816	c do l=1,klev
1817	c fmc(ig,l+1)=f(ig)*f_star(ig,l+1)
1818	fmc(ig,l+1)=fmc(ig,l)+alim(ig,l)+entr(ig,l)-detr(ig,l)
1819	c print*,'??!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1820	c s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1821	c s 'f+1=',fmc(ig,l+1)
1822	if (fmc(ig,l+1).lt.0.) then
1823	print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1824	fmc(ig,l+1)=fmc(ig,l)
1825	detr(ig,l)=alim(ig,l)+entr(ig,l)
1826	c fmc(ig,l+1)=0.
1827	c print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1828	endif
1829	c if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1830	c f_old=fmc(ig,l+1)
1831	c fmc(ig,l+1)=fmc(ig,l)
1832	c detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1833	c endif
1834
1835	c if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1836	c f_old=fmc(ig,l+1)
1837	c fmc(ig,l+1)=fmc(ig,l)
1838	c detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l)
1839	c endif
1840	crajout du test sur alpha croissant
1841	cif test
1842	c if (1.eq.0) then
1843
1844	if (l.eq.klev) then
1845	print*,'THERMCELL PB ig=',ig,' l=',l
1846	stop
1847	endif
1848	! if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1849	! s (l.ge.lentr(ig)).and.
1850	if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1851	s (l.ge.lentr(ig)) ) then
1852	if ( ((fmc(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.
1853	s (fmc(ig,l)/(rhobarz(ig,l)*zw2(ig,l))))) then
1854	f_old=fmc(ig,l+1)
1855	fmc(ig,l+1)=fmc(ig,l)rhobarz(ig,l+1)zw2(ig,l+1)
1856	s /(rhobarz(ig,l)*zw2(ig,l))
1857	detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1858	c detr(ig,l)=(fmc(ig,l+1)-fmc(ig,l))/(0.4-1.)
1859	c entr(ig,l)=0.4*detr(ig,l)
1860	c entr(ig,l)=fmc(ig,l+1)-fmc(ig,l)+detr(ig,l)
1861	endif
1862	endif
1863	if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1864	f_old=fmc(ig,l+1)
1865	fmc(ig,l+1)=fmc(ig,l)
1866	detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1867	endif
1868	if (detr(ig,l).gt.fmc(ig,l)) then
1869	detr(ig,l)=fmc(ig,l)
1870	entr(ig,l)=fmc(ig,l+1)-alim(ig,l)
1871	endif
1872	if (fmc(ig,l+1).lt.0.) then
1873	detr(ig,l)=detr(ig,l)+fmc(ig,l+1)
1874	fmc(ig,l+1)=0.
1875	print*,'fmc2<0',l+1,lmax(ig)
1876	endif
1877
1878	ctest pour ne pas avoir f=0 et d=e/=0
1879	c if (fmc(ig,l+1).lt.1.e-10) then
1880	c detr(ig,l+1)=0.
1881	c entr(ig,l+1)=0.
1882	c zqla(ig,l+1)=0.
1883	c zw2(ig,l+1)=0.
1884	c lmax(ig)=l+1
1885	c zmax(ig)=zlev(ig,lmax(ig))
1886	c endif
1887	if (zw2(ig,l+1).gt.1.e-10) then
1888	if ((((fmc(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.
1889	s 1.)) then
1890	f_old=fmc(ig,l+1)
1891	fmc(ig,l+1)=rhobarz(ig,l+1)*zw2(ig,l+1)
1892	zw2(ig,l+1)=0.
1893	zqla(ig,l+1)=0.
1894	detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1895	lmax(ig)=l+1
1896	zmax(ig)=zlev(ig,lmax(ig))
1897	print*,'alpha>1',l+1,lmax(ig)
1898	endif
1899	endif
1900	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
1901	cendif test
1902	c endif
1903	enddo
1904	enddo
1905	do ig=1,ngrid
1906	c if (fmc(ig,lmax(ig)+1).ne.0.) then
1907	fmc(ig,lmax(ig)+1)=0.
1908	entr(ig,lmax(ig))=0.
1909	detr(ig,lmax(ig))=fmc(ig,lmax(ig))+entr(ig,lmax(ig))
1910	s +alim(ig,lmax(ig))
1911	c endif
1912	enddo
1913	ctest sur le signe de fmc
1914	do ig=1,ngrid
1915	do l=1,klev+1
1916	if (fmc(ig,l).lt.0.) then
1917	print*,'fm1<0!!!','ig=',ig,'l=',l,'a=',alim(ig,l-1),'e='
1918	s ,entr(ig,l-1),'f=',fmc(ig,l-1),'d=',detr(ig,l-1),'f+1=',fmc(ig,l)
1919	endif
1920	enddo
1921	enddo
1922	ctest de verification
1923	do ig=1,ngrid
1924	do l=1,lmax(ig)
1925	if ((abs(fmc(ig,l+1)-fmc(ig,l)-alim(ig,l)-entr(ig,l)+detr(ig,l)))
1926	s .gt.1.e-4) then
1927	c print*,'pbcm!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1928	c s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1929	c s 'f+1=',fmc(ig,l+1)
1930	endif
1931	if (detr(ig,l).lt.0.) then
1932	print*,'detrdemi<0!!!'
1933	endif
1934	enddo
1935	enddo
1936	c
1937	cRC
1938	cCR def de zmix continu (profil parabolique des vitesses)
1939	do ig=1,ngrid
1940	if (lmix(ig).gt.1.) then
1941	c test
1942	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
1943	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
1944	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
1945	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
1946	s then
1947	c
1948	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
1949	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
1950	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
1951	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
1952	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
1953	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
1954	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
1955	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
1956	else
1957	zmix(ig)=zlev(ig,lmix(ig))
1958	print*,'pb zmix'
1959	endif
1960	else
1961	zmix(ig)=0.
1962	endif
1963	ctest
1964	if ((zmax(ig)-zmix(ig)).le.0.) then
1965	zmix(ig)=0.9*zmax(ig)
1966	c print*,'pb zmix>zmax'
1967	endif
1968	enddo
1969	do ig=1,klon
1970	zmix0(ig)=zmix(ig)
1971	enddo
1972	c
1973	c calcul du nouveau lmix correspondant
1974	do ig=1,ngrid
1975	do l=1,klev
1976	if (zmix(ig).ge.zlev(ig,l).and.
1977	s zmix(ig).lt.zlev(ig,l+1)) then
1978	lmix(ig)=l
1979	endif
1980	enddo
1981	enddo
1982	c
1983	cne devrait pas arriver!!!!!
1984	do ig=1,ngrid
1985	do l=1,klev
1986	if (detr(ig,l).gt.(fmc(ig,l)+alim(ig,l))+entr(ig,l)) then
1987	print*,'detr2>fmc2!!!','ig=',ig,'l=',l,'d=',detr(ig,l),
1988	s 'f=',fmc(ig,l),'lmax=',lmax(ig)
1989	c detr(ig,l)=fmc(ig,l)+alim(ig,l)+entr(ig,l)
1990	c entr(ig,l)=0.
1991	c fmc(ig,l+1)=0.
1992	c zw2(ig,l+1)=0.
1993	c zqla(ig,l+1)=0.
1994	print*,'pb!fm=0 et f_star>0',l,lmax(ig)
1995	c lmax(ig)=l
1996	endif
1997	enddo
1998	enddo
1999	do ig=1,ngrid
2000	do l=lmax(ig)+1,klev+1
2001	c fmc(ig,l)=0.
2002	c detr(ig,l)=0.
2003	c entr(ig,l)=0.
2004	c zw2(ig,l)=0.
2005	c zqla(ig,l)=0.
2006	enddo
2007	enddo
2008
2009	cCalcul du detrainement lors du premier passage
2010	c print*,'9 OK convect8'
2011	c print*,'WA1 ',wa_moy
2012
2013	c determination de l'indice du debut de la mixed layer ou w decroit
2014
2015	c calcul de la largeur de chaque ascendance dans le cas conservatif.
2016	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
2017	c d'une couche est égale à la hauteur de la couche alimentante.
2018	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
2019	c de la vitesse d'entrainement horizontal dans la couche alimentante.
2020
2021	do l=2,nlay
2022	do ig=1,ngrid
2023	if (l.le.lmax(ig).and.(test(ig).eq.1)) then
2024	zw=max(wa_moy(ig,l),1.e-10)
2025	larg_cons(ig,l)=zmax(ig)*r_aspect
2026	s fmc(ig,l)/(rhobarz(ig,l)zw)
2027	endif
2028	enddo
2029	enddo
2030
2031	do l=2,nlay
2032	do ig=1,ngrid
2033	if (l.le.lmax(ig).and.(test(ig).eq.1)) then
2034	c if (idetr.eq.0) then
2035	c cette option est finalement en dur.
2036	if ((l_mix*zlev(ig,l)).lt.0.)then
2037	print,'pb l_mixzlev<0'
2038	endif
2039	cCR: test: nouvelle def de lambda
2040	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2041	if (zw2(ig,l).gt.1.e-10) then
2042	larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
2043	else
2044	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2045	endif
2046	c else if (idetr.eq.1) then
2047	c larg_detr(ig,l)=larg_cons(ig,l)
2048	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
2049	c else if (idetr.eq.2) then
2050	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2051	c s *sqrt(wa_moy(ig,l))
2052	c else if (idetr.eq.4) then
2053	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2054	c s *wa_moy(ig,l)
2055	c endif
2056	endif
2057	enddo
2058	enddo
2059
2060	c print*,'10 OK convect8'
2061	c print*,'WA2 ',wa_moy
2062	c cal1cul de la fraction de la maille concernée par l'ascendance en tenant
2063	c compte de l'epluchage du thermique.
2064	c
2065	c
2066	do l=2,nlay
2067	do ig=1,ngrid
2068	if(larg_cons(ig,l).gt.1..and.(test(ig).eq.1)) then
2069	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
2070	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
2071	s /(r_aspect*zmax(ig))
2072	c test
2073	fraca(ig,l)=max(fraca(ig,l),0.)
2074	fraca(ig,l)=min(fraca(ig,l),0.5)
2075	fracd(ig,l)=1.-fraca(ig,l)
2076	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
2077	else
2078	c wa_moy(ig,l)=0.
2079	fraca(ig,l)=0.
2080	fracc(ig,l)=0.
2081	fracd(ig,l)=1.
2082	endif
2083	enddo
2084	enddo
2085	cCR: calcul de fracazmix
2086	do ig=1,ngrid
2087	if (test(ig).eq.1) then
2088	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
2089	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
2090	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
2091	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
2092	endif
2093	enddo
2094	c
2095	do l=2,nlay
2096	do ig=1,ngrid
2097	if(larg_cons(ig,l).gt.1..and.(test(ig).eq.1)) then
2098	if (l.gt.lmix(ig)) then
2099	ctest
2100	if (zmax(ig)-zmix(ig).lt.1.e-10) then
2101	c print*,'pb xxx'
2102	xxx(ig,l)=(lmax(ig)+1.-l)/(lmax(ig)+1.-lmix(ig))
2103	else
2104	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
2105	endif
2106	if (idetr.eq.0) then
2107	fraca(ig,l)=fracazmix(ig)
2108	else if (idetr.eq.1) then
2109	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
2110	else if (idetr.eq.2) then
2111	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
2112	else
2113	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
2114	endif
2115	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
2116	fraca(ig,l)=max(fraca(ig,l),0.)
2117	fraca(ig,l)=min(fraca(ig,l),0.5)
2118	fracd(ig,l)=1.-fraca(ig,l)
2119	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
2120	endif
2121	endif
2122	enddo
2123	enddo
2124
2125	print*,'fin calcul fraca'
2126	c print*,'11 OK convect8'
2127	c print*,'Ea3 ',wa_moy
2128	c------------------------------------------------------------------
2129	c Calcul de fracd, wd
2130	c somme wa - wd = 0
2131	c------------------------------------------------------------------
2132
2133
2134	do ig=1,ngrid
2135	fm(ig,1)=0.
2136	fm(ig,nlay+1)=0.
2137	enddo
2138
2139	do l=2,nlay
2140	do ig=1,ngrid
2141	if (test(ig).eq.1) then
2142	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
2143	cCR:test
2144	if (alim(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
2145	s .and.l.gt.lmix(ig)) then
2146	fm(ig,l)=fm(ig,l-1)
2147	c write(1,*)'ajustement fm, l',l
2148	endif
2149	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
2150	cRC
2151	endif
2152	enddo
2153	do ig=1,ngrid
2154	if(fracd(ig,l).lt.0.1.and.(test(ig).eq.1)) then
2155	stop'fracd trop petit'
2156	else
2157	c vitesse descendante "diagnostique"
2158	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
2159	endif
2160	enddo
2161	enddo
2162
2163	do l=1,nlay+1
2164	do ig=1,ngrid
2165	if (test(ig).eq.0) then
2166	fm(ig,l)=fmc(ig,l)
2167	endif
2168	enddo
2169	enddo
2170
2171	cfin du first
2172	do l=1,nlay
2173	do ig=1,ngrid
2174	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
2175	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
2176	enddo
2177	enddo
2178
2179	print*,'12 OK convect8'
2180	c print*,'WA4 ',wa_moy
2181	cc------------------------------------------------------------------
2182	c calcul du transport vertical
2183	c------------------------------------------------------------------
2184
2185	go to 4444
2186	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
2187	do l=2,nlay-1
2188	do ig=1,ngrid
2189	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
2190	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
2191	print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
2192	s ,fm(ig,l+1)*ptimestep
2193	s ,' M=',masse(ig,l),masse(ig,l+1)
2194	endif
2195	enddo
2196	enddo
2197
2198	do l=1,nlay
2199	do ig=1,ngrid
2200	if((alim(ig,l)+entr(ig,l))*ptimestep.gt.masse(ig,l)) then
2201	print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
2202	s ,(entr(ig,l)+alim(ig,l))*ptimestep
2203	s ,' M=',masse(ig,l)
2204	endif
2205	enddo
2206	enddo
2207
2208	do l=1,nlay
2209	do ig=1,ngrid
2210	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
2211	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
2212	c s ,' FM=',fm(ig,l)
2213	endif
2214	if(.not.masse(ig,l).ge.1.e-10
2215	s .or..not.masse(ig,l).le.1.e4) then
2216	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
2217	c s ,' M=',masse(ig,l)
2218	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
2219	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
2220	c print*,'zlev(ig,l+1),zlev(ig,l)'
2221	c s ,zlev(ig,l+1),zlev(ig,l)
2222	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
2223	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
2224	endif
2225	if(.not.alim(ig,l).ge.0..or..not.alim(ig,l).le.10.) then
2226	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
2227	c s ,' E=',entr(ig,l)
2228	endif
2229	enddo
2230	enddo
2231
2232	4444 continue
2233
2234	cCR:redefinition du entr
2235	cCR:test:on ne change pas la def du entr mais la def du fm
2236	do l=1,nlay
2237	do ig=1,ngrid
2238	if (test(ig).eq.1) then
2239	detr(ig,l)=fm(ig,l)+alim(ig,l)-fm(ig,l+1)
2240	if (detr(ig,l).lt.0.) then
2241	c entr(ig,l)=entr(ig,l)-detr(ig,l)
2242	fm(ig,l+1)=fm(ig,l)+alim(ig,l)
2243	detr(ig,l)=0.
2244	c write(11,*)'l,ig,entr',l,ig,entr(ig,l)
2245	c print*,'WARNING !!! detrainement negatif ',ig,l
2246	endif
2247	endif
2248	enddo
2249	enddo
2250	cRC
2251
2252	if (w2di.eq.1) then
2253	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
2254	entr0=entr0+ptimestep*(alim+entr-entr0)/float(tho)
2255	else
2256	fm0=fm
2257	entr0=alim+entr
2258	detr0=detr
2259	alim0=alim
2260	c zoa=zqta
2261	c entr0=alim
2262	endif
2263
2264	if (1.eq.1) then
2265	c call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2266	c . ,zh,zdhadj,zha)
2267	c call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2268	c . ,zo,pdoadj,zoa)
2269	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse,
2270	. zthl,zdthladj,zta)
2271	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse,
2272	. po,pdoadj,zoa)
2273	else
2274	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
2275	. ,zh,zdhadj,zha)
2276	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
2277	. ,zo,pdoadj,zoa)
2278	endif
2279
2280	if (1.eq.0) then
2281	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
2282	. ,fraca,zmax
2283	. ,zu,zv,pduadj,pdvadj,zua,zva)
2284	else
2285	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2286	. ,zu,pduadj,zua)
2287	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2288	. ,zv,pdvadj,zva)
2289	endif
2290
2291	cCalcul des moments
2292	c do l=1,nlay
2293	c do ig=1,ngrid
2294	c zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
2295	c zf2=zf/(1.-zf)
2296	c thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
2297	c wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
2298	c enddo
2299	c enddo
2300
2301
2302
2303
2304
2305
2306	c print*,'13 OK convect8'
2307	c print*,'WA5 ',wa_moy
2308	do l=1,nlay
2309	do ig=1,ngrid
2310	c pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
2311	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
2312	enddo
2313	enddo
2314
2315
2316	c do l=1,nlay
2317	c do ig=1,ngrid
2318	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
2319	c print*,'WARN!!! ig=',ig,' l=',l
2320	c s ,' pdtadj=',pdtadj(ig,l)
2321	c endif
2322	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
2323	c print*,'WARN!!! ig=',ig,' l=',l
2324	c s ,' pdoadj=',pdoadj(ig,l)
2325	c endif
2326	c enddo
2327	c enddo
2328
2329	print*,'14 OK convect8'
2330	c------------------------------------------------------------------
2331	c Calculs pour les sorties
2332	c------------------------------------------------------------------
2333	ccalcul de fraca pour les sorties
2334	do l=2,klev
2335	do ig=1,klon
2336	if (zw2(ig,l).gt.1.e-10) then
2337	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
2338	else
2339	fraca(ig,l)=0.
2340	endif
2341	enddo
2342	enddo
2343	if(sorties) then
2344	do l=1,nlay
2345	do ig=1,ngrid
2346	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
2347	zld(ig,l)=fracd(ig,l)*zmax(ig)
2348	if(1.-fracd(ig,l).gt.1.e-10)
2349	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
2350	enddo
2351	enddo
2352	c CR calcul du niveau de condensation
2353	c initialisation
2354	do ig=1,ngrid
2355	nivcon(ig)=0.
2356	zcon(ig)=0.
2357	enddo
2358	do k=nlay,1,-1
2359	do ig=1,ngrid
2360	if (zqla(ig,k).gt.1e-10) then
2361	nivcon(ig)=k
2362	zcon(ig)=zlev(ig,k)
2363	endif
2364	c if (zcon(ig).gt.1.e-10) then
2365	c nuage=.true.
2366	c else
2367	c nuage=.false.
2368	c endif
2369	enddo
2370	enddo
2371
2372	do l=1,nlay
2373	do ig=1,ngrid
2374	zf=fraca(ig,l)
2375	zf2=zf/(1.-zf)
2376	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))*2
2377	wth2(ig,l)=zf2(zw2(ig,l))*2
2378	c print*,'wth2=',wth2(ig,l)
2379	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l))
2380	s zw2(ig,l)zw2(ig,l)*zw2(ig,l)
2381	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
2382	ctest: on calcul q2/po=ratqsc
2383	c if (nuage) then
2384	ratqscth(ig,l)=sqrt(q2(ig,l))/(po(ig,l)*1000.)
2385	c else
2386	c ratqscth(ig,l)=0.
2387	c endif
2388	enddo
2389	enddo
2390	ccalcul du ratqscdiff
2391	sum=0.
2392	sumdiff=0.
2393	ratqsdiff(:,:)=0.
2394	do ig=1,ngrid
2395	do l=1,lentr(ig)
2396	sum=sum+alim_star(ig,l)zqta(ig,l)1000.
2397	enddo
2398	enddo
2399	do ig=1,ngrid
2400	do l=1,lentr(ig)
2401	zf=fraca(ig,l)
2402	zf2=zf/(1.-zf)
2403	sumdiff=sumdiff+alim_star(ig,l)
2404	s (zqta(ig,l)1000.-sum)**2
2405	c ratqsdiff=ratqsdiff+alim_star(ig,l)*
2406	c s (zqta(ig,l)1000.-po(ig,l)1000.)**2
2407	enddo
2408	enddo
2409	do l=1,klev
2410	do ig=1,ngrid
2411	ratqsdiff(ig,l)=sqrt(sumdiff)/(po(ig,l)*1000.)
2412	c write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
2413	enddo
2414	enddo
2415	cdeja fait
2416	c do l=1,nlay
2417	c do ig=1,ngrid
2418	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
2419	c if (detr(ig,l).lt.0.) then
2420	c entr(ig,l)=entr(ig,l)-detr(ig,l)
2421	c detr(ig,l)=0.
2422	c print*,'WARNING !!! detrainement negatif ',ig,l
2423	c endif
2424	c enddo
2425	c enddo
2426
2427	c print*,'15 OK convect8'
2428
2429	isplit=isplit+1
2430
2431
2432	c #define und
2433	goto 123
2434	#ifdef und
2435	CALL writeg1d(1,nlay,wd,'wd ','wd ')
2436	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
2437	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
2438	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
2439	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
2440	CALL writeg1d(1,nlay,zla,'la ','la ')
2441	CALL writeg1d(1,nlay,zld,'ld ','ld ')
2442	CALL writeg1d(1,nlay,pt,'pt ','pt ')
2443	CALL writeg1d(1,nlay,zh,'zh ','zh ')
2444	CALL writeg1d(1,nlay,zha,'zha ','zha ')
2445	CALL writeg1d(1,nlay,zu,'zu ','zu ')
2446	CALL writeg1d(1,nlay,zv,'zv ','zv ')
2447	CALL writeg1d(1,nlay,zo,'zo ','zo ')
2448	CALL writeg1d(1,nlay,wh,'wh ','wh ')
2449	CALL writeg1d(1,nlay,wu,'wu ','wu ')
2450	CALL writeg1d(1,nlay,wv,'wv ','wv ')
2451	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
2452	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
2453	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
2454	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
2455	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
2456	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
2457	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
2458	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
2459
2460	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
2461	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
2462	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
2463
2464	c recalcul des flux en diagnostique...
2465	c print*,'PAS DE TEMPS ',ptimestep
2466	call dt2F(pplev,pplay,pt,pdtadj,wh)
2467	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
2468	#endif
2469	123 continue
2470	! #define troisD
2471	#ifdef troisD
2472	c if (sorties) then
2473	print*,'Debut des wrgradsfi'
2474
2475	c print*,'16 OK convect8'
2476	c call wrgradsfi(1,nlay,wd,'wd ','wd ')
2477	c call wrgradsfi(1,nlay,zwa,'wa ','wa ')
2478	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
2479	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
2480	c call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
2481	c call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
2482	c print*,'WA6 ',wa_moy
2483	c call wrgradsfi(1,nlay,zla,'la ','la ')
2484	c call wrgradsfi(1,nlay,zld,'ld ','ld ')
2485	call wrgradsfi(1,nlay,pt,'pt ','pt ')
2486	call wrgradsfi(1,nlay,zh,'zh ','zh ')
2487	call wrgradsfi(1,nlay,zha,'zha ','zha ')
2488	call wrgradsfi(1,nlay,zua,'zua ','zua ')
2489	call wrgradsfi(1,nlay,zva,'zva ','zva ')
2490	call wrgradsfi(1,nlay,zu,'zu ','zu ')
2491	call wrgradsfi(1,nlay,zv,'zv ','zv ')
2492	call wrgradsfi(1,nlay,zo,'zo ','zo ')
2493	call wrgradsfi(1,nlay,wh,'wh ','wh ')
2494	call wrgradsfi(1,nlay,wu,'wu ','wu ')
2495	call wrgradsfi(1,nlay,wv,'wv ','wv ')
2496	call wrgradsfi(1,nlay,wo,'wo ','wo ')
2497	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
2498	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
2499	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
2500	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
2501	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
2502	call wrgradsfi(1,nlay,entr,'entr ','entr ')
2503	call wrgradsfi(1,nlay,detr,'detr ','detr ')
2504	call wrgradsfi(1,nlay,fm,'fm ','fm ')
2505	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
2506	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
2507	call wrgradsfi(1,nlay,w_est,'w_est ','w_est ')
2508	con sort les moments
2509	call wrgradsfi(1,nlay,thetath2,'zh2 ','zh2 ')
2510	call wrgradsfi(1,nlay,wth2,'w2 ','w2 ')
2511	call wrgradsfi(1,nlay,wth3,'w3 ','w3 ')
2512	call wrgradsfi(1,nlay,q2,'q2 ','q2 ')
2513	call wrgradsfi(1,nlay,dtheta,'dT ','dT ')
2514	c
2515	call wrgradsfi(1,nlay,zw_sec,'zw_sec ','zw_sec ')
2516	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
2517	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
2518	call wrgradsfi(1,nlay,nu,'nu ','nu ')
2519
2520	call wrgradsfi(1,nlay,zo,'zo ','zo ')
2521	call wrgradsfi(1,nlay,zoa,'zoa ','zoa ')
2522	c call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
2523	c call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
2524
2525	cAM:nouveaux diagnostiques
2526	call wrgradsfi(1,nlay,zthl,'zthl ','zthl ')
2527	call wrgradsfi(1,nlay,zta,'zta ','zta ')
2528	call wrgradsfi(1,nlay,zl,'zl ','zl ')
2529	call wrgradsfi(1,nlay,zdthladj,'zdthladj ',
2530	s 'zdthladj ')
2531	call wrgradsfi(1,nlay,ztla,'ztla ','ztla ')
2532	call wrgradsfi(1,nlay,zqta,'zqta ','zqta ')
2533	call wrgradsfi(1,nlay,zqla,'zqla ','zqla ')
2534	call wrgradsfi(1,nlay,deltaz,'deltaz ','deltaz ')
2535	cCR:nouveaux diagnostiques
2536	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
2537	call wrgradsfi(1,nlay,detr_star ,'detr_star ','detr_star ')
2538	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
2539	call wrgradsfi(1,nlay,zqsat ,'zqsat ','zqsat ')
2540	call wrgradsfi(1,nlay,zqsatth ,'qsath ','qsath ')
2541	call wrgradsfi(1,nlay,alim_star ,'alim_star ','alim_star ')
2542	call wrgradsfi(1,nlay,alim ,'alim ','alim ')
2543	call wrgradsfi(1,1,f,'f ','f ')
2544	call wrgradsfi(1,1,alim_star_tot,'a_s_t ','a_s_t ')
2545	call wrgradsfi(1,1,alim_star2,'a_2 ','a_2 ')
2546	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
2547	call wrgradsfi(1,1,zmax_sec,'z_sec ','z_sec ')
2548	c call wrgradsfi(1,1,zmax_sec2,'zz_se ','zz_se ')
2549	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
2550	call wrgradsfi(1,1,nivcon,'nivcon ','nivcon ')
2551	call wrgradsfi(1,1,zcon,'zcon ','zcon ')
2552	zsortie1d(:)=lmax(:)
2553	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
2554	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
2555	call wrgradsfi(1,1,wmax_sec,'w_sec ','w_sec ')
2556	zsortie1d(:)=lmix(:)
2557	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
2558	zsortie1d(:)=lentr(:)
2559	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
2560
2561	c print*,'17 OK convect8'
2562
2563	do k=1,klev/10
2564	write(str2,'(i2.2)') k
2565	str10='wa'//str2
2566	do l=1,nlay
2567	do ig=1,ngrid
2568	zsortie(ig,l)=wa(ig,k,l)
2569	enddo
2570	enddo
2571	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
2572	do l=1,nlay
2573	do ig=1,ngrid
2574	zsortie(ig,l)=larg_part(ig,k,l)
2575	enddo
2576	enddo
2577	str10='la'//str2
2578	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
2579	enddo
2580
2581
2582	c print*,'18 OK convect8'
2583	c endif
2584	print*,'Fin des wrgradsfi'
2585	#endif
2586
2587	endif
2588
2589	c if(wa_moy(1,4).gt.1.e-10) stop
2590
2591	print*,'19 OK convect8'
2592	return
2593	end
2594	SUBROUTINE thermcell_eau(ngrid,nlay,ptimestep
2595	s ,pplay,pplev,pphi
2596	s ,pu,pv,pt,po
2597	s ,pduadj,pdvadj,pdtadj,pdoadj
2598	s ,fm0,entr0
2599	c s ,pu_therm,pv_therm
2600	s ,r_aspect,l_mix,w2di,tho)
2601
2602	USE dimphy
2603	IMPLICIT NONE
2604
2605	c=======================================================================
2606	c
2607	c Calcul du transport verticale dans la couche limite en presence
2608	c de "thermiques" explicitement representes
2609	c
2610	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
2611	c
2612	c le thermique est supposé homogène et dissipé par mélange avec
2613	c son environnement. la longueur l_mix contrôle l'efficacité du
2614	c mélange
2615	c
2616	c Le calcul du transport des différentes espèces se fait en prenant
2617	c en compte:
2618	c 1. un flux de masse montant
2619	c 2. un flux de masse descendant
2620	c 3. un entrainement
2621	c 4. un detrainement
2622	c
2623	c=======================================================================
2624
2625	c-----------------------------------------------------------------------
2626	c declarations:
2627	c -------------
2628
2629	#include "dimensions.h"
2630	cccc#include "dimphy.h"
2631	#include "YOMCST.h"
2632	#include "YOETHF.h"
2633	#include "FCTTRE.h"
2634
2635	c arguments:
2636	c ----------
2637
2638	INTEGER ngrid,nlay,w2di,tho
2639	real ptimestep,l_mix,r_aspect
2640	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
2641	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
2642	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
2643	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
2644	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
2645	real pphi(ngrid,nlay)
2646
2647	integer idetr
2648	save idetr
2649	data idetr/3/
2650
2651	c local:
2652	c ------
2653
2654	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
2655	real zsortie1d(klon)
2656	c CR: on remplace lmax(klon,klev+1)
2657	INTEGER lmax(klon),lmin(klon),lentr(klon)
2658	real linter(klon)
2659	real zmix(klon), fracazmix(klon)
2660	c RC
2661	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
2662
2663	real zlev(klon,klev+1),zlay(klon,klev)
2664	REAL zh(klon,klev),zdhadj(klon,klev)
2665	real zthl(klon,klev),zdthladj(klon,klev)
2666	REAL ztv(klon,klev)
2667	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
2668	real zl(klon,klev)
2669	REAL wh(klon,klev+1)
2670	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
2671	real zla(klon,klev+1)
2672	real zwa(klon,klev+1)
2673	real zld(klon,klev+1)
2674	real zwd(klon,klev+1)
2675	real zsortie(klon,klev)
2676	real zva(klon,klev)
2677	real zua(klon,klev)
2678	real zoa(klon,klev)
2679
2680	real zta(klon,klev)
2681	real zha(klon,klev)
2682	real wa_moy(klon,klev+1)
2683	real fraca(klon,klev+1)
2684	real fracc(klon,klev+1)
2685	real zf,zf2
2686	real thetath2(klon,klev),wth2(klon,klev)
2687	! common/comtherm/thetath2,wth2
2688
2689	real count_time
2690	integer isplit,nsplit,ialt
2691	parameter (nsplit=10)
2692	data isplit/0/
2693	save isplit
2694
2695	logical sorties
2696	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
2697	real zpspsk(klon,klev)
2698
2699	c real wmax(klon,klev),wmaxa(klon)
2700	real wmax(klon),wmaxa(klon)
2701	real wa(klon,klev,klev+1)
2702	real wd(klon,klev+1)
2703	real larg_part(klon,klev,klev+1)
2704	real fracd(klon,klev+1)
2705	real xxx(klon,klev+1)
2706	real larg_cons(klon,klev+1)
2707	real larg_detr(klon,klev+1)
2708	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
2709	real pu_therm(klon,klev),pv_therm(klon,klev)
2710	real fm(klon,klev+1),entr(klon,klev)
2711	real fmc(klon,klev+1)
2712
2713	real zcor,zdelta,zcvm5,qlbef
2714	real Tbef(klon),qsatbef(klon)
2715	real dqsat_dT,DT,num,denom
2716	REAL REPS,RLvCp,DDT0
2717	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
2718
2719	PARAMETER (DDT0=.01)
2720
2721	cCR:nouvelles variables
2722	real f_star(klon,klev+1),entr_star(klon,klev)
2723	real entr_star_tot(klon),entr_star2(klon)
2724	real f(klon), f0(klon)
2725	real zlevinter(klon)
2726	logical first
2727	data first /.false./
2728	save first
2729	cRC
2730
2731	character*2 str2
2732	character*10 str10
2733
2734	LOGICAL vtest(klon),down
2735	LOGICAL Zsat(klon)
2736
2737	EXTERNAL SCOPY
2738
2739	integer ncorrec,ll
2740	save ncorrec
2741	data ncorrec/0/
2742	c
2743
2744	c-----------------------------------------------------------------------
2745	c initialisation:
2746	c ---------------
2747	c
2748	sorties=.true.
2749	IF(ngrid.NE.klon) THEN
2750	PRINT*
2751	PRINT*,'STOP dans convadj'
2752	PRINT*,'ngrid =',ngrid
2753	PRINT*,'klon =',klon
2754	ENDIF
2755	c
2756	c Initialisation
2757	RLvCp = RLVTT/RCPD
2758	REPS = RD/RV
2759	c
2760	c-----------------------------------------------------------------------
2761	cAM Calcul de T,q,ql a partir de Tl et qT
2762	c ---------------------------------------------------
2763	c
2764	c Pr Tprec=Tl calcul de qsat
2765	c Si qsat>qT T=Tl, q=qT
2766	c Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
2767	c On cherche DDT < DDT0
2768	c
2769	c defaut
2770	DO ll=1,nlay
2771	DO ig=1,ngrid
2772	zo(ig,ll)=po(ig,ll)
2773	zl(ig,ll)=0.
2774	zh(ig,ll)=pt(ig,ll)
2775	EndDO
2776	EndDO
2777	do ig=1,ngrid
2778	Zsat(ig)=.false.
2779	enddo
2780	c
2781	c
2782	DO ll=1,nlay
2783	c les points insatures sont definitifs
2784	DO ig=1,ngrid
2785	Tbef(ig)=pt(ig,ll)
2786	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
2787	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
2788	qsatbef(ig)=MIN(0.5,qsatbef(ig))
2789	zcor=1./(1.-retv*qsatbef(ig))
2790	qsatbef(ig)=qsatbef(ig)*zcor
2791	Zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig)) .gt. 0.00001)
2792	EndDO
2793
2794	DO ig=1,ngrid
2795	if (Zsat(ig)) then
2796	qlbef=max(0.,po(ig,ll)-qsatbef(ig))
2797	c si sature: ql est surestime, d'ou la sous-relax
2798	DT = 0.5RLvCpqlbef
2799	c on pourra enchainer 2 ou 3 calculs sans Do while
2800	do while (DT.gt.DDT0)
2801	c il faut verifier si c,a conserve quand on repasse en insature ...
2802	Tbef(ig)=Tbef(ig)+DT
2803	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
2804	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,ll)
2805	qsatbef(ig)=MIN(0.5,qsatbef(ig))
2806	zcor=1./(1.-retv*qsatbef(ig))
2807	qsatbef(ig)=qsatbef(ig)*zcor
2808	c on veut le signe de qlbef
2809	qlbef=po(ig,ll)-qsatbef(ig)
2810	c dqsat_dT
2811	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
2812	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
2813	zcor=1./(1.-retv*qsatbef(ig))
2814	dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
2815	num=-Tbef(ig)+pt(ig,ll)+RLvCp*qlbef
2816	denom=1.+RLvCp*dqsat_dT
2817	DT=num/denom
2818	enddo
2819	c on ecrit de maniere conservative (sat ou non)
2820	zl(ig,ll) = max(0.,qlbef)
2821	c T = Tl +Lv/Cp ql
2822	zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
2823	zo(ig,ll) = po(ig,ll)-zl(ig,ll)
2824	endif
2825	EndDO
2826	EndDO
2827	cAM fin
2828	c
2829	c-----------------------------------------------------------------------
2830	c incrementation eventuelle de tendances precedentes:
2831	c ---------------------------------------------------
2832
2833	print*,'0 OK convect8'
2834
2835	DO 1010 l=1,nlay
2836	DO 1015 ig=1,ngrid
2837	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
2838	c zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
2839	zu(ig,l)=pu(ig,l)
2840	zv(ig,l)=pv(ig,l)
2841	c zo(ig,l)=po(ig,l)
2842	c ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
2843	cAM attention zh est maintenant le profil de T et plus le profil de theta !
2844	c
2845	c T-> Theta
2846	ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
2847	cAM Theta_v
2848	ztv(ig,l)=ztv(ig,l)(1.+RETV(zo(ig,l))
2849	s -zl(ig,l))
2850	cAM Thetal
2851	zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
2852	c
2853	1015 CONTINUE
2854	1010 CONTINUE
2855
2856	c print*,'1 OK convect8'
2857	c --------------------
2858	c
2859	c
2860	c + + + + + + + + + + +
2861	c
2862	c
2863	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
2864	c wh,wt,wo ...
2865	c
2866	c + + + + + + + + + + + zh,zu,zv,zo,rho
2867	c
2868	c
2869	c -------------------- zlev(1)
2870	c \\\\\\\\\\\\\\\\\\\\
2871	c
2872	c
2873
2874	c-----------------------------------------------------------------------
2875	c Calcul des altitudes des couches
2876	c-----------------------------------------------------------------------
2877
2878	do l=2,nlay
2879	do ig=1,ngrid
2880	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
2881	enddo
2882	enddo
2883	do ig=1,ngrid
2884	zlev(ig,1)=0.
2885	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
2886	enddo
2887	do l=1,nlay
2888	do ig=1,ngrid
2889	zlay(ig,l)=pphi(ig,l)/RG
2890	enddo
2891	enddo
2892
2893	c print*,'2 OK convect8'
2894	c-----------------------------------------------------------------------
2895	c Calcul des densites
2896	c-----------------------------------------------------------------------
2897
2898	do l=1,nlay
2899	do ig=1,ngrid
2900	c rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
2901	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDztv(ig,l))
2902	enddo
2903	enddo
2904
2905	do l=2,nlay
2906	do ig=1,ngrid
2907	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
2908	enddo
2909	enddo
2910
2911	do k=1,nlay
2912	do l=1,nlay+1
2913	do ig=1,ngrid
2914	wa(ig,k,l)=0.
2915	enddo
2916	enddo
2917	enddo
2918
2919	c print*,'3 OK convect8'
2920	c------------------------------------------------------------------
2921	c Calcul de w2, quarre de w a partir de la cape
2922	c a partir de w2, on calcule wa, vitesse de l'ascendance
2923	c
2924	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
2925	c w2 est stoke dans wa
2926	c
2927	c ATTENTION: dans convect8, on n'utilise le calcule des wa
2928	c independants par couches que pour calculer l'entrainement
2929	c a la base et la hauteur max de l'ascendance.
2930	c
2931	c Indicages:
2932	c l'ascendance provenant du niveau k traverse l'interface l avec
2933	c une vitesse wa(k,l).
2934	c
2935	c --------------------
2936	c
2937	c + + + + + + + + + +
2938	c
2939	c wa(k,l) ---- -------------------- l
2940	c /\
2941	c /\|\|\ + + + + + + + + + +
2942	c \|\|
2943	c \|\| --------------------
2944	c \|\|
2945	c \|\| + + + + + + + + + +
2946	c \|\|
2947	c \|\| --------------------
2948	c \|\|__
2949	c \|___ + + + + + + + + + + k
2950	c
2951	c --------------------
2952	c
2953	c
2954	c
2955	c------------------------------------------------------------------
2956
2957	cCR: ponderation entrainement des couches instables
2958	cdef des entr_star tels que entr=f*entr_star
2959	do l=1,klev
2960	do ig=1,ngrid
2961	entr_star(ig,l)=0.
2962	enddo
2963	enddo
2964	c determination de la longueur de la couche d entrainement
2965	do ig=1,ngrid
2966	lentr(ig)=1
2967	enddo
2968
2969	con ne considere que les premieres couches instables
2970	do k=nlay-1,1,-1
2971	do ig=1,ngrid
2972	if (ztv(ig,k).gt.ztv(ig,k+1).and.
2973	s ztv(ig,k+1).lt.ztv(ig,k+2)) then
2974	lentr(ig)=k
2975	endif
2976	enddo
2977	enddo
2978
2979	c determination du lmin: couche d ou provient le thermique
2980	do ig=1,ngrid
2981	lmin(ig)=1
2982	enddo
2983	do ig=1,ngrid
2984	do l=nlay,2,-1
2985	if (ztv(ig,l-1).gt.ztv(ig,l)) then
2986	lmin(ig)=l-1
2987	endif
2988	enddo
2989	enddo
2990	c
2991	c definition de l'entrainement des couches
2992	do l=1,klev-1
2993	do ig=1,ngrid
2994	if (ztv(ig,l).gt.ztv(ig,l+1).and.
2995	s l.ge.lmin(ig).and.l.le.lentr(ig)) then
2996	entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
2997	s (zlev(ig,l+1)-zlev(ig,l))
2998	endif
2999	enddo
3000	enddo
3001	c pas de thermique si couche 1 stable
3002	do ig=1,ngrid
3003	if (lmin(ig).gt.1) then
3004	do l=1,klev
3005	entr_star(ig,l)=0.
3006	enddo
3007	endif
3008	enddo
3009	c calcul de l entrainement total
3010	do ig=1,ngrid
3011	entr_star_tot(ig)=0.
3012	enddo
3013	do ig=1,ngrid
3014	do k=1,klev
3015	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
3016	enddo
3017	enddo
3018	c
3019	do k=1,klev
3020	do ig=1,ngrid
3021	ztva(ig,k)=ztv(ig,k)
3022	enddo
3023	enddo
3024	cRC
3025	cAM:initialisations
3026	do k=1,nlay
3027	do ig=1,ngrid
3028	ztva(ig,k)=ztv(ig,k)
3029	ztla(ig,k)=zthl(ig,k)
3030	zqla(ig,k)=0.
3031	zqta(ig,k)=po(ig,k)
3032	Zsat(ig) =.false.
3033	enddo
3034	enddo
3035	c
3036	c print*,'7 OK convect8'
3037	do k=1,klev+1
3038	do ig=1,ngrid
3039	zw2(ig,k)=0.
3040	fmc(ig,k)=0.
3041	cCR
3042	f_star(ig,k)=0.
3043	cRC
3044	larg_cons(ig,k)=0.
3045	larg_detr(ig,k)=0.
3046	wa_moy(ig,k)=0.
3047	enddo
3048	enddo
3049
3050	c print*,'8 OK convect8'
3051	do ig=1,ngrid
3052	linter(ig)=1.
3053	lmaxa(ig)=1
3054	lmix(ig)=1
3055	wmaxa(ig)=0.
3056	enddo
3057
3058	cCR:
3059	do l=1,nlay-2
3060	do ig=1,ngrid
3061	if (ztv(ig,l).gt.ztv(ig,l+1)
3062	s .and.entr_star(ig,l).gt.1.e-10
3063	s .and.zw2(ig,l).lt.1e-10) then
3064	cAM
3065	ztla(ig,l)=zthl(ig,l)
3066	zqta(ig,l)=po(ig,l)
3067	zqla(ig,l)=zl(ig,l)
3068	cAM
3069	f_star(ig,l+1)=entr_star(ig,l)
3070	ctest:calcul de dteta
3071	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
3072	s *(zlev(ig,l+1)-zlev(ig,l))
3073	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
3074	larg_detr(ig,l)=0.
3075	else if ((zw2(ig,l).ge.1e-10).and.
3076	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
3077	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
3078	c
3079	cAM on melange Tl et qt du thermique
3080	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+entr_star(ig,l)
3081	s *zthl(ig,l))/f_star(ig,l+1)
3082	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+entr_star(ig,l)
3083	s *po(ig,l))/f_star(ig,l+1)
3084	c
3085	c ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
3086	c s *ztv(ig,l))/f_star(ig,l+1)
3087	c
3088	cAM on en deduit thetav et ql du thermique
3089	Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
3090	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
3091	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
3092	qsatbef(ig)=MIN(0.5,qsatbef(ig))
3093	zcor=1./(1.-retv*qsatbef(ig))
3094	qsatbef(ig)=qsatbef(ig)*zcor
3095	Zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig)) .gt. 0.00001)
3096	endif
3097	enddo
3098	DO ig=1,ngrid
3099	if (Zsat(ig)) then
3100	qlbef=max(0.,zqta(ig,l)-qsatbef(ig))
3101	DT = 0.5RLvCpqlbef
3102	do while (DT.gt.DDT0)
3103	Tbef(ig)=Tbef(ig)+DT
3104	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
3105	qsatbef(ig)= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig,l)
3106	qsatbef(ig)=MIN(0.5,qsatbef(ig))
3107	zcor=1./(1.-retv*qsatbef(ig))
3108	qsatbef(ig)=qsatbef(ig)*zcor
3109	qlbef=zqta(ig,l)-qsatbef(ig)
3110
3111	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
3112	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
3113	zcor=1./(1.-retv*qsatbef(ig))
3114	dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef(ig),zcor)
3115	num=-Tbef(ig)+ztla(ig,l)zpspsk(ig,l)+RLvCpqlbef
3116	denom=1.+RLvCp*dqsat_dT
3117	DT=num/denom
3118	enddo
3119	zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
3120	endif
3121	c on ecrit de maniere conservative (sat ou non)
3122	c T = Tl +Lv/Cp ql
3123	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
3124	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
3125	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)
3126	s -zqla(ig,l))-zqla(ig,l))
3127
3128	enddo
3129	DO ig=1,ngrid
3130	if (zw2(ig,l).ge.1.e-10.and.
3131	s f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
3132	c mise a jour de la vitesse ascendante (l'air entraine de la couche
3133	c consideree commence avec une vitesse nulle).
3134	c
3135	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
3136	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
3137	s *(zlev(ig,l+1)-zlev(ig,l))
3138	endif
3139	c determination de zmax continu par interpolation lineaire
3140	if (zw2(ig,l+1).lt.0.) then
3141	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
3142	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
3143	zw2(ig,l+1)=0.
3144	lmaxa(ig)=l
3145	else
3146	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
3147	endif
3148	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
3149	c lmix est le niveau de la couche ou w (wa_moy) est maximum
3150	lmix(ig)=l+1
3151	wmaxa(ig)=wa_moy(ig,l+1)
3152	endif
3153	enddo
3154	enddo
3155	c
3156	c Calcul de la couche correspondant a la hauteur du thermique
3157	do ig=1,ngrid
3158	lmax(ig)=lentr(ig)
3159	enddo
3160	do ig=1,ngrid
3161	do l=nlay,lentr(ig)+1,-1
3162	if (zw2(ig,l).le.1.e-10) then
3163	lmax(ig)=l-1
3164	endif
3165	enddo
3166	enddo
3167	c pas de thermique si couche 1 stable
3168	do ig=1,ngrid
3169	if (lmin(ig).gt.1) then
3170	lmax(ig)=1
3171	lmin(ig)=1
3172	endif
3173	enddo
3174	c
3175	c Determination de zw2 max
3176	do ig=1,ngrid
3177	wmax(ig)=0.
3178	enddo
3179
3180	do l=1,nlay
3181	do ig=1,ngrid
3182	if (l.le.lmax(ig)) then
3183	zw2(ig,l)=sqrt(zw2(ig,l))
3184	wmax(ig)=max(wmax(ig),zw2(ig,l))
3185	else
3186	zw2(ig,l)=0.
3187	endif
3188	enddo
3189	enddo
3190
3191	c Longueur caracteristique correspondant a la hauteur des thermiques.
3192	do ig=1,ngrid
3193	zmax(ig)=500.
3194	zlevinter(ig)=zlev(ig,1)
3195	enddo
3196	do ig=1,ngrid
3197	c calcul de zlevinter
3198	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
3199	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
3200	s -zlev(ig,lmax(ig)))
3201	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
3202	enddo
3203
3204	c Fermeture,determination de f
3205	do ig=1,ngrid
3206	entr_star2(ig)=0.
3207	enddo
3208	do ig=1,ngrid
3209	if (entr_star_tot(ig).LT.1.e-10) then
3210	f(ig)=0.
3211	else
3212	do k=lmin(ig),lentr(ig)
3213	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
3214	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
3215	enddo
3216	c Nouvelle fermeture
3217	f(ig)=wmax(ig)/(zmax(ig)r_aspectentr_star2(ig))
3218	s *entr_star_tot(ig)
3219	ctest
3220	if (first) then
3221	f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
3222	s *wmax(ig))
3223	endif
3224	endif
3225	f0(ig)=f(ig)
3226	first=.true.
3227	enddo
3228
3229	c Calcul de l'entrainement
3230	do k=1,klev
3231	do ig=1,ngrid
3232	entr(ig,k)=f(ig)*entr_star(ig,k)
3233	enddo
3234	enddo
3235	c Calcul des flux
3236	do ig=1,ngrid
3237	do l=1,lmax(ig)-1
3238	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
3239	enddo
3240	enddo
3241
3242	cRC
3243
3244
3245	c print*,'9 OK convect8'
3246	c print*,'WA1 ',wa_moy
3247
3248	c determination de l'indice du debut de la mixed layer ou w decroit
3249
3250	c calcul de la largeur de chaque ascendance dans le cas conservatif.
3251	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
3252	c d'une couche est égale à la hauteur de la couche alimentante.
3253	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
3254	c de la vitesse d'entrainement horizontal dans la couche alimentante.
3255
3256	do l=2,nlay
3257	do ig=1,ngrid
3258	if (l.le.lmaxa(ig)) then
3259	zw=max(wa_moy(ig,l),1.e-10)
3260	larg_cons(ig,l)=zmax(ig)*r_aspect
3261	s fmc(ig,l)/(rhobarz(ig,l)zw)
3262	endif
3263	enddo
3264	enddo
3265
3266	do l=2,nlay
3267	do ig=1,ngrid
3268	if (l.le.lmaxa(ig)) then
3269	c if (idetr.eq.0) then
3270	c cette option est finalement en dur.
3271	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3272	c else if (idetr.eq.1) then
3273	c larg_detr(ig,l)=larg_cons(ig,l)
3274	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
3275	c else if (idetr.eq.2) then
3276	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3277	c s *sqrt(wa_moy(ig,l))
3278	c else if (idetr.eq.4) then
3279	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3280	c s *wa_moy(ig,l)
3281	c endif
3282	endif
3283	enddo
3284	enddo
3285
3286	c print*,'10 OK convect8'
3287	c print*,'WA2 ',wa_moy
3288	c calcul de la fraction de la maille concernée par l'ascendance en tenant
3289	c compte de l'epluchage du thermique.
3290	c
3291	cCR def de zmix continu (profil parabolique des vitesses)
3292	do ig=1,ngrid
3293	if (lmix(ig).gt.1.) then
3294	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
3295	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
3296	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
3297	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
3298	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
3299	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
3300	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
3301	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3302	else
3303	zmix(ig)=0.
3304	endif
3305	enddo
3306	c
3307	c calcul du nouveau lmix correspondant
3308	do ig=1,ngrid
3309	do l=1,klev
3310	if (zmix(ig).ge.zlev(ig,l).and.
3311	s zmix(ig).lt.zlev(ig,l+1)) then
3312	lmix(ig)=l
3313	endif
3314	enddo
3315	enddo
3316	c
3317	do l=2,nlay
3318	do ig=1,ngrid
3319	if(larg_cons(ig,l).gt.1.) then
3320	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
3321	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
3322	s /(r_aspect*zmax(ig))
3323	c test
3324	fraca(ig,l)=max(fraca(ig,l),0.)
3325	fraca(ig,l)=min(fraca(ig,l),0.5)
3326	fracd(ig,l)=1.-fraca(ig,l)
3327	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
3328	else
3329	c wa_moy(ig,l)=0.
3330	fraca(ig,l)=0.
3331	fracc(ig,l)=0.
3332	fracd(ig,l)=1.
3333	endif
3334	enddo
3335	enddo
3336	cCR: calcul de fracazmix
3337	do ig=1,ngrid
3338	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
3339	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
3340	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
3341	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3342	enddo
3343	c
3344	do l=2,nlay
3345	do ig=1,ngrid
3346	if(larg_cons(ig,l).gt.1.) then
3347	if (l.gt.lmix(ig)) then
3348	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3349	if (idetr.eq.0) then
3350	fraca(ig,l)=fracazmix(ig)
3351	else if (idetr.eq.1) then
3352	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
3353	else if (idetr.eq.2) then
3354	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
3355	else
3356	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
3357	endif
3358	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3359	fraca(ig,l)=max(fraca(ig,l),0.)
3360	fraca(ig,l)=min(fraca(ig,l),0.5)
3361	fracd(ig,l)=1.-fraca(ig,l)
3362	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
3363	endif
3364	endif
3365	enddo
3366	enddo
3367
3368	c print*,'11 OK convect8'
3369	c print*,'Ea3 ',wa_moy
3370	c------------------------------------------------------------------
3371	c Calcul de fracd, wd
3372	c somme wa - wd = 0
3373	c------------------------------------------------------------------
3374
3375
3376	do ig=1,ngrid
3377	fm(ig,1)=0.
3378	fm(ig,nlay+1)=0.
3379	enddo
3380
3381	do l=2,nlay
3382	do ig=1,ngrid
3383	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
3384	cCR:test
3385	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
3386	s .and.l.gt.lmix(ig)) then
3387	fm(ig,l)=fm(ig,l-1)
3388	c write(1,*)'ajustement fm, l',l
3389	endif
3390	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3391	cRC
3392	enddo
3393	do ig=1,ngrid
3394	if(fracd(ig,l).lt.0.1) then
3395	stop'fracd trop petit'
3396	else
3397	c vitesse descendante "diagnostique"
3398	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
3399	endif
3400	enddo
3401	enddo
3402
3403	do l=1,nlay
3404	do ig=1,ngrid
3405	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3406	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
3407	enddo
3408	enddo
3409
3410	c print*,'12 OK convect8'
3411	c print*,'WA4 ',wa_moy
3412	cc------------------------------------------------------------------
3413	c calcul du transport vertical
3414	c------------------------------------------------------------------
3415
3416	go to 4444
3417	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3418	do l=2,nlay-1
3419	do ig=1,ngrid
3420	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
3421	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
3422	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
3423	c s ,fm(ig,l+1)*ptimestep
3424	c s ,' M=',masse(ig,l),masse(ig,l+1)
3425	endif
3426	enddo
3427	enddo
3428
3429	do l=1,nlay
3430	do ig=1,ngrid
3431	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
3432	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
3433	c s ,entr(ig,l)*ptimestep
3434	c s ,' M=',masse(ig,l)
3435	endif
3436	enddo
3437	enddo
3438
3439	do l=1,nlay
3440	do ig=1,ngrid
3441	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
3442	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
3443	c s ,' FM=',fm(ig,l)
3444	endif
3445	if(.not.masse(ig,l).ge.1.e-10
3446	s .or..not.masse(ig,l).le.1.e4) then
3447	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
3448	c s ,' M=',masse(ig,l)
3449	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3450	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3451	c print*,'zlev(ig,l+1),zlev(ig,l)'
3452	c s ,zlev(ig,l+1),zlev(ig,l)
3453	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3454	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3455	endif
3456	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
3457	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
3458	c s ,' E=',entr(ig,l)
3459	endif
3460	enddo
3461	enddo
3462
3463	4444 continue
3464
3465	if (w2di.eq.1) then
3466	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
3467	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
3468	else
3469	fm0=fm
3470	entr0=entr
3471	endif
3472
3473	if (1.eq.1) then
3474	c call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3475	c . ,zh,zdhadj,zha)
3476	c call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3477	c . ,zo,pdoadj,zoa)
3478	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3479	. ,zthl,zdthladj,zta)
3480	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3481	. ,po,pdoadj,zoa)
3482	else
3483	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
3484	. ,zh,zdhadj,zha)
3485	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
3486	. ,zo,pdoadj,zoa)
3487	endif
3488
3489	if (1.eq.0) then
3490	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
3491	. ,fraca,zmax
3492	. ,zu,zv,pduadj,pdvadj,zua,zva)
3493	else
3494	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3495	. ,zu,pduadj,zua)
3496	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3497	. ,zv,pdvadj,zva)
3498	endif
3499
3500	do l=1,nlay
3501	do ig=1,ngrid
3502	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
3503	zf2=zf/(1.-zf)
3504	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
3505	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
3506	enddo
3507	enddo
3508
3509
3510
3511	c print*,'13 OK convect8'
3512	c print*,'WA5 ',wa_moy
3513	do l=1,nlay
3514	do ig=1,ngrid
3515	c pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
3516	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
3517	enddo
3518	enddo
3519
3520
3521	c do l=1,nlay
3522	c do ig=1,ngrid
3523	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
3524	c print*,'WARN!!! ig=',ig,' l=',l
3525	c s ,' pdtadj=',pdtadj(ig,l)
3526	c endif
3527	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
3528	c print*,'WARN!!! ig=',ig,' l=',l
3529	c s ,' pdoadj=',pdoadj(ig,l)
3530	c endif
3531	c enddo
3532	c enddo
3533
3534	c print*,'14 OK convect8'
3535	c------------------------------------------------------------------
3536	c Calculs pour les sorties
3537	c------------------------------------------------------------------
3538
3539	if(sorties) then
3540	do l=1,nlay
3541	do ig=1,ngrid
3542	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
3543	zld(ig,l)=fracd(ig,l)*zmax(ig)
3544	if(1.-fracd(ig,l).gt.1.e-10)
3545	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
3546	enddo
3547	enddo
3548
3549	do l=1,nlay
3550	do ig=1,ngrid
3551	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
3552	if (detr(ig,l).lt.0.) then
3553	entr(ig,l)=entr(ig,l)-detr(ig,l)
3554	detr(ig,l)=0.
3555	c print*,'WARNING !!! detrainement negatif ',ig,l
3556	endif
3557	enddo
3558	enddo
3559
3560	c print*,'15 OK convect8'
3561
3562	isplit=isplit+1
3563
3564
3565	c #define und
3566	goto 123
3567	#ifdef und
3568	CALL writeg1d(1,nlay,wd,'wd ','wd ')
3569	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
3570	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
3571	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
3572	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
3573	CALL writeg1d(1,nlay,zla,'la ','la ')
3574	CALL writeg1d(1,nlay,zld,'ld ','ld ')
3575	CALL writeg1d(1,nlay,pt,'pt ','pt ')
3576	CALL writeg1d(1,nlay,zh,'zh ','zh ')
3577	CALL writeg1d(1,nlay,zha,'zha ','zha ')
3578	CALL writeg1d(1,nlay,zu,'zu ','zu ')
3579	CALL writeg1d(1,nlay,zv,'zv ','zv ')
3580	CALL writeg1d(1,nlay,zo,'zo ','zo ')
3581	CALL writeg1d(1,nlay,wh,'wh ','wh ')
3582	CALL writeg1d(1,nlay,wu,'wu ','wu ')
3583	CALL writeg1d(1,nlay,wv,'wv ','wv ')
3584	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
3585	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
3586	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
3587	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
3588	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
3589	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
3590	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
3591	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
3592
3593	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
3594	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
3595	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
3596
3597	c recalcul des flux en diagnostique...
3598	c print*,'PAS DE TEMPS ',ptimestep
3599	call dt2F(pplev,pplay,pt,pdtadj,wh)
3600	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
3601	#endif
3602	123 continue
3603	! #define troisD
3604	#ifdef troisD
3605	c if (sorties) then
3606	print*,'Debut des wrgradsfi'
3607
3608	c print*,'16 OK convect8'
3609	call wrgradsfi(1,nlay,wd,'wd ','wd ')
3610	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
3611	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
3612	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
3613	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
3614	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
3615	c print*,'WA6 ',wa_moy
3616	call wrgradsfi(1,nlay,zla,'la ','la ')
3617	call wrgradsfi(1,nlay,zld,'ld ','ld ')
3618	call wrgradsfi(1,nlay,pt,'pt ','pt ')
3619	call wrgradsfi(1,nlay,zh,'zh ','zh ')
3620	call wrgradsfi(1,nlay,zha,'zha ','zha ')
3621	call wrgradsfi(1,nlay,zua,'zua ','zua ')
3622	call wrgradsfi(1,nlay,zva,'zva ','zva ')
3623	call wrgradsfi(1,nlay,zu,'zu ','zu ')
3624	call wrgradsfi(1,nlay,zv,'zv ','zv ')
3625	call wrgradsfi(1,nlay,zo,'zo ','zo ')
3626	call wrgradsfi(1,nlay,wh,'wh ','wh ')
3627	call wrgradsfi(1,nlay,wu,'wu ','wu ')
3628	call wrgradsfi(1,nlay,wv,'wv ','wv ')
3629	call wrgradsfi(1,nlay,wo,'wo ','wo ')
3630	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
3631	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
3632	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
3633	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
3634	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
3635	call wrgradsfi(1,nlay,entr,'entr ','entr ')
3636	call wrgradsfi(1,nlay,detr,'detr ','detr ')
3637	call wrgradsfi(1,nlay,fm,'fm ','fm ')
3638	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
3639	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
3640	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
3641	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
3642
3643	call wrgradsfi(1,nlay,zo,'zo ','zo ')
3644	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
3645	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
3646
3647	cAM:nouveaux diagnostiques
3648	call wrgradsfi(1,nlay,zthl,'zthl ','zthl ')
3649	call wrgradsfi(1,nlay,zta,'zta ','zta ')
3650	call wrgradsfi(1,nlay,zl,'zl ','zl ')
3651	call wrgradsfi(1,nlay,zdthladj,'zdthladj ',
3652	s 'zdthladj ')
3653	call wrgradsfi(1,nlay,ztla,'ztla ','ztla ')
3654	call wrgradsfi(1,nlay,zqta,'zqta ','zqta ')
3655	call wrgradsfi(1,nlay,zqla,'zqla ','zqla ')
3656	cCR:nouveaux diagnostiques
3657	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
3658	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
3659	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
3660	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
3661	zsortie1d(:)=lmax(:)
3662	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
3663	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
3664	zsortie1d(:)=lmix(:)
3665	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
3666	zsortie1d(:)=lentr(:)
3667	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
3668
3669	c print*,'17 OK convect8'
3670
3671	do k=1,klev/10
3672	write(str2,'(i2.2)') k
3673	str10='wa'//str2
3674	do l=1,nlay
3675	do ig=1,ngrid
3676	zsortie(ig,l)=wa(ig,k,l)
3677	enddo
3678	enddo
3679	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
3680	do l=1,nlay
3681	do ig=1,ngrid
3682	zsortie(ig,l)=larg_part(ig,k,l)
3683	enddo
3684	enddo
3685	str10='la'//str2
3686	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
3687	enddo
3688
3689
3690	c print*,'18 OK convect8'
3691	c endif
3692	print*,'Fin des wrgradsfi'
3693	#endif
3694
3695	endif
3696
3697	c if(wa_moy(1,4).gt.1.e-10) stop
3698
3699	c print*,'19 OK convect8'
3700	return
3701	end
3702
3703	SUBROUTINE thermcell(ngrid,nlay,ptimestep
3704	s ,pplay,pplev,pphi
3705	s ,pu,pv,pt,po
3706	s ,pduadj,pdvadj,pdtadj,pdoadj
3707	s ,fm0,entr0
3708	c s ,pu_therm,pv_therm
3709	s ,r_aspect,l_mix,w2di,tho)
3710
3711	USE dimphy
3712	IMPLICIT NONE
3713
3714	c=======================================================================
3715	c
3716	c Calcul du transport verticale dans la couche limite en presence
3717	c de "thermiques" explicitement representes
3718	c
3719	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
3720	c
3721	c le thermique est supposé homogène et dissipé par mélange avec
3722	c son environnement. la longueur l_mix contrôle l'efficacité du
3723	c mélange
3724	c
3725	c Le calcul du transport des différentes espèces se fait en prenant
3726	c en compte:
3727	c 1. un flux de masse montant
3728	c 2. un flux de masse descendant
3729	c 3. un entrainement
3730	c 4. un detrainement
3731	c
3732	c=======================================================================
3733
3734	c-----------------------------------------------------------------------
3735	c declarations:
3736	c -------------
3737
3738	#include "dimensions.h"
3739	cccc#include "dimphy.h"
3740	#include "YOMCST.h"
3741
3742	c arguments:
3743	c ----------
3744
3745	INTEGER ngrid,nlay,w2di,tho
3746	real ptimestep,l_mix,r_aspect
3747	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
3748	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
3749	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
3750	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
3751	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
3752	real pphi(ngrid,nlay)
3753
3754	integer idetr
3755	save idetr
3756	data idetr/3/
3757
3758	c local:
3759	c ------
3760
3761	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
3762	real zsortie1d(klon)
3763	c CR: on remplace lmax(klon,klev+1)
3764	INTEGER lmax(klon),lmin(klon),lentr(klon)
3765	real linter(klon)
3766	real zmix(klon), fracazmix(klon)
3767	c RC
3768	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
3769
3770	real zlev(klon,klev+1),zlay(klon,klev)
3771	REAL zh(klon,klev),zdhadj(klon,klev)
3772	REAL ztv(klon,klev)
3773	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
3774	REAL wh(klon,klev+1)
3775	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
3776	real zla(klon,klev+1)
3777	real zwa(klon,klev+1)
3778	real zld(klon,klev+1)
3779	real zwd(klon,klev+1)
3780	real zsortie(klon,klev)
3781	real zva(klon,klev)
3782	real zua(klon,klev)
3783	real zoa(klon,klev)
3784
3785	real zha(klon,klev)
3786	real wa_moy(klon,klev+1)
3787	real fraca(klon,klev+1)
3788	real fracc(klon,klev+1)
3789	real zf,zf2
3790	real thetath2(klon,klev),wth2(klon,klev)
3791	! common/comtherm/thetath2,wth2
3792
3793	real count_time
3794	integer isplit,nsplit,ialt
3795	parameter (nsplit=10)
3796	data isplit/0/
3797	save isplit
3798
3799	logical sorties
3800	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
3801	real zpspsk(klon,klev)
3802
3803	c real wmax(klon,klev),wmaxa(klon)
3804	real wmax(klon),wmaxa(klon)
3805	real wa(klon,klev,klev+1)
3806	real wd(klon,klev+1)
3807	real larg_part(klon,klev,klev+1)
3808	real fracd(klon,klev+1)
3809	real xxx(klon,klev+1)
3810	real larg_cons(klon,klev+1)
3811	real larg_detr(klon,klev+1)
3812	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
3813	real pu_therm(klon,klev),pv_therm(klon,klev)
3814	real fm(klon,klev+1),entr(klon,klev)
3815	real fmc(klon,klev+1)
3816
3817	cCR:nouvelles variables
3818	real f_star(klon,klev+1),entr_star(klon,klev)
3819	real entr_star_tot(klon),entr_star2(klon)
3820	real f(klon), f0(klon)
3821	real zlevinter(klon)
3822	logical first
3823	data first /.false./
3824	save first
3825	cRC
3826
3827	character*2 str2
3828	character*10 str10
3829
3830	LOGICAL vtest(klon),down
3831
3832	EXTERNAL SCOPY
3833
3834	integer ncorrec,ll
3835	save ncorrec
3836	data ncorrec/0/
3837
3838	c
3839	c-----------------------------------------------------------------------
3840	c initialisation:
3841	c ---------------
3842	c
3843	sorties=.true.
3844	IF(ngrid.NE.klon) THEN
3845	PRINT*
3846	PRINT*,'STOP dans convadj'
3847	PRINT*,'ngrid =',ngrid
3848	PRINT*,'klon =',klon
3849	ENDIF
3850	c
3851	c-----------------------------------------------------------------------
3852	c incrementation eventuelle de tendances precedentes:
3853	c ---------------------------------------------------
3854
3855	print*,'0 OK convect8'
3856
3857	DO 1010 l=1,nlay
3858	DO 1015 ig=1,ngrid
3859	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
3860	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
3861	zu(ig,l)=pu(ig,l)
3862	zv(ig,l)=pv(ig,l)
3863	zo(ig,l)=po(ig,l)
3864	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
3865	1015 CONTINUE
3866	1010 CONTINUE
3867
3868	print*,'1 OK convect8'
3869	c --------------------
3870	c
3871	c
3872	c + + + + + + + + + + +
3873	c
3874	c
3875	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
3876	c wh,wt,wo ...
3877	c
3878	c + + + + + + + + + + + zh,zu,zv,zo,rho
3879	c
3880	c
3881	c -------------------- zlev(1)
3882	c \\\\\\\\\\\\\\\\\\\\
3883	c
3884	c
3885
3886	c-----------------------------------------------------------------------
3887	c Calcul des altitudes des couches
3888	c-----------------------------------------------------------------------
3889
3890	do l=2,nlay
3891	do ig=1,ngrid
3892	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
3893	enddo
3894	enddo
3895	do ig=1,ngrid
3896	zlev(ig,1)=0.
3897	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
3898	enddo
3899	do l=1,nlay
3900	do ig=1,ngrid
3901	zlay(ig,l)=pphi(ig,l)/RG
3902	enddo
3903	enddo
3904
3905	c print*,'2 OK convect8'
3906	c-----------------------------------------------------------------------
3907	c Calcul des densites
3908	c-----------------------------------------------------------------------
3909
3910	do l=1,nlay
3911	do ig=1,ngrid
3912	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
3913	enddo
3914	enddo
3915
3916	do l=2,nlay
3917	do ig=1,ngrid
3918	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
3919	enddo
3920	enddo
3921
3922	do k=1,nlay
3923	do l=1,nlay+1
3924	do ig=1,ngrid
3925	wa(ig,k,l)=0.
3926	enddo
3927	enddo
3928	enddo
3929
3930	c print*,'3 OK convect8'
3931	c------------------------------------------------------------------
3932	c Calcul de w2, quarre de w a partir de la cape
3933	c a partir de w2, on calcule wa, vitesse de l'ascendance
3934	c
3935	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
3936	c w2 est stoke dans wa
3937	c
3938	c ATTENTION: dans convect8, on n'utilise le calcule des wa
3939	c independants par couches que pour calculer l'entrainement
3940	c a la base et la hauteur max de l'ascendance.
3941	c
3942	c Indicages:
3943	c l'ascendance provenant du niveau k traverse l'interface l avec
3944	c une vitesse wa(k,l).
3945	c
3946	c --------------------
3947	c
3948	c + + + + + + + + + +
3949	c
3950	c wa(k,l) ---- -------------------- l
3951	c /\
3952	c /\|\|\ + + + + + + + + + +
3953	c \|\|
3954	c \|\| --------------------
3955	c \|\|
3956	c \|\| + + + + + + + + + +
3957	c \|\|
3958	c \|\| --------------------
3959	c \|\|__
3960	c \|___ + + + + + + + + + + k
3961	c
3962	c --------------------
3963	c
3964	c
3965	c
3966	c------------------------------------------------------------------
3967
3968	cCR: ponderation entrainement des couches instables
3969	cdef des entr_star tels que entr=f*entr_star
3970	do l=1,klev
3971	do ig=1,ngrid
3972	entr_star(ig,l)=0.
3973	enddo
3974	enddo
3975	c determination de la longueur de la couche d entrainement
3976	do ig=1,ngrid
3977	lentr(ig)=1
3978	enddo
3979
3980	con ne considere que les premieres couches instables
3981	do k=nlay-2,1,-1
3982	do ig=1,ngrid
3983	if (ztv(ig,k).gt.ztv(ig,k+1).and.
3984	s ztv(ig,k+1).le.ztv(ig,k+2)) then
3985	lentr(ig)=k
3986	endif
3987	enddo
3988	enddo
3989
3990	c determination du lmin: couche d ou provient le thermique
3991	do ig=1,ngrid
3992	lmin(ig)=1
3993	enddo
3994	do ig=1,ngrid
3995	do l=nlay,2,-1
3996	if (ztv(ig,l-1).gt.ztv(ig,l)) then
3997	lmin(ig)=l-1
3998	endif
3999	enddo
4000	enddo
4001	c
4002	c definition de l'entrainement des couches
4003	do l=1,klev-1
4004	do ig=1,ngrid
4005	if (ztv(ig,l).gt.ztv(ig,l+1).and.
4006	s l.ge.lmin(ig).and.l.le.lentr(ig)) then
4007	entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
4008	s (zlev(ig,l+1)-zlev(ig,l))
4009	endif
4010	enddo
4011	enddo
4012	c pas de thermique si couches 1->5 stables
4013	do ig=1,ngrid
4014	if (lmin(ig).gt.5) then
4015	do l=1,klev
4016	entr_star(ig,l)=0.
4017	enddo
4018	endif
4019	enddo
4020	c calcul de l entrainement total
4021	do ig=1,ngrid
4022	entr_star_tot(ig)=0.
4023	enddo
4024	do ig=1,ngrid
4025	do k=1,klev
4026	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
4027	enddo
4028	enddo
4029	c
4030	print*,'fin calcul entr_star'
4031	do k=1,klev
4032	do ig=1,ngrid
4033	ztva(ig,k)=ztv(ig,k)
4034	enddo
4035	enddo
4036	cRC
4037	c print*,'7 OK convect8'
4038	do k=1,klev+1
4039	do ig=1,ngrid
4040	zw2(ig,k)=0.
4041	fmc(ig,k)=0.
4042	cCR
4043	f_star(ig,k)=0.
4044	cRC
4045	larg_cons(ig,k)=0.
4046	larg_detr(ig,k)=0.
4047	wa_moy(ig,k)=0.
4048	enddo
4049	enddo
4050
4051	c print*,'8 OK convect8'
4052	do ig=1,ngrid
4053	linter(ig)=1.
4054	lmaxa(ig)=1
4055	lmix(ig)=1
4056	wmaxa(ig)=0.
4057	enddo
4058
4059	cCR:
4060	do l=1,nlay-2
4061	do ig=1,ngrid
4062	if (ztv(ig,l).gt.ztv(ig,l+1)
4063	s .and.entr_star(ig,l).gt.1.e-10
4064	s .and.zw2(ig,l).lt.1e-10) then
4065	f_star(ig,l+1)=entr_star(ig,l)
4066	ctest:calcul de dteta
4067	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
4068	s *(zlev(ig,l+1)-zlev(ig,l))
4069	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
4070	larg_detr(ig,l)=0.
4071	else if ((zw2(ig,l).ge.1e-10).and.
4072	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
4073	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
4074	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
4075	s *ztv(ig,l))/f_star(ig,l+1)
4076	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
4077	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
4078	s *(zlev(ig,l+1)-zlev(ig,l))
4079	endif
4080	c determination de zmax continu par interpolation lineaire
4081	if (zw2(ig,l+1).lt.0.) then
4082	ctest
4083	if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
4084	print*,'pb linter'
4085	endif
4086	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
4087	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
4088	zw2(ig,l+1)=0.
4089	lmaxa(ig)=l
4090	else
4091	if (zw2(ig,l+1).lt.0.) then
4092	print*,'pb1 zw2<0'
4093	endif
4094	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
4095	endif
4096	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
4097	c lmix est le niveau de la couche ou w (wa_moy) est maximum
4098	lmix(ig)=l+1
4099	wmaxa(ig)=wa_moy(ig,l+1)
4100	endif
4101	enddo
4102	enddo
4103	print*,'fin calcul zw2'
4104	c
4105	c Calcul de la couche correspondant a la hauteur du thermique
4106	do ig=1,ngrid
4107	lmax(ig)=lentr(ig)
4108	enddo
4109	do ig=1,ngrid
4110	do l=nlay,lentr(ig)+1,-1
4111	if (zw2(ig,l).le.1.e-10) then
4112	lmax(ig)=l-1
4113	endif
4114	enddo
4115	enddo
4116	c pas de thermique si couches 1->5 stables
4117	do ig=1,ngrid
4118	if (lmin(ig).gt.5) then
4119	lmax(ig)=1
4120	lmin(ig)=1
4121	endif
4122	enddo
4123	c
4124	c Determination de zw2 max
4125	do ig=1,ngrid
4126	wmax(ig)=0.
4127	enddo
4128
4129	do l=1,nlay
4130	do ig=1,ngrid
4131	if (l.le.lmax(ig)) then
4132	if (zw2(ig,l).lt.0.)then
4133	print*,'pb2 zw2<0'
4134	endif
4135	zw2(ig,l)=sqrt(zw2(ig,l))
4136	wmax(ig)=max(wmax(ig),zw2(ig,l))
4137	else
4138	zw2(ig,l)=0.
4139	endif
4140	enddo
4141	enddo
4142
4143	c Longueur caracteristique correspondant a la hauteur des thermiques.
4144	do ig=1,ngrid
4145	zmax(ig)=0.
4146	zlevinter(ig)=zlev(ig,1)
4147	enddo
4148	do ig=1,ngrid
4149	c calcul de zlevinter
4150	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
4151	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
4152	s -zlev(ig,lmax(ig)))
4153	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
4154	enddo
4155
4156	print*,'avant fermeture'
4157	c Fermeture,determination de f
4158	do ig=1,ngrid
4159	entr_star2(ig)=0.
4160	enddo
4161	do ig=1,ngrid
4162	if (entr_star_tot(ig).LT.1.e-10) then
4163	f(ig)=0.
4164	else
4165	do k=lmin(ig),lentr(ig)
4166	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
4167	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
4168	enddo
4169	c Nouvelle fermeture
4170	f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
4171	s entr_star2(ig))entr_star_tot(ig)
4172	ctest
4173	c if (first) then
4174	c f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
4175	c s *wmax(ig))
4176	c endif
4177	endif
4178	c f0(ig)=f(ig)
4179	c first=.true.
4180	enddo
4181	print*,'apres fermeture'
4182
4183	c Calcul de l'entrainement
4184	do k=1,klev
4185	do ig=1,ngrid
4186	entr(ig,k)=f(ig)*entr_star(ig,k)
4187	enddo
4188	enddo
4189	c Calcul des flux
4190	do ig=1,ngrid
4191	do l=1,lmax(ig)-1
4192	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
4193	enddo
4194	enddo
4195
4196	cRC
4197
4198
4199	c print*,'9 OK convect8'
4200	c print*,'WA1 ',wa_moy
4201
4202	c determination de l'indice du debut de la mixed layer ou w decroit
4203
4204	c calcul de la largeur de chaque ascendance dans le cas conservatif.
4205	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
4206	c d'une couche est égale à la hauteur de la couche alimentante.
4207	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
4208	c de la vitesse d'entrainement horizontal dans la couche alimentante.
4209
4210	do l=2,nlay
4211	do ig=1,ngrid
4212	if (l.le.lmaxa(ig)) then
4213	zw=max(wa_moy(ig,l),1.e-10)
4214	larg_cons(ig,l)=zmax(ig)*r_aspect
4215	s fmc(ig,l)/(rhobarz(ig,l)zw)
4216	endif
4217	enddo
4218	enddo
4219
4220	do l=2,nlay
4221	do ig=1,ngrid
4222	if (l.le.lmaxa(ig)) then
4223	c if (idetr.eq.0) then
4224	c cette option est finalement en dur.
4225	if ((l_mix*zlev(ig,l)).lt.0.)then
4226	print,'pb l_mixzlev<0'
4227	endif
4228	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4229	c else if (idetr.eq.1) then
4230	c larg_detr(ig,l)=larg_cons(ig,l)
4231	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
4232	c else if (idetr.eq.2) then
4233	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4234	c s *sqrt(wa_moy(ig,l))
4235	c else if (idetr.eq.4) then
4236	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4237	c s *wa_moy(ig,l)
4238	c endif
4239	endif
4240	enddo
4241	enddo
4242
4243	c print*,'10 OK convect8'
4244	c print*,'WA2 ',wa_moy
4245	c calcul de la fraction de la maille concernée par l'ascendance en tenant
4246	c compte de l'epluchage du thermique.
4247	c
4248	cCR def de zmix continu (profil parabolique des vitesses)
4249	do ig=1,ngrid
4250	if (lmix(ig).gt.1.) then
4251	c test
4252	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
4253	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
4254	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
4255	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
4256	s then
4257	c
4258	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
4259	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
4260	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
4261	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
4262	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
4263	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
4264	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
4265	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
4266	else
4267	zmix(ig)=zlev(ig,lmix(ig))
4268	print*,'pb zmix'
4269	endif
4270	else
4271	zmix(ig)=0.
4272	endif
4273	ctest
4274	if ((zmax(ig)-zmix(ig)).lt.0.) then
4275	zmix(ig)=0.99*zmax(ig)
4276	c print*,'pb zmix>zmax'
4277	endif
4278	enddo
4279	c
4280	c calcul du nouveau lmix correspondant
4281	do ig=1,ngrid
4282	do l=1,klev
4283	if (zmix(ig).ge.zlev(ig,l).and.
4284	s zmix(ig).lt.zlev(ig,l+1)) then
4285	lmix(ig)=l
4286	endif
4287	enddo
4288	enddo
4289	c
4290	do l=2,nlay
4291	do ig=1,ngrid
4292	if(larg_cons(ig,l).gt.1.) then
4293	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
4294	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
4295	s /(r_aspect*zmax(ig))
4296	c test
4297	fraca(ig,l)=max(fraca(ig,l),0.)
4298	fraca(ig,l)=min(fraca(ig,l),0.5)
4299	fracd(ig,l)=1.-fraca(ig,l)
4300	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
4301	else
4302	c wa_moy(ig,l)=0.
4303	fraca(ig,l)=0.
4304	fracc(ig,l)=0.
4305	fracd(ig,l)=1.
4306	endif
4307	enddo
4308	enddo
4309	cCR: calcul de fracazmix
4310	do ig=1,ngrid
4311	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
4312	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
4313	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
4314	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
4315	enddo
4316	c
4317	do l=2,nlay
4318	do ig=1,ngrid
4319	if(larg_cons(ig,l).gt.1.) then
4320	if (l.gt.lmix(ig)) then
4321	ctest
4322	if (zmax(ig)-zmix(ig).lt.1.e-10) then
4323	c print*,'pb xxx'
4324	xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
4325	else
4326	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
4327	endif
4328	if (idetr.eq.0) then
4329	fraca(ig,l)=fracazmix(ig)
4330	else if (idetr.eq.1) then
4331	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
4332	else if (idetr.eq.2) then
4333	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
4334	else
4335	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
4336	endif
4337	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
4338	fraca(ig,l)=max(fraca(ig,l),0.)
4339	fraca(ig,l)=min(fraca(ig,l),0.5)
4340	fracd(ig,l)=1.-fraca(ig,l)
4341	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
4342	endif
4343	endif
4344	enddo
4345	enddo
4346
4347	print*,'fin calcul fraca'
4348	c print*,'11 OK convect8'
4349	c print*,'Ea3 ',wa_moy
4350	c------------------------------------------------------------------
4351	c Calcul de fracd, wd
4352	c somme wa - wd = 0
4353	c------------------------------------------------------------------
4354
4355
4356	do ig=1,ngrid
4357	fm(ig,1)=0.
4358	fm(ig,nlay+1)=0.
4359	enddo
4360
4361	do l=2,nlay
4362	do ig=1,ngrid
4363	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
4364	cCR:test
4365	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
4366	s .and.l.gt.lmix(ig)) then
4367	fm(ig,l)=fm(ig,l-1)
4368	c write(1,*)'ajustement fm, l',l
4369	endif
4370	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
4371	cRC
4372	enddo
4373	do ig=1,ngrid
4374	if(fracd(ig,l).lt.0.1) then
4375	stop'fracd trop petit'
4376	else
4377	c vitesse descendante "diagnostique"
4378	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
4379	endif
4380	enddo
4381	enddo
4382
4383	do l=1,nlay
4384	do ig=1,ngrid
4385	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
4386	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
4387	enddo
4388	enddo
4389
4390	print*,'12 OK convect8'
4391	c print*,'WA4 ',wa_moy
4392	cc------------------------------------------------------------------
4393	c calcul du transport vertical
4394	c------------------------------------------------------------------
4395
4396	go to 4444
4397	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
4398	do l=2,nlay-1
4399	do ig=1,ngrid
4400	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
4401	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
4402	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
4403	c s ,fm(ig,l+1)*ptimestep
4404	c s ,' M=',masse(ig,l),masse(ig,l+1)
4405	endif
4406	enddo
4407	enddo
4408
4409	do l=1,nlay
4410	do ig=1,ngrid
4411	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
4412	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
4413	c s ,entr(ig,l)*ptimestep
4414	c s ,' M=',masse(ig,l)
4415	endif
4416	enddo
4417	enddo
4418
4419	do l=1,nlay
4420	do ig=1,ngrid
4421	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
4422	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
4423	c s ,' FM=',fm(ig,l)
4424	endif
4425	if(.not.masse(ig,l).ge.1.e-10
4426	s .or..not.masse(ig,l).le.1.e4) then
4427	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
4428	c s ,' M=',masse(ig,l)
4429	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
4430	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
4431	c print*,'zlev(ig,l+1),zlev(ig,l)'
4432	c s ,zlev(ig,l+1),zlev(ig,l)
4433	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
4434	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
4435	endif
4436	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
4437	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
4438	c s ,' E=',entr(ig,l)
4439	endif
4440	enddo
4441	enddo
4442
4443	4444 continue
4444
4445	cCR:redefinition du entr
4446	do l=1,nlay
4447	do ig=1,ngrid
4448	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4449	if (detr(ig,l).lt.0.) then
4450	entr(ig,l)=entr(ig,l)-detr(ig,l)
4451	detr(ig,l)=0.
4452	c print*,'WARNING !!! detrainement negatif ',ig,l
4453	endif
4454	enddo
4455	enddo
4456	cRC
4457	if (w2di.eq.1) then
4458	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
4459	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
4460	else
4461	fm0=fm
4462	entr0=entr
4463	endif
4464
4465	if (1.eq.1) then
4466	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4467	. ,zh,zdhadj,zha)
4468	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4469	. ,zo,pdoadj,zoa)
4470	else
4471	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
4472	. ,zh,zdhadj,zha)
4473	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
4474	. ,zo,pdoadj,zoa)
4475	endif
4476
4477	if (1.eq.0) then
4478	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
4479	. ,fraca,zmax
4480	. ,zu,zv,pduadj,pdvadj,zua,zva)
4481	else
4482	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4483	. ,zu,pduadj,zua)
4484	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
4485	. ,zv,pdvadj,zva)
4486	endif
4487
4488	do l=1,nlay
4489	do ig=1,ngrid
4490	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
4491	zf2=zf/(1.-zf)
4492	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
4493	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
4494	enddo
4495	enddo
4496
4497
4498
4499	c print*,'13 OK convect8'
4500	c print*,'WA5 ',wa_moy
4501	do l=1,nlay
4502	do ig=1,ngrid
4503	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
4504	enddo
4505	enddo
4506
4507
4508	c do l=1,nlay
4509	c do ig=1,ngrid
4510	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
4511	c print*,'WARN!!! ig=',ig,' l=',l
4512	c s ,' pdtadj=',pdtadj(ig,l)
4513	c endif
4514	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
4515	c print*,'WARN!!! ig=',ig,' l=',l
4516	c s ,' pdoadj=',pdoadj(ig,l)
4517	c endif
4518	c enddo
4519	c enddo
4520
4521	print*,'14 OK convect8'
4522	c------------------------------------------------------------------
4523	c Calculs pour les sorties
4524	c------------------------------------------------------------------
4525
4526	if(sorties) then
4527	do l=1,nlay
4528	do ig=1,ngrid
4529	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
4530	zld(ig,l)=fracd(ig,l)*zmax(ig)
4531	if(1.-fracd(ig,l).gt.1.e-10)
4532	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
4533	enddo
4534	enddo
4535
4536	cdeja fait
4537	c do l=1,nlay
4538	c do ig=1,ngrid
4539	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4540	c if (detr(ig,l).lt.0.) then
4541	c entr(ig,l)=entr(ig,l)-detr(ig,l)
4542	c detr(ig,l)=0.
4543	c print*,'WARNING !!! detrainement negatif ',ig,l
4544	c endif
4545	c enddo
4546	c enddo
4547
4548	c print*,'15 OK convect8'
4549
4550	isplit=isplit+1
4551
4552
4553	c #define und
4554	goto 123
4555	#ifdef und
4556	CALL writeg1d(1,nlay,wd,'wd ','wd ')
4557	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
4558	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
4559	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
4560	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
4561	CALL writeg1d(1,nlay,zla,'la ','la ')
4562	CALL writeg1d(1,nlay,zld,'ld ','ld ')
4563	CALL writeg1d(1,nlay,pt,'pt ','pt ')
4564	CALL writeg1d(1,nlay,zh,'zh ','zh ')
4565	CALL writeg1d(1,nlay,zha,'zha ','zha ')
4566	CALL writeg1d(1,nlay,zu,'zu ','zu ')
4567	CALL writeg1d(1,nlay,zv,'zv ','zv ')
4568	CALL writeg1d(1,nlay,zo,'zo ','zo ')
4569	CALL writeg1d(1,nlay,wh,'wh ','wh ')
4570	CALL writeg1d(1,nlay,wu,'wu ','wu ')
4571	CALL writeg1d(1,nlay,wv,'wv ','wv ')
4572	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
4573	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
4574	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
4575	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
4576	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
4577	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
4578	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
4579	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
4580
4581	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
4582	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
4583	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
4584
4585	c recalcul des flux en diagnostique...
4586	c print*,'PAS DE TEMPS ',ptimestep
4587	call dt2F(pplev,pplay,pt,pdtadj,wh)
4588	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
4589	#endif
4590	123 continue
4591	#define troisD
4592	#ifdef troisD
4593	c if (sorties) then
4594	print*,'Debut des wrgradsfi'
4595
4596	c print*,'16 OK convect8'
4597	call wrgradsfi(1,nlay,wd,'wd ','wd ')
4598	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
4599	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
4600	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
4601	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
4602	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
4603	c print*,'WA6 ',wa_moy
4604	call wrgradsfi(1,nlay,zla,'la ','la ')
4605	call wrgradsfi(1,nlay,zld,'ld ','ld ')
4606	call wrgradsfi(1,nlay,pt,'pt ','pt ')
4607	call wrgradsfi(1,nlay,zh,'zh ','zh ')
4608	call wrgradsfi(1,nlay,zha,'zha ','zha ')
4609	call wrgradsfi(1,nlay,zua,'zua ','zua ')
4610	call wrgradsfi(1,nlay,zva,'zva ','zva ')
4611	call wrgradsfi(1,nlay,zu,'zu ','zu ')
4612	call wrgradsfi(1,nlay,zv,'zv ','zv ')
4613	call wrgradsfi(1,nlay,zo,'zo ','zo ')
4614	call wrgradsfi(1,nlay,wh,'wh ','wh ')
4615	call wrgradsfi(1,nlay,wu,'wu ','wu ')
4616	call wrgradsfi(1,nlay,wv,'wv ','wv ')
4617	call wrgradsfi(1,nlay,wo,'wo ','wo ')
4618	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
4619	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
4620	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
4621	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
4622	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
4623	call wrgradsfi(1,nlay,entr,'entr ','entr ')
4624	call wrgradsfi(1,nlay,detr,'detr ','detr ')
4625	call wrgradsfi(1,nlay,fm,'fm ','fm ')
4626	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
4627	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
4628	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
4629	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
4630
4631	call wrgradsfi(1,nlay,zo,'zo ','zo ')
4632	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
4633	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
4634
4635	cCR:nouveaux diagnostiques
4636	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
4637	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
4638	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
4639	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
4640	zsortie1d(:)=lmax(:)
4641	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
4642	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
4643	zsortie1d(:)=lmix(:)
4644	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
4645	zsortie1d(:)=lentr(:)
4646	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
4647
4648	c print*,'17 OK convect8'
4649
4650	do k=1,klev/10
4651	write(str2,'(i2.2)') k
4652	str10='wa'//str2
4653	do l=1,nlay
4654	do ig=1,ngrid
4655	zsortie(ig,l)=wa(ig,k,l)
4656	enddo
4657	enddo
4658	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
4659	do l=1,nlay
4660	do ig=1,ngrid
4661	zsortie(ig,l)=larg_part(ig,k,l)
4662	enddo
4663	enddo
4664	str10='la'//str2
4665	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
4666	enddo
4667
4668
4669	c print*,'18 OK convect8'
4670	c endif
4671	print*,'Fin des wrgradsfi'
4672	#endif
4673
4674	endif
4675
4676	c if(wa_moy(1,4).gt.1.e-10) stop
4677
4678	print*,'19 OK convect8'
4679	return
4680	end
4681
4682	subroutine dqthermcell(ngrid,nlay,ptimestep,fm,entr,
4683	. masse,q,dq,qa)
4684	USE dimphy
4685	implicit none
4686
4687	c=======================================================================
4688	c
4689	c Calcul du transport verticale dans la couche limite en presence
4690	c de "thermiques" explicitement representes
4691	c calcul du dq/dt une fois qu'on connait les ascendances
4692	c
4693	c=======================================================================
4694
4695	#include "dimensions.h"
4696	cccc#include "dimphy.h"
4697
4698	integer ngrid,nlay
4699
4700	real ptimestep
4701	real masse(ngrid,nlay),fm(ngrid,nlay+1)
4702	real entr(ngrid,nlay)
4703	real q(ngrid,nlay)
4704	real dq(ngrid,nlay)
4705
4706	real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
4707
4708	integer ig,k
4709
4710	c calcul du detrainement
4711
4712	do k=1,nlay
4713	do ig=1,ngrid
4714	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
4715	ctest
4716	if (detr(ig,k).lt.0.) then
4717	entr(ig,k)=entr(ig,k)-detr(ig,k)
4718	detr(ig,k)=0.
4719	c print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
4720	c s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
4721	endif
4722	if (fm(ig,k+1).lt.0.) then
4723	c print*,'fm2<0!!!'
4724	endif
4725	if (entr(ig,k).lt.0.) then
4726	c print*,'entr2<0!!!'
4727	endif
4728	enddo
4729	enddo
4730
4731	c calcul de la valeur dans les ascendances
4732	do ig=1,ngrid
4733	qa(ig,1)=q(ig,1)
4734	enddo
4735
4736	do k=2,nlay
4737	do ig=1,ngrid
4738	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
4739	s 1.e-5*masse(ig,k)) then
4740	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k))
4741	s /(fm(ig,k+1)+detr(ig,k))
4742	else
4743	qa(ig,k)=q(ig,k)
4744	endif
4745	if (qa(ig,k).lt.0.) then
4746	c print*,'qa<0!!!'
4747	endif
4748	if (q(ig,k).lt.0.) then
4749	c print*,'q<0!!!'
4750	endif
4751	enddo
4752	enddo
4753
4754	do k=2,nlay
4755	do ig=1,ngrid
4756	c wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
4757	wqd(ig,k)=fm(ig,k)*q(ig,k)
4758	if (wqd(ig,k).lt.0.) then
4759	c print*,'wqd<0!!!'
4760	endif
4761	enddo
4762	enddo
4763	do ig=1,ngrid
4764	wqd(ig,1)=0.
4765	wqd(ig,nlay+1)=0.
4766	enddo
4767
4768	do k=1,nlay
4769	do ig=1,ngrid
4770	dq(ig,k)=(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k)
4771	s -wqd(ig,k)+wqd(ig,k+1))
4772	s /masse(ig,k)
4773	c if (dq(ig,k).lt.0.) then
4774	c print*,'dq<0!!!'
4775	c endif
4776	enddo
4777	enddo
4778
4779	return
4780	end
4781	subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
4782	. ,fraca,larga
4783	. ,u,v,du,dv,ua,va)
4784	USE dimphy
4785	implicit none
4786
4787	c=======================================================================
4788	c
4789	c Calcul du transport verticale dans la couche limite en presence
4790	c de "thermiques" explicitement representes
4791	c calcul du dq/dt une fois qu'on connait les ascendances
4792	c
4793	c=======================================================================
4794
4795	#include "dimensions.h"
4796	cccc#include "dimphy.h"
4797
4798	integer ngrid,nlay
4799
4800	real ptimestep
4801	real masse(ngrid,nlay),fm(ngrid,nlay+1)
4802	real fraca(ngrid,nlay+1)
4803	real larga(ngrid)
4804	real entr(ngrid,nlay)
4805	real u(ngrid,nlay)
4806	real ua(ngrid,nlay)
4807	real du(ngrid,nlay)
4808	real v(ngrid,nlay)
4809	real va(ngrid,nlay)
4810	real dv(ngrid,nlay)
4811
4812	real qa(klon,klev),detr(klon,klev)
4813	real wvd(klon,klev+1),wud(klon,klev+1)
4814	real gamma0,gamma(klon,klev+1)
4815	real dua,dva
4816	integer iter
4817
4818	integer ig,k
4819
4820	c calcul du detrainement
4821
4822	do k=1,nlay
4823	do ig=1,ngrid
4824	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
4825	enddo
4826	enddo
4827
4828	c calcul de la valeur dans les ascendances
4829	do ig=1,ngrid
4830	ua(ig,1)=u(ig,1)
4831	va(ig,1)=v(ig,1)
4832	enddo
4833
4834	do k=2,nlay
4835	do ig=1,ngrid
4836	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
4837	s 1.e-5*masse(ig,k)) then
4838	c On itère sur la valeur du coeff de freinage.
4839	c gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4840	gamma0=masse(ig,k)
4841	s sqrt( 0.5(fraca(ig,k+1)+fraca(ig,k)) )
4842	s *0.5/larga(ig)
4843	c gamma0=0.
4844	c la première fois on multiplie le coefficient de freinage
4845	c par le module du vent dans la couche en dessous.
4846	dua=ua(ig,k-1)-u(ig,k-1)
4847	dva=va(ig,k-1)-v(ig,k-1)
4848	do iter=1,5
4849	gamma(ig,k)=gamma0sqrt(dua2+dva*2)
4850	ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
4851	s +(entr(ig,k)+gamma(ig,k))*u(ig,k))
4852	s /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4853	va(ig,k)=(fm(ig,k)*va(ig,k-1)
4854	s +(entr(ig,k)+gamma(ig,k))*v(ig,k))
4855	s /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4856	c print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4857	dua=ua(ig,k)-u(ig,k)
4858	dva=va(ig,k)-v(ig,k)
4859	enddo
4860	else
4861	ua(ig,k)=u(ig,k)
4862	va(ig,k)=v(ig,k)
4863	gamma(ig,k)=0.
4864	endif
4865	enddo
4866	enddo
4867
4868	do k=2,nlay
4869	do ig=1,ngrid
4870	wud(ig,k)=fm(ig,k)*u(ig,k)
4871	wvd(ig,k)=fm(ig,k)*v(ig,k)
4872	enddo
4873	enddo
4874	do ig=1,ngrid
4875	wud(ig,1)=0.
4876	wud(ig,nlay+1)=0.
4877	wvd(ig,1)=0.
4878	wvd(ig,nlay+1)=0.
4879	enddo
4880
4881	do k=1,nlay
4882	do ig=1,ngrid
4883	du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
4884	s -(entr(ig,k)+gamma(ig,k))*u(ig,k)
4885	s -wud(ig,k)+wud(ig,k+1))
4886	s /masse(ig,k)
4887	dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
4888	s -(entr(ig,k)+gamma(ig,k))*v(ig,k)
4889	s -wvd(ig,k)+wvd(ig,k+1))
4890	s /masse(ig,k)
4891	enddo
4892	enddo
4893
4894	return
4895	end
4896	subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
4897	. ,q,dq,qa)
4898	USE dimphy
4899	implicit none
4900
4901	c=======================================================================
4902	c
4903	c Calcul du transport verticale dans la couche limite en presence
4904	c de "thermiques" explicitement representes
4905	c calcul du dq/dt une fois qu'on connait les ascendances
4906	c
4907	c=======================================================================
4908
4909	#include "dimensions.h"
4910	cccc#include "dimphy.h"
4911
4912	integer ngrid,nlay
4913
4914	real ptimestep
4915	real masse(ngrid,nlay),fm(ngrid,nlay+1)
4916	real entr(ngrid,nlay),frac(ngrid,nlay)
4917	real q(ngrid,nlay)
4918	real dq(ngrid,nlay)
4919
4920	real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
4921	real qe(klon,klev),zf,zf2
4922
4923	integer ig,k
4924
4925	c calcul du detrainement
4926
4927	do k=1,nlay
4928	do ig=1,ngrid
4929	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
4930	enddo
4931	enddo
4932
4933	c calcul de la valeur dans les ascendances
4934	do ig=1,ngrid
4935	qa(ig,1)=q(ig,1)
4936	qe(ig,1)=q(ig,1)
4937	enddo
4938
4939	do k=2,nlay
4940	do ig=1,ngrid
4941	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
4942	s 1.e-5*masse(ig,k)) then
4943	zf=0.5*(frac(ig,k)+frac(ig,k+1))
4944	zf2=1./(1.-zf)
4945	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+zf2entr(ig,k)*q(ig,k))
4946	s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2)
4947	qe(ig,k)=(q(ig,k)-zfqa(ig,k))zf2
4948	else
4949	qa(ig,k)=q(ig,k)
4950	qe(ig,k)=q(ig,k)
4951	endif
4952	enddo
4953	enddo
4954
4955	do k=2,nlay
4956	do ig=1,ngrid
4957	c wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
4958	wqd(ig,k)=fm(ig,k)*qe(ig,k)
4959	enddo
4960	enddo
4961	do ig=1,ngrid
4962	wqd(ig,1)=0.
4963	wqd(ig,nlay+1)=0.
4964	enddo
4965
4966	do k=1,nlay
4967	do ig=1,ngrid
4968	dq(ig,k)=(detr(ig,k)qa(ig,k)-entr(ig,k)qe(ig,k)
4969	s -wqd(ig,k)+wqd(ig,k+1))
4970	s /masse(ig,k)
4971	enddo
4972	enddo
4973
4974	return
4975	end
4976	subroutine dvthermcell2(ngrid,nlay,ptimestep,fm,entr,masse
4977	. ,fraca,larga
4978	. ,u,v,du,dv,ua,va)
4979	use dimphy
4980	implicit none
4981
4982	c=======================================================================
4983	c
4984	c Calcul du transport verticale dans la couche limite en presence
4985	c de "thermiques" explicitement representes
4986	c calcul du dq/dt une fois qu'on connait les ascendances
4987	c
4988	c=======================================================================
4989
4990	#include "dimensions.h"
4991	cccc#include "dimphy.h"
4992
4993	integer ngrid,nlay
4994
4995	real ptimestep
4996	real masse(ngrid,nlay),fm(ngrid,nlay+1)
4997	real fraca(ngrid,nlay+1)
4998	real larga(ngrid)
4999	real entr(ngrid,nlay)
5000	real u(ngrid,nlay)
5001	real ua(ngrid,nlay)
5002	real du(ngrid,nlay)
5003	real v(ngrid,nlay)
5004	real va(ngrid,nlay)
5005	real dv(ngrid,nlay)
5006
5007	real qa(klon,klev),detr(klon,klev),zf,zf2
5008	real wvd(klon,klev+1),wud(klon,klev+1)
5009	real gamma0,gamma(klon,klev+1)
5010	real ue(klon,klev),ve(klon,klev)
5011	real dua,dva
5012	integer iter
5013
5014	integer ig,k
5015
5016	c calcul du detrainement
5017
5018	do k=1,nlay
5019	do ig=1,ngrid
5020	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
5021	enddo
5022	enddo
5023
5024	c calcul de la valeur dans les ascendances
5025	do ig=1,ngrid
5026	ua(ig,1)=u(ig,1)
5027	va(ig,1)=v(ig,1)
5028	ue(ig,1)=u(ig,1)
5029	ve(ig,1)=v(ig,1)
5030	enddo
5031
5032	do k=2,nlay
5033	do ig=1,ngrid
5034	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
5035	s 1.e-5*masse(ig,k)) then
5036	c On itère sur la valeur du coeff de freinage.
5037	c gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
5038	gamma0=masse(ig,k)
5039	s sqrt( 0.5(fraca(ig,k+1)+fraca(ig,k)) )
5040	s *0.5/larga(ig)
5041	s *1.
5042	c s *0.5
5043	c gamma0=0.
5044	zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
5045	zf=0.
5046	zf2=1./(1.-zf)
5047	c la première fois on multiplie le coefficient de freinage
5048	c par le module du vent dans la couche en dessous.
5049	dua=ua(ig,k-1)-u(ig,k-1)
5050	dva=va(ig,k-1)-v(ig,k-1)
5051	do iter=1,5
5052	c On choisit une relaxation lineaire.
5053	gamma(ig,k)=gamma0
5054	c On choisit une relaxation quadratique.
5055	gamma(ig,k)=gamma0sqrt(dua2+dva*2)
5056	ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
5057	s +(zf2entr(ig,k)+gamma(ig,k))u(ig,k))
5058	s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2
5059	s +gamma(ig,k))
5060	va(ig,k)=(fm(ig,k)*va(ig,k-1)
5061	s +(zf2entr(ig,k)+gamma(ig,k))v(ig,k))
5062	s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2
5063	s +gamma(ig,k))
5064	c print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
5065	dua=ua(ig,k)-u(ig,k)
5066	dva=va(ig,k)-v(ig,k)
5067	ue(ig,k)=(u(ig,k)-zfua(ig,k))zf2
5068	ve(ig,k)=(v(ig,k)-zfva(ig,k))zf2
5069	enddo
5070	else
5071	ua(ig,k)=u(ig,k)
5072	va(ig,k)=v(ig,k)
5073	ue(ig,k)=u(ig,k)
5074	ve(ig,k)=v(ig,k)
5075	gamma(ig,k)=0.
5076	endif
5077	enddo
5078	enddo
5079
5080	do k=2,nlay
5081	do ig=1,ngrid
5082	wud(ig,k)=fm(ig,k)*ue(ig,k)
5083	wvd(ig,k)=fm(ig,k)*ve(ig,k)
5084	enddo
5085	enddo
5086	do ig=1,ngrid
5087	wud(ig,1)=0.
5088	wud(ig,nlay+1)=0.
5089	wvd(ig,1)=0.
5090	wvd(ig,nlay+1)=0.
5091	enddo
5092
5093	do k=1,nlay
5094	do ig=1,ngrid
5095	du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
5096	s -(entr(ig,k)+gamma(ig,k))*ue(ig,k)
5097	s -wud(ig,k)+wud(ig,k+1))
5098	s /masse(ig,k)
5099	dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
5100	s -(entr(ig,k)+gamma(ig,k))*ve(ig,k)
5101	s -wvd(ig,k)+wvd(ig,k+1))
5102	s /masse(ig,k)
5103	enddo
5104	enddo
5105
5106	return
5107	end
5108	SUBROUTINE thermcell_sec(ngrid,nlay,ptimestep
5109	s ,pplay,pplev,pphi,zlev
5110	s ,pu,pv,pt,po
5111	s ,pduadj,pdvadj,pdtadj,pdoadj
5112	s ,fm0,entr0
5113	c s ,pu_therm,pv_therm
5114	s ,r_aspect,l_mix,w2di,tho)
5115
5116	use dimphy
5117	IMPLICIT NONE
5118
5119	c=======================================================================
5120	c
5121	c Calcul du transport verticale dans la couche limite en presence
5122	c de "thermiques" explicitement representes
5123	c
5124	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
5125	c
5126	c le thermique est supposé homogène et dissipé par mélange avec
5127	c son environnement. la longueur l_mix contrôle l'efficacité du
5128	c mélange
5129	c
5130	c Le calcul du transport des différentes espèces se fait en prenant
5131	c en compte:
5132	c 1. un flux de masse montant
5133	c 2. un flux de masse descendant
5134	c 3. un entrainement
5135	c 4. un detrainement
5136	c
5137	c=======================================================================
5138
5139	c-----------------------------------------------------------------------
5140	c declarations:
5141	c -------------
5142
5143	#include "dimensions.h"
5144	cccc#include "dimphy.h"
5145	#include "YOMCST.h"
5146
5147	c arguments:
5148	c ----------
5149
5150	INTEGER ngrid,nlay,w2di,tho
5151	real ptimestep,l_mix,r_aspect
5152	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
5153	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
5154	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
5155	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
5156	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
5157	real pphi(ngrid,nlay)
5158
5159	integer idetr
5160	save idetr
5161	data idetr/3/
5162
5163	c local:
5164	c ------
5165
5166	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
5167	real zsortie1d(klon)
5168	c CR: on remplace lmax(klon,klev+1)
5169	INTEGER lmax(klon),lmin(klon),lentr(klon)
5170	real linter(klon)
5171	real zmix(klon), fracazmix(klon)
5172	c RC
5173	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
5174
5175	real zlev(klon,klev+1),zlay(klon,klev)
5176	REAL zh(klon,klev),zdhadj(klon,klev)
5177	REAL ztv(klon,klev)
5178	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
5179	REAL wh(klon,klev+1)
5180	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
5181	real zla(klon,klev+1)
5182	real zwa(klon,klev+1)
5183	real zld(klon,klev+1)
5184	real zwd(klon,klev+1)
5185	real zsortie(klon,klev)
5186	real zva(klon,klev)
5187	real zua(klon,klev)
5188	real zoa(klon,klev)
5189
5190	real zha(klon,klev)
5191	real wa_moy(klon,klev+1)
5192	real fraca(klon,klev+1)
5193	real fracc(klon,klev+1)
5194	real zf,zf2
5195	real thetath2(klon,klev),wth2(klon,klev)
5196	! common/comtherm/thetath2,wth2
5197
5198	real count_time
5199	integer isplit,nsplit,ialt
5200	parameter (nsplit=10)
5201	data isplit/0/
5202	save isplit
5203
5204	logical sorties
5205	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
5206	real zpspsk(klon,klev)
5207
5208	c real wmax(klon,klev),wmaxa(klon)
5209	real wmax(klon),wmaxa(klon)
5210	real wa(klon,klev,klev+1)
5211	real wd(klon,klev+1)
5212	real larg_part(klon,klev,klev+1)
5213	real fracd(klon,klev+1)
5214	real xxx(klon,klev+1)
5215	real larg_cons(klon,klev+1)
5216	real larg_detr(klon,klev+1)
5217	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
5218	real pu_therm(klon,klev),pv_therm(klon,klev)
5219	real fm(klon,klev+1),entr(klon,klev)
5220	real fmc(klon,klev+1)
5221
5222	cCR:nouvelles variables
5223	real f_star(klon,klev+1),entr_star(klon,klev)
5224	real entr_star_tot(klon),entr_star2(klon)
5225	real f(klon), f0(klon)
5226	real zlevinter(klon)
5227	logical first
5228	data first /.false./
5229	save first
5230	cRC
5231
5232	character*2 str2
5233	character*10 str10
5234
5235	LOGICAL vtest(klon),down
5236
5237	EXTERNAL SCOPY
5238
5239	integer ncorrec,ll
5240	save ncorrec
5241	data ncorrec/0/
5242
5243	c
5244	c-----------------------------------------------------------------------
5245	c initialisation:
5246	c ---------------
5247	c
5248	sorties=.true.
5249	IF(ngrid.NE.klon) THEN
5250	PRINT*
5251	PRINT*,'STOP dans convadj'
5252	PRINT*,'ngrid =',ngrid
5253	PRINT*,'klon =',klon
5254	ENDIF
5255	c
5256	c-----------------------------------------------------------------------
5257	c incrementation eventuelle de tendances precedentes:
5258	c ---------------------------------------------------
5259
5260	c print*,'0 OK convect8'
5261
5262	DO 1010 l=1,nlay
5263	DO 1015 ig=1,ngrid
5264	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
5265	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
5266	zu(ig,l)=pu(ig,l)
5267	zv(ig,l)=pv(ig,l)
5268	zo(ig,l)=po(ig,l)
5269	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
5270	1015 CONTINUE
5271	1010 CONTINUE
5272
5273	c print*,'1 OK convect8'
5274	c --------------------
5275	c
5276	c
5277	c + + + + + + + + + + +
5278	c
5279	c
5280	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
5281	c wh,wt,wo ...
5282	c
5283	c + + + + + + + + + + + zh,zu,zv,zo,rho
5284	c
5285	c
5286	c -------------------- zlev(1)
5287	c \\\\\\\\\\\\\\\\\\\\
5288	c
5289	c
5290
5291	c-----------------------------------------------------------------------
5292	c Calcul des altitudes des couches
5293	c-----------------------------------------------------------------------
5294
5295	do l=2,nlay
5296	do ig=1,ngrid
5297	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
5298	enddo
5299	enddo
5300	do ig=1,ngrid
5301	zlev(ig,1)=0.
5302	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
5303	enddo
5304	do l=1,nlay
5305	do ig=1,ngrid
5306	zlay(ig,l)=pphi(ig,l)/RG
5307	enddo
5308	enddo
5309
5310	c print*,'2 OK convect8'
5311	c-----------------------------------------------------------------------
5312	c Calcul des densites
5313	c-----------------------------------------------------------------------
5314
5315	do l=1,nlay
5316	do ig=1,ngrid
5317	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
5318	enddo
5319	enddo
5320
5321	do l=2,nlay
5322	do ig=1,ngrid
5323	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
5324	enddo
5325	enddo
5326
5327	do k=1,nlay
5328	do l=1,nlay+1
5329	do ig=1,ngrid
5330	wa(ig,k,l)=0.
5331	enddo
5332	enddo
5333	enddo
5334
5335	c print*,'3 OK convect8'
5336	c------------------------------------------------------------------
5337	c Calcul de w2, quarre de w a partir de la cape
5338	c a partir de w2, on calcule wa, vitesse de l'ascendance
5339	c
5340	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
5341	c w2 est stoke dans wa
5342	c
5343	c ATTENTION: dans convect8, on n'utilise le calcule des wa
5344	c independants par couches que pour calculer l'entrainement
5345	c a la base et la hauteur max de l'ascendance.
5346	c
5347	c Indicages:
5348	c l'ascendance provenant du niveau k traverse l'interface l avec
5349	c une vitesse wa(k,l).
5350	c
5351	c --------------------
5352	c
5353	c + + + + + + + + + +
5354	c
5355	c wa(k,l) ---- -------------------- l
5356	c /\
5357	c /\|\|\ + + + + + + + + + +
5358	c \|\|
5359	c \|\| --------------------
5360	c \|\|
5361	c \|\| + + + + + + + + + +
5362	c \|\|
5363	c \|\| --------------------
5364	c \|\|__
5365	c \|___ + + + + + + + + + + k
5366	c
5367	c --------------------
5368	c
5369	c
5370	c
5371	c------------------------------------------------------------------
5372
5373	cCR: ponderation entrainement des couches instables
5374	cdef des entr_star tels que entr=f*entr_star
5375	do l=1,klev
5376	do ig=1,ngrid
5377	entr_star(ig,l)=0.
5378	enddo
5379	enddo
5380	c determination de la longueur de la couche d entrainement
5381	do ig=1,ngrid
5382	lentr(ig)=1
5383	enddo
5384
5385	con ne considere que les premieres couches instables
5386	do k=nlay-2,1,-1
5387	do ig=1,ngrid
5388	if (ztv(ig,k).gt.ztv(ig,k+1).and.
5389	s ztv(ig,k+1).le.ztv(ig,k+2)) then
5390	lentr(ig)=k
5391	endif
5392	enddo
5393	enddo
5394
5395	c determination du lmin: couche d ou provient le thermique
5396	do ig=1,ngrid
5397	lmin(ig)=1
5398	enddo
5399	do ig=1,ngrid
5400	do l=nlay,2,-1
5401	if (ztv(ig,l-1).gt.ztv(ig,l)) then
5402	lmin(ig)=l-1
5403	endif
5404	enddo
5405	enddo
5406	c
5407	c definition de l'entrainement des couches
5408	do l=1,klev-1
5409	do ig=1,ngrid
5410	if (ztv(ig,l).gt.ztv(ig,l+1).and.
5411	s l.ge.lmin(ig).and.l.le.lentr(ig)) then
5412	entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
5413	c s (zlev(ig,l+1)-zlev(ig,l))
5414	s *sqrt(zlev(ig,l+1))
5415	endif
5416	enddo
5417	enddo
5418	c pas de thermique si couche 1 stable
5419	do ig=1,ngrid
5420	if (lmin(ig).gt.1) then
5421	do l=1,klev
5422	entr_star(ig,l)=0.
5423	enddo
5424	endif
5425	enddo
5426	c calcul de l entrainement total
5427	do ig=1,ngrid
5428	entr_star_tot(ig)=0.
5429	enddo
5430	do ig=1,ngrid
5431	do k=1,klev
5432	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
5433	enddo
5434	enddo
5435	c
5436	c print*,'fin calcul entr_star'
5437	do k=1,klev
5438	do ig=1,ngrid
5439	ztva(ig,k)=ztv(ig,k)
5440	enddo
5441	enddo
5442	cRC
5443	c print*,'7 OK convect8'
5444	do k=1,klev+1
5445	do ig=1,ngrid
5446	zw2(ig,k)=0.
5447	fmc(ig,k)=0.
5448	cCR
5449	f_star(ig,k)=0.
5450	cRC
5451	larg_cons(ig,k)=0.
5452	larg_detr(ig,k)=0.
5453	wa_moy(ig,k)=0.
5454	enddo
5455	enddo
5456
5457	c print*,'8 OK convect8'
5458	do ig=1,ngrid
5459	linter(ig)=1.
5460	lmaxa(ig)=1
5461	lmix(ig)=1
5462	wmaxa(ig)=0.
5463	enddo
5464
5465	cCR:
5466	do l=1,nlay-2
5467	do ig=1,ngrid
5468	if (ztv(ig,l).gt.ztv(ig,l+1)
5469	s .and.entr_star(ig,l).gt.1.e-10
5470	s .and.zw2(ig,l).lt.1e-10) then
5471	f_star(ig,l+1)=entr_star(ig,l)
5472	ctest:calcul de dteta
5473	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
5474	s *(zlev(ig,l+1)-zlev(ig,l))
5475	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
5476	larg_detr(ig,l)=0.
5477	else if ((zw2(ig,l).ge.1e-10).and.
5478	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
5479	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
5480	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
5481	s *ztv(ig,l))/f_star(ig,l+1)
5482	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
5483	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
5484	s *(zlev(ig,l+1)-zlev(ig,l))
5485	endif
5486	c determination de zmax continu par interpolation lineaire
5487	if (zw2(ig,l+1).lt.0.) then
5488	ctest
5489	if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
5490	c print*,'pb linter'
5491	endif
5492	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
5493	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
5494	zw2(ig,l+1)=0.
5495	lmaxa(ig)=l
5496	else
5497	if (zw2(ig,l+1).lt.0.) then
5498	c print*,'pb1 zw2<0'
5499	endif
5500	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
5501	endif
5502	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
5503	c lmix est le niveau de la couche ou w (wa_moy) est maximum
5504	lmix(ig)=l+1
5505	wmaxa(ig)=wa_moy(ig,l+1)
5506	endif
5507	enddo
5508	enddo
5509	c print*,'fin calcul zw2'
5510	c
5511	c Calcul de la couche correspondant a la hauteur du thermique
5512	do ig=1,ngrid
5513	lmax(ig)=lentr(ig)
5514	enddo
5515	do ig=1,ngrid
5516	do l=nlay,lentr(ig)+1,-1
5517	if (zw2(ig,l).le.1.e-10) then
5518	lmax(ig)=l-1
5519	endif
5520	enddo
5521	enddo
5522	c pas de thermique si couche 1 stable
5523	do ig=1,ngrid
5524	if (lmin(ig).gt.1) then
5525	lmax(ig)=1
5526	lmin(ig)=1
5527	endif
5528	enddo
5529	c
5530	c Determination de zw2 max
5531	do ig=1,ngrid
5532	wmax(ig)=0.
5533	enddo
5534
5535	do l=1,nlay
5536	do ig=1,ngrid
5537	if (l.le.lmax(ig)) then
5538	if (zw2(ig,l).lt.0.)then
5539	c print*,'pb2 zw2<0'
5540	endif
5541	zw2(ig,l)=sqrt(zw2(ig,l))
5542	wmax(ig)=max(wmax(ig),zw2(ig,l))
5543	else
5544	zw2(ig,l)=0.
5545	endif
5546	enddo
5547	enddo
5548
5549	c Longueur caracteristique correspondant a la hauteur des thermiques.
5550	do ig=1,ngrid
5551	zmax(ig)=0.
5552	zlevinter(ig)=zlev(ig,1)
5553	enddo
5554	do ig=1,ngrid
5555	c calcul de zlevinter
5556	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
5557	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
5558	s -zlev(ig,lmax(ig)))
5559	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
5560	enddo
5561
5562	c print*,'avant fermeture'
5563	c Fermeture,determination de f
5564	do ig=1,ngrid
5565	entr_star2(ig)=0.
5566	enddo
5567	do ig=1,ngrid
5568	if (entr_star_tot(ig).LT.1.e-10) then
5569	f(ig)=0.
5570	else
5571	do k=lmin(ig),lentr(ig)
5572	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
5573	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
5574	enddo
5575	c Nouvelle fermeture
5576	f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
5577	s entr_star2(ig))entr_star_tot(ig)
5578	ctest
5579	c if (first) then
5580	c f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
5581	c s *wmax(ig))
5582	c endif
5583	endif
5584	c f0(ig)=f(ig)
5585	c first=.true.
5586	enddo
5587	c print*,'apres fermeture'
5588
5589	c Calcul de l'entrainement
5590	do k=1,klev
5591	do ig=1,ngrid
5592	entr(ig,k)=f(ig)*entr_star(ig,k)
5593	enddo
5594	enddo
5595	cCR:test pour entrainer moins que la masse
5596	do ig=1,ngrid
5597	do l=1,lentr(ig)
5598	if ((entr(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
5599	entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
5600	s -0.9*masse(ig,l)/ptimestep
5601	entr(ig,l)=0.9*masse(ig,l)/ptimestep
5602	endif
5603	enddo
5604	enddo
5605	cCR: fin test
5606	c Calcul des flux
5607	do ig=1,ngrid
5608	do l=1,lmax(ig)-1
5609	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
5610	enddo
5611	enddo
5612
5613	cRC
5614
5615
5616	c print*,'9 OK convect8'
5617	c print*,'WA1 ',wa_moy
5618
5619	c determination de l'indice du debut de la mixed layer ou w decroit
5620
5621	c calcul de la largeur de chaque ascendance dans le cas conservatif.
5622	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
5623	c d'une couche est égale à la hauteur de la couche alimentante.
5624	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
5625	c de la vitesse d'entrainement horizontal dans la couche alimentante.
5626
5627	do l=2,nlay
5628	do ig=1,ngrid
5629	if (l.le.lmaxa(ig)) then
5630	zw=max(wa_moy(ig,l),1.e-10)
5631	larg_cons(ig,l)=zmax(ig)*r_aspect
5632	s fmc(ig,l)/(rhobarz(ig,l)zw)
5633	endif
5634	enddo
5635	enddo
5636
5637	do l=2,nlay
5638	do ig=1,ngrid
5639	if (l.le.lmaxa(ig)) then
5640	c if (idetr.eq.0) then
5641	c cette option est finalement en dur.
5642	if ((l_mix*zlev(ig,l)).lt.0.)then
5643	c print,'pb l_mixzlev<0'
5644	endif
5645	cCR: test: nouvelle def de lambda
5646	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5647	if (zw2(ig,l).gt.1.e-10) then
5648	larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
5649	else
5650	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5651	endif
5652	cRC
5653	c else if (idetr.eq.1) then
5654	c larg_detr(ig,l)=larg_cons(ig,l)
5655	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
5656	c else if (idetr.eq.2) then
5657	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5658	c s *sqrt(wa_moy(ig,l))
5659	c else if (idetr.eq.4) then
5660	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5661	c s *wa_moy(ig,l)
5662	c endif
5663	endif
5664	enddo
5665	enddo
5666
5667	c print*,'10 OK convect8'
5668	c print*,'WA2 ',wa_moy
5669	c calcul de la fraction de la maille concernée par l'ascendance en tenant
5670	c compte de l'epluchage du thermique.
5671	c
5672	cCR def de zmix continu (profil parabolique des vitesses)
5673	do ig=1,ngrid
5674	if (lmix(ig).gt.1.) then
5675	c test
5676	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
5677	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
5678	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
5679	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
5680	s then
5681	c
5682	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
5683	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
5684	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
5685	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
5686	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
5687	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
5688	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
5689	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
5690	else
5691	zmix(ig)=zlev(ig,lmix(ig))
5692	c print*,'pb zmix'
5693	endif
5694	else
5695	zmix(ig)=0.
5696	endif
5697	ctest
5698	if ((zmax(ig)-zmix(ig)).lt.0.) then
5699	zmix(ig)=0.99*zmax(ig)
5700	c print*,'pb zmix>zmax'
5701	endif
5702	enddo
5703	c
5704	c calcul du nouveau lmix correspondant
5705	do ig=1,ngrid
5706	do l=1,klev
5707	if (zmix(ig).ge.zlev(ig,l).and.
5708	s zmix(ig).lt.zlev(ig,l+1)) then
5709	lmix(ig)=l
5710	endif
5711	enddo
5712	enddo
5713	c
5714	do l=2,nlay
5715	do ig=1,ngrid
5716	if(larg_cons(ig,l).gt.1.) then
5717	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
5718	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
5719	s /(r_aspect*zmax(ig))
5720	c test
5721	fraca(ig,l)=max(fraca(ig,l),0.)
5722	fraca(ig,l)=min(fraca(ig,l),0.5)
5723	fracd(ig,l)=1.-fraca(ig,l)
5724	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
5725	else
5726	c wa_moy(ig,l)=0.
5727	fraca(ig,l)=0.
5728	fracc(ig,l)=0.
5729	fracd(ig,l)=1.
5730	endif
5731	enddo
5732	enddo
5733	cCR: calcul de fracazmix
5734	do ig=1,ngrid
5735	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
5736	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
5737	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
5738	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
5739	enddo
5740	c
5741	do l=2,nlay
5742	do ig=1,ngrid
5743	if(larg_cons(ig,l).gt.1.) then
5744	if (l.gt.lmix(ig)) then
5745	ctest
5746	if (zmax(ig)-zmix(ig).lt.1.e-10) then
5747	c print*,'pb xxx'
5748	xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
5749	else
5750	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
5751	endif
5752	if (idetr.eq.0) then
5753	fraca(ig,l)=fracazmix(ig)
5754	else if (idetr.eq.1) then
5755	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
5756	else if (idetr.eq.2) then
5757	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
5758	else
5759	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
5760	endif
5761	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
5762	fraca(ig,l)=max(fraca(ig,l),0.)
5763	fraca(ig,l)=min(fraca(ig,l),0.5)
5764	fracd(ig,l)=1.-fraca(ig,l)
5765	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
5766	endif
5767	endif
5768	enddo
5769	enddo
5770
5771	c print*,'fin calcul fraca'
5772	c print*,'11 OK convect8'
5773	c print*,'Ea3 ',wa_moy
5774	c------------------------------------------------------------------
5775	c Calcul de fracd, wd
5776	c somme wa - wd = 0
5777	c------------------------------------------------------------------
5778
5779
5780	do ig=1,ngrid
5781	fm(ig,1)=0.
5782	fm(ig,nlay+1)=0.
5783	enddo
5784
5785	do l=2,nlay
5786	do ig=1,ngrid
5787	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
5788	cCR:test
5789	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
5790	s .and.l.gt.lmix(ig)) then
5791	fm(ig,l)=fm(ig,l-1)
5792	c write(1,*)'ajustement fm, l',l
5793	endif
5794	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
5795	cRC
5796	enddo
5797	do ig=1,ngrid
5798	if(fracd(ig,l).lt.0.1) then
5799	stop'fracd trop petit'
5800	else
5801	c vitesse descendante "diagnostique"
5802	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
5803	endif
5804	enddo
5805	enddo
5806
5807	do l=1,nlay
5808	do ig=1,ngrid
5809	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
5810	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
5811	enddo
5812	enddo
5813
5814	c print*,'12 OK convect8'
5815	c print*,'WA4 ',wa_moy
5816	cc------------------------------------------------------------------
5817	c calcul du transport vertical
5818	c------------------------------------------------------------------
5819
5820	go to 4444
5821	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
5822	do l=2,nlay-1
5823	do ig=1,ngrid
5824	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
5825	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
5826	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
5827	c s ,fm(ig,l+1)*ptimestep
5828	c s ,' M=',masse(ig,l),masse(ig,l+1)
5829	endif
5830	enddo
5831	enddo
5832
5833	do l=1,nlay
5834	do ig=1,ngrid
5835	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
5836	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
5837	c s ,entr(ig,l)*ptimestep
5838	c s ,' M=',masse(ig,l)
5839	endif
5840	enddo
5841	enddo
5842
5843	do l=1,nlay
5844	do ig=1,ngrid
5845	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
5846	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
5847	c s ,' FM=',fm(ig,l)
5848	endif
5849	if(.not.masse(ig,l).ge.1.e-10
5850	s .or..not.masse(ig,l).le.1.e4) then
5851	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
5852	c s ,' M=',masse(ig,l)
5853	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
5854	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
5855	c print*,'zlev(ig,l+1),zlev(ig,l)'
5856	c s ,zlev(ig,l+1),zlev(ig,l)
5857	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
5858	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
5859	endif
5860	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
5861	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
5862	c s ,' E=',entr(ig,l)
5863	endif
5864	enddo
5865	enddo
5866
5867	4444 continue
5868
5869	cCR:redefinition du entr
5870	do l=1,nlay
5871	do ig=1,ngrid
5872	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
5873	if (detr(ig,l).lt.0.) then
5874	entr(ig,l)=entr(ig,l)-detr(ig,l)
5875	detr(ig,l)=0.
5876	c print*,'WARNING !!! detrainement negatif ',ig,l
5877	endif
5878	enddo
5879	enddo
5880	cRC
5881	if (w2di.eq.1) then
5882	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
5883	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
5884	else
5885	fm0=fm
5886	entr0=entr
5887	endif
5888
5889	if (1.eq.1) then
5890	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5891	. ,zh,zdhadj,zha)
5892	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5893	. ,zo,pdoadj,zoa)
5894	else
5895	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
5896	. ,zh,zdhadj,zha)
5897	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
5898	. ,zo,pdoadj,zoa)
5899	endif
5900
5901	if (1.eq.0) then
5902	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
5903	. ,fraca,zmax
5904	. ,zu,zv,pduadj,pdvadj,zua,zva)
5905	else
5906	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5907	. ,zu,pduadj,zua)
5908	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
5909	. ,zv,pdvadj,zva)
5910	endif
5911
5912	do l=1,nlay
5913	do ig=1,ngrid
5914	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
5915	zf2=zf/(1.-zf)
5916	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
5917	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
5918	enddo
5919	enddo
5920
5921
5922
5923	c print*,'13 OK convect8'
5924	c print*,'WA5 ',wa_moy
5925	do l=1,nlay
5926	do ig=1,ngrid
5927	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
5928	enddo
5929	enddo
5930
5931
5932	c do l=1,nlay
5933	c do ig=1,ngrid
5934	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
5935	c print*,'WARN!!! ig=',ig,' l=',l
5936	c s ,' pdtadj=',pdtadj(ig,l)
5937	c endif
5938	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
5939	c print*,'WARN!!! ig=',ig,' l=',l
5940	c s ,' pdoadj=',pdoadj(ig,l)
5941	c endif
5942	c enddo
5943	c enddo
5944
5945	c print*,'14 OK convect8'
5946	c------------------------------------------------------------------
5947	c Calculs pour les sorties
5948	c------------------------------------------------------------------
5949
5950	if(sorties) then
5951	do l=1,nlay
5952	do ig=1,ngrid
5953	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
5954	zld(ig,l)=fracd(ig,l)*zmax(ig)
5955	if(1.-fracd(ig,l).gt.1.e-10)
5956	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
5957	enddo
5958	enddo
5959
5960	cdeja fait
5961	c do l=1,nlay
5962	c do ig=1,ngrid
5963	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
5964	c if (detr(ig,l).lt.0.) then
5965	c entr(ig,l)=entr(ig,l)-detr(ig,l)
5966	c detr(ig,l)=0.
5967	c print*,'WARNING !!! detrainement negatif ',ig,l
5968	c endif
5969	c enddo
5970	c enddo
5971
5972	c print*,'15 OK convect8'
5973
5974	isplit=isplit+1
5975
5976
5977	c #define und
5978	goto 123
5979	#ifdef und
5980	CALL writeg1d(1,nlay,wd,'wd ','wd ')
5981	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
5982	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
5983	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
5984	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
5985	CALL writeg1d(1,nlay,zla,'la ','la ')
5986	CALL writeg1d(1,nlay,zld,'ld ','ld ')
5987	CALL writeg1d(1,nlay,pt,'pt ','pt ')
5988	CALL writeg1d(1,nlay,zh,'zh ','zh ')
5989	CALL writeg1d(1,nlay,zha,'zha ','zha ')
5990	CALL writeg1d(1,nlay,zu,'zu ','zu ')
5991	CALL writeg1d(1,nlay,zv,'zv ','zv ')
5992	CALL writeg1d(1,nlay,zo,'zo ','zo ')
5993	CALL writeg1d(1,nlay,wh,'wh ','wh ')
5994	CALL writeg1d(1,nlay,wu,'wu ','wu ')
5995	CALL writeg1d(1,nlay,wv,'wv ','wv ')
5996	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
5997	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
5998	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
5999	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
6000	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
6001	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
6002	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
6003	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
6004
6005	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
6006	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
6007	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
6008
6009	c recalcul des flux en diagnostique...
6010	c print*,'PAS DE TEMPS ',ptimestep
6011	call dt2F(pplev,pplay,pt,pdtadj,wh)
6012	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
6013	#endif
6014	123 continue
6015	! #define troisD
6016	#ifdef troisD
6017	c if (sorties) then
6018	print*,'Debut des wrgradsfi'
6019
6020	c print*,'16 OK convect8'
6021	call wrgradsfi(1,nlay,wd,'wd ','wd ')
6022	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
6023	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
6024	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
6025	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
6026	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
6027	c print*,'WA6 ',wa_moy
6028	call wrgradsfi(1,nlay,zla,'la ','la ')
6029	call wrgradsfi(1,nlay,zld,'ld ','ld ')
6030	call wrgradsfi(1,nlay,pt,'pt ','pt ')
6031	call wrgradsfi(1,nlay,zh,'zh ','zh ')
6032	call wrgradsfi(1,nlay,zha,'zha ','zha ')
6033	call wrgradsfi(1,nlay,zua,'zua ','zua ')
6034	call wrgradsfi(1,nlay,zva,'zva ','zva ')
6035	call wrgradsfi(1,nlay,zu,'zu ','zu ')
6036	call wrgradsfi(1,nlay,zv,'zv ','zv ')
6037	call wrgradsfi(1,nlay,zo,'zo ','zo ')
6038	call wrgradsfi(1,nlay,wh,'wh ','wh ')
6039	call wrgradsfi(1,nlay,wu,'wu ','wu ')
6040	call wrgradsfi(1,nlay,wv,'wv ','wv ')
6041	call wrgradsfi(1,nlay,wo,'wo ','wo ')
6042	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
6043	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
6044	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
6045	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
6046	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
6047	call wrgradsfi(1,nlay,entr,'entr ','entr ')
6048	call wrgradsfi(1,nlay,detr,'detr ','detr ')
6049	call wrgradsfi(1,nlay,fm,'fm ','fm ')
6050	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
6051	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
6052	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
6053	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
6054
6055	call wrgradsfi(1,nlay,zo,'zo ','zo ')
6056	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
6057	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
6058
6059	cCR:nouveaux diagnostiques
6060	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
6061	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
6062	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
6063	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
6064	zsortie1d(:)=lmax(:)
6065	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
6066	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
6067	zsortie1d(:)=lmix(:)
6068	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
6069	zsortie1d(:)=lentr(:)
6070	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
6071
6072	c print*,'17 OK convect8'
6073
6074	do k=1,klev/10
6075	write(str2,'(i2.2)') k
6076	str10='wa'//str2
6077	do l=1,nlay
6078	do ig=1,ngrid
6079	zsortie(ig,l)=wa(ig,k,l)
6080	enddo
6081	enddo
6082	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
6083	do l=1,nlay
6084	do ig=1,ngrid
6085	zsortie(ig,l)=larg_part(ig,k,l)
6086	enddo
6087	enddo
6088	str10='la'//str2
6089	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
6090	enddo
6091
6092
6093	c print*,'18 OK convect8'
6094	c endif
6095	print*,'Fin des wrgradsfi'
6096	#endif
6097
6098	endif
6099
6100	c if(wa_moy(1,4).gt.1.e-10) stop
6101
6102	c print*,'19 OK convect8'
6103	return
6104	end
6105

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

Download in other formats:

Original Format