Context Navigation

thermcell_old.F @ 901

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

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

Download in other formats:

Original Format