Context Navigation

thermcell_old.F @ 1978

Last change on this file since 1978 was 1978, checked in by fhourdin, 11 years ago
Suppression de l'appel à wrgradsfi Removing calls to wrgradsfi
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 162.9 KB

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

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

Download in other formats:

Original Format