Context Navigation

lmdz_thermcell_old.f90

Last change on this file was 5434, checked in by fhourdin, 15 hours ago

Superessing CPP in lmdz_*

Not possible for lmdz_thermcell_main because of isotopes

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

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

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

Download in other formats:

Original Format