Context Navigation

thermcell_old.F @ 1

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

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

Download in other formats:

Original Format