Context Navigation

thermcell_old.F90 @ 3819

Last change on this file since 3819 was 3819, checked in by ymipsl, 10 years ago

Removed all iim et jjm depedency. Replaced by nbp_lon and nbp_lat.
Supress gr_fi_ecrit, replaced by grid1dTo2d_glo

File size: 149.6 KB

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

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

Download in other formats:

Original Format