Context Navigation

thermcell_old.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 190.2 KB

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

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

Download in other formats:

Original Format