Context Navigation

thermcell_old.F @ 5440

Last change on this file since 5440 was 987, checked in by Laurent Fairhead, 16 years ago
Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9 et ORCHIDEE YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 189.4 KB

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

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

Download in other formats:

Original Format