Context Navigation

thermcell_old.F90 @ 3654

Last change on this file since 3654 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 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	! print*,'ZMAX ZMAX ZMAX ',zmax
330	! call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX ')
331
332	! print*,'OKl336'
333	! Calcul de l'entrainement.
334	! Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
335	! de la couche d'alimentation en partant du principe que la vitesse
336	! maximum dans l'ascendance est la vitesse d'entrainement horizontale.
337	DO k = 1, nlay
338	DO ig = 1, ngrid
339	zzz = rho(ig, k)wmax(ig, k)(zlev(ig,k+1)-zlev(ig,k))/ &
340	(zmax(ig)*r_aspect)
341	IF (w2di==2) THEN
342	entr(ig, k) = entr(ig, k) + ptimestep*(zzz-entr(ig,k))/tho
343	ELSE
344	entr(ig, k) = zzz
345	END IF
346	ztva(ig, k) = ztv(ig, k)
347	END DO
348	END DO
349
350
351	! print*,'7 OK convect8'
352	DO k = 1, klev + 1
353	DO ig = 1, ngrid
354	zw2(ig, k) = 0.
355	fmc(ig, k) = 0.
356	larg_cons(ig, k) = 0.
357	larg_detr(ig, k) = 0.
358	wa_moy(ig, k) = 0.
359	END DO
360	END DO
361
362	! print*,'8 OK convect8'
363	DO ig = 1, ngrid
364	lmaxa(ig) = 1
365	lmix(ig) = 1
366	wmaxa(ig) = 0.
367	END DO
368
369
370	! print*,'OKl372'
371	DO l = 1, nlay - 2
372	DO ig = 1, ngrid
373	! if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)) then
374	! print*,'COUCOU ',l,zw2(ig,l),ztv(ig,l),ztv(ig,l+1)
375	IF (zw2(ig,l)<1.E-10 .AND. ztv(ig,l)>ztv(ig,l+1) .AND. &
376	entr(ig,l)>1.E-10) THEN
377	! print*,'COUCOU cas 1'
378	! Initialisation de l'ascendance
379	! lmix(ig)=1
380	ztva(ig, l) = ztv(ig, l)
381	fmc(ig, l) = 0.
382	fmc(ig, l+1) = entr(ig, l)
383	zw2(ig, l) = 0.
384	! if (.not.ztv(ig,l+1).gt.150.) then
385	! print*,'ig,l+1,ztv(ig,l+1)'
386	! print*, ig,l+1,ztv(ig,l+1)
387	! endif
388	zw2(ig, l+1) = 2.rg(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
389	(zlev(ig,l+1)-zlev(ig,l))
390	larg_detr(ig, l) = 0.
391	ELSE IF (zw2(ig,l)>=1.E-10 .AND. fmc(ig,l)+entr(ig,l)>1.E-10) THEN
392	! Incrementation...
393	fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
394	! if (.not.fmc(ig,l+1).gt.1.e-15) then
395	! print*,'ig,l+1,fmc(ig,l+1)'
396	! print*, ig,l+1,fmc(ig,l+1)
397	! print*,'Fmc ',(fmc(ig,ll),ll=1,klev+1)
398	! print*,'W2 ',(zw2(ig,ll),ll=1,klev+1)
399	! print*,'Tv ',(ztv(ig,ll),ll=1,klev)
400	! print*,'Entr ',(entr(ig,ll),ll=1,klev)
401	! endif
402	ztva(ig, l) = (fmc(ig,l)ztva(ig,l-1)+entr(ig,l)ztv(ig,l))/ &
403	fmc(ig, l+1)
404	! mise a jour de la vitesse ascendante (l'air entraine de la couche
405	! consideree commence avec une vitesse nulle).
406	zw2(ig, l+1) = zw2(ig, l)(fmc(ig,l)/fmc(ig,l+1))*2 + &
407	2.rg(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
408	END IF
409	IF (zw2(ig,l+1)<0.) THEN
410	zw2(ig, l+1) = 0.
411	lmaxa(ig) = l
412	ELSE
413	wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
414	END IF
415	IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
416	! lmix est le niveau de la couche ou w (wa_moy) est maximum
417	lmix(ig) = l + 1
418	wmaxa(ig) = wa_moy(ig, l+1)
419	END IF
420	! print*,'COUCOU cas 2 LMIX=',lmix(ig),wa_moy(ig,l+1),wmaxa(ig)
421	END DO
422	END DO
423
424	! print*,'9 OK convect8'
425	! print*,'WA1 ',wa_moy
426
427	! determination de l'indice du debut de la mixed layer ou w decroit
428
429	! calcul de la largeur de chaque ascendance dans le cas conservatif.
430	! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
431	! d'une couche est égale à la hauteur de la couche alimentante.
432	! La vitesse maximale dans l'ascendance est aussi prise comme estimation
433	! de la vitesse d'entrainement horizontal dans la couche alimentante.
434
435	! print*,'OKl439'
436	DO l = 2, nlay
437	DO ig = 1, ngrid
438	IF (l<=lmaxa(ig)) THEN
439	zw = max(wa_moy(ig,l), 1.E-10)
440	larg_cons(ig, l) = zmax(ig)r_aspectfmc(ig, l)/(rhobarz(ig,l)*zw)
441	END IF
442	END DO
443	END DO
444
445	DO l = 2, nlay
446	DO ig = 1, ngrid
447	IF (l<=lmaxa(ig)) THEN
448	! if (idetr.eq.0) then
449	! cette option est finalement en dur.
450	larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
451	! else if (idetr.eq.1) then
452	! larg_detr(ig,l)=larg_cons(ig,l)
453	! s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
454	! else if (idetr.eq.2) then
455	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
456	! s *sqrt(wa_moy(ig,l))
457	! else if (idetr.eq.4) then
458	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
459	! s *wa_moy(ig,l)
460	! endif
461	END IF
462	END DO
463	END DO
464
465	! print*,'10 OK convect8'
466	! print*,'WA2 ',wa_moy
467	! calcul de la fraction de la maille concernée par l'ascendance en tenant
468	! compte de l'epluchage du thermique.
469
470	DO l = 2, nlay
471	DO ig = 1, ngrid
472	IF (larg_cons(ig,l)>1.) THEN
473	! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
474	fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
475	IF (l>lmix(ig)) THEN
476	xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
477	IF (idetr==0) THEN
478	fraca(ig, l) = fraca(ig, lmix(ig))
479	ELSE IF (idetr==1) THEN
480	fraca(ig, l) = fraca(ig, lmix(ig))*xxx(ig, l)
481	ELSE IF (idetr==2) THEN
482	fraca(ig, l) = fraca(ig, lmix(ig))(1.-(1.-xxx(ig,l))*2)
483	ELSE
484	fraca(ig, l) = fraca(ig, lmix(ig))xxx(ig, l)*2
485	END IF
486	END IF
487	! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
488	fraca(ig, l) = max(fraca(ig,l), 0.)
489	fraca(ig, l) = min(fraca(ig,l), 0.5)
490	fracd(ig, l) = 1. - fraca(ig, l)
491	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
492	ELSE
493	! wa_moy(ig,l)=0.
494	fraca(ig, l) = 0.
495	fracc(ig, l) = 0.
496	fracd(ig, l) = 1.
497	END IF
498	END DO
499	END DO
500
501	! print*,'11 OK convect8'
502	! print*,'Ea3 ',wa_moy
503	! ------------------------------------------------------------------
504	! Calcul de fracd, wd
505	! somme wa - wd = 0
506	! ------------------------------------------------------------------
507
508
509	DO ig = 1, ngrid
510	fm(ig, 1) = 0.
511	fm(ig, nlay+1) = 0.
512	END DO
513
514	DO l = 2, nlay
515	DO ig = 1, ngrid
516	fm(ig, l) = fraca(ig, l)wa_moy(ig, l)rhobarz(ig, l)
517	END DO
518	DO ig = 1, ngrid
519	IF (fracd(ig,l)<0.1) THEN
520	abort_message = 'fracd trop petit'
521	CALL abort_physic(modname, abort_message, 1)
522	ELSE
523	! vitesse descendante "diagnostique"
524	wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
525	END IF
526	END DO
527	END DO
528
529	DO l = 1, nlay
530	DO ig = 1, ngrid
531	! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
532	masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
533	END DO
534	END DO
535
536	! print*,'12 OK convect8'
537	! print*,'WA4 ',wa_moy
538	! c------------------------------------------------------------------
539	! calcul du transport vertical
540	! ------------------------------------------------------------------
541
542	GO TO 4444
543	! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
544	DO l = 2, nlay - 1
545	DO ig = 1, ngrid
546	IF (fm(ig,l+1)ptimestep>masse(ig,l) .AND. fm(ig,l+1)ptimestep>masse( &
547	ig,l+1)) THEN
548	! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
549	! s ,fm(ig,l+1)*ptimestep
550	! s ,' M=',masse(ig,l),masse(ig,l+1)
551	END IF
552	END DO
553	END DO
554
555	DO l = 1, nlay
556	DO ig = 1, ngrid
557	IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
558	! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
559	! s ,entr(ig,l)*ptimestep
560	! s ,' M=',masse(ig,l)
561	END IF
562	END DO
563	END DO
564
565	DO l = 1, nlay
566	DO ig = 1, ngrid
567	IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
568	! print*,'WARN!!! fm exagere ig=',ig,' l=',l
569	! s ,' FM=',fm(ig,l)
570	END IF
571	IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
572	! print*,'WARN!!! masse exagere ig=',ig,' l=',l
573	! s ,' M=',masse(ig,l)
574	! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
575	! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
576	! print*,'zlev(ig,l+1),zlev(ig,l)'
577	! s ,zlev(ig,l+1),zlev(ig,l)
578	! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
579	! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
580	END IF
581	IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
582	! print*,'WARN!!! entr exagere ig=',ig,' l=',l
583	! s ,' E=',entr(ig,l)
584	END IF
585	END DO
586	END DO
587
588	4444 CONTINUE
589	! print*,'OK 444 '
590
591	IF (w2di==1) THEN
592	fm0 = fm0 + ptimestep*(fm-fm0)/tho
593	entr0 = entr0 + ptimestep*(entr-entr0)/tho
594	ELSE
595	fm0 = fm
596	entr0 = entr
597	END IF
598
599	IF (flagdq==0) THEN
600	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
601	zha)
602	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
603	zoa)
604	PRINT *, 'THERMALS OPT 1'
605	ELSE IF (flagdq==1) THEN
606	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
607	zdhadj, zha)
608	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
609	pdoadj, zoa)
610	PRINT *, 'THERMALS OPT 2'
611	ELSE
612	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zh, &
613	zdhadj, zha, lev_out)
614	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zo, &
615	pdoadj, zoa, lev_out)
616	PRINT *, 'THERMALS OPT 3', dqimpl
617	END IF
618
619	PRINT *, 'TH VENT ', dvdq
620	IF (dvdq==0) THEN
621	! print*,'TH VENT OK ',dvdq
622	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
623	zua)
624	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
625	zva)
626	ELSE IF (dvdq==1) THEN
627	CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
628	zu, zv, pduadj, pdvadj, zua, zva)
629	ELSE IF (dvdq==2) THEN
630	CALL thermcell_dv2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, &
631	zmax, zu, zv, pduadj, pdvadj, zua, zva, lev_out)
632	ELSE IF (dvdq==3) THEN
633	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zu, &
634	pduadj, zua, lev_out)
635	CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zv, &
636	pdvadj, zva, lev_out)
637	END IF
638
639	! CALL writefield_phy('duadj',pduadj,klev)
640
641	DO l = 1, nlay
642	DO ig = 1, ngrid
643	zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
644	zf2 = zf/(1.-zf)
645	thetath2(ig, l) = zf2(zha(ig,l)-zh(ig,l))*2
646	wth2(ig, l) = zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
647	END DO
648	END DO
649
650
651
652	! print*,'13 OK convect8'
653	! print*,'WA5 ',wa_moy
654	DO l = 1, nlay
655	DO ig = 1, ngrid
656	pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
657	END DO
658	END DO
659
660
661	! do l=1,nlay
662	! do ig=1,ngrid
663	! if(abs(pdtadj(ig,l))*86400..gt.500.) then
664	! print*,'WARN!!! ig=',ig,' l=',l
665	! s ,' pdtadj=',pdtadj(ig,l)
666	! endif
667	! if(abs(pdoadj(ig,l))*86400..gt.1.) then
668	! print*,'WARN!!! ig=',ig,' l=',l
669	! s ,' pdoadj=',pdoadj(ig,l)
670	! endif
671	! enddo
672	! enddo
673
674	! print*,'14 OK convect8'
675	! ------------------------------------------------------------------
676	! Calculs pour les sorties
677	! ------------------------------------------------------------------
678
679	IF (sorties) THEN
680	DO l = 1, nlay
681	DO ig = 1, ngrid
682	zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
683	zld(ig, l) = fracd(ig, l)*zmax(ig)
684	IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
685	(1.-fracd(ig,l))
686	END DO
687	END DO
688
689	DO l = 1, nlay
690	DO ig = 1, ngrid
691	detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
692	IF (detr(ig,l)<0.) THEN
693	entr(ig, l) = entr(ig, l) - detr(ig, l)
694	detr(ig, l) = 0.
695	! print*,'WARNING !!! detrainement negatif ',ig,l
696	END IF
697	END DO
698	END DO
699	END IF
700
701	! print*,'15 OK convect8'
702
703
704	! if(wa_moy(1,4).gt.1.e-10) stop
705
706	! print*,'19 OK convect8'
707	RETURN
708	END SUBROUTINE thermcell_2002
709
710	SUBROUTINE thermcell_cld(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, &
711	debut, pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0, zqla, &
712	lmax, zmax_sec, wmax_sec, zw_sec, lmix_sec, ratqscth, ratqsdiff & ! s
713	! ,pu_therm,pv_therm
714	, r_aspect, l_mix, w2di, tho)
715
716	USE dimphy
717	IMPLICIT NONE
718
719	! =======================================================================
720
721	! Calcul du transport verticale dans la couche limite en presence
722	! de "thermiques" explicitement representes
723
724	! Réécriture à partir d'un listing papier à Habas, le 14/02/00
725
726	! le thermique est supposé homogène et dissipé par mélange avec
727	! son environnement. la longueur l_mix contrôle l'efficacité du
728	! mélange
729
730	! Le calcul du transport des différentes espèces se fait en prenant
731	! en compte:
732	! 1. un flux de masse montant
733	! 2. un flux de masse descendant
734	! 3. un entrainement
735	! 4. un detrainement
736
737	! =======================================================================
738
739	! -----------------------------------------------------------------------
740	! declarations:
741	! -------------
742
743	include "YOMCST.h"
744	include "YOETHF.h"
745	include "FCTTRE.h"
746
747	! arguments:
748	! ----------
749
750	INTEGER ngrid, nlay, w2di
751	REAL tho
752	REAL ptimestep, l_mix, r_aspect
753	REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
754	REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
755	REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
756	REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
757	REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
758	REAL pphi(ngrid, nlay)
759
760	INTEGER idetr
761	SAVE idetr
762	DATA idetr/3/
763	!$OMP THREADPRIVATE(idetr)
764
765	! local:
766	! ------
767
768	INTEGER ig, k, l, lmaxa(klon), lmix(klon)
769	REAL zsortie1d(klon)
770	! CR: on remplace lmax(klon,klev+1)
771	INTEGER lmax(klon), lmin(klon), lentr(klon)
772	REAL linter(klon)
773	REAL zmix(klon), fracazmix(klon)
774	REAL alpha
775	SAVE alpha
776	DATA alpha/1./
777	!$OMP THREADPRIVATE(alpha)
778
779	! RC
780	REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
781	REAL zmax_sec(klon)
782	REAL zmax_sec2(klon)
783	REAL zw_sec(klon, klev+1)
784	INTEGER lmix_sec(klon)
785	REAL w_est(klon, klev+1)
786	! on garde le zmax du pas de temps precedent
787	! real zmax0(klon)
788	! save zmax0
789	! real zmix0(klon)
790	! save zmix0
791	REAL, SAVE, ALLOCATABLE :: zmax0(:), zmix0(:)
792	!$OMP THREADPRIVATE(zmax0, zmix0)
793
794	REAL zlev(klon, klev+1), zlay(klon, klev)
795	REAL deltaz(klon, klev)
796	REAL zh(klon, klev), zdhadj(klon, klev)
797	REAL zthl(klon, klev), zdthladj(klon, klev)
798	REAL ztv(klon, klev)
799	REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
800	REAL zl(klon, klev)
801	REAL wh(klon, klev+1)
802	REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
803	REAL zla(klon, klev+1)
804	REAL zwa(klon, klev+1)
805	REAL zld(klon, klev+1)
806	REAL zwd(klon, klev+1)
807	REAL zsortie(klon, klev)
808	REAL zva(klon, klev)
809	REAL zua(klon, klev)
810	REAL zoa(klon, klev)
811
812	REAL zta(klon, klev)
813	REAL zha(klon, klev)
814	REAL wa_moy(klon, klev+1)
815	REAL fraca(klon, klev+1)
816	REAL fracc(klon, klev+1)
817	REAL zf, zf2
818	REAL thetath2(klon, klev), wth2(klon, klev), wth3(klon, klev)
819	REAL q2(klon, klev)
820	REAL dtheta(klon, klev)
821	! common/comtherm/thetath2,wth2
822
823	REAL ratqscth(klon, klev)
824	REAL sum
825	REAL sumdiff
826	REAL ratqsdiff(klon, klev)
827	REAL count_time
828	INTEGER ialt
829
830	LOGICAL sorties
831	REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
832	REAL zpspsk(klon, klev)
833
834	! real wmax(klon,klev),wmaxa(klon)
835	REAL wmax(klon), wmaxa(klon)
836	REAL wmax_sec(klon)
837	REAL wmax_sec2(klon)
838	REAL wa(klon, klev, klev+1)
839	REAL wd(klon, klev+1)
840	REAL larg_part(klon, klev, klev+1)
841	REAL fracd(klon, klev+1)
842	REAL xxx(klon, klev+1)
843	REAL larg_cons(klon, klev+1)
844	REAL larg_detr(klon, klev+1)
845	REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
846	REAL massetot(klon, klev)
847	REAL detr0(klon, klev)
848	REAL alim0(klon, klev)
849	REAL pu_therm(klon, klev), pv_therm(klon, klev)
850	REAL fm(klon, klev+1), entr(klon, klev)
851	REAL fmc(klon, klev+1)
852
853	REAL zcor, zdelta, zcvm5, qlbef
854	REAL tbef(klon), qsatbef(klon)
855	REAL dqsat_dt, dt, num, denom
856	REAL reps, rlvcp, ddt0
857	REAL ztla(klon, klev), zqla(klon, klev), zqta(klon, klev)
858	! CR niveau de condensation
859	REAL nivcon(klon)
860	REAL zcon(klon)
861	REAL zqsat(klon, klev)
862	REAL zqsatth(klon, klev)
863	PARAMETER (ddt0=.01)
864
865
866	! CR:nouvelles variables
867	REAL f_star(klon, klev+1), entr_star(klon, klev)
868	REAL detr_star(klon, klev)
869	REAL alim_star_tot(klon), alim_star2(klon)
870	REAL entr_star_tot(klon)
871	REAL detr_star_tot(klon)
872	REAL alim_star(klon, klev)
873	REAL alim(klon, klev)
874	REAL nu(klon, klev)
875	REAL nu_e(klon, klev)
876	REAL nu_min
877	REAL nu_max
878	REAL nu_r
879	REAL f(klon)
880	! real f(klon), f0(klon)
881	! save f0
882	REAL, SAVE, ALLOCATABLE :: f0(:)
883	!$OMP THREADPRIVATE(f0)
884
885	REAL f_old
886	REAL zlevinter(klon)
887	LOGICAL, SAVE :: first = .TRUE.
888	!$OMP THREADPRIVATE(first)
889	! data first /.false./
890	! save first
891	LOGICAL nuage
892	! save nuage
893	LOGICAL boucle
894	LOGICAL therm
895	LOGICAL debut
896	LOGICAL rale
897	INTEGER test(klon)
898	INTEGER signe_zw2
899	! RC
900
901	CHARACTER *2 str2
902	CHARACTER *10 str10
903
904	CHARACTER (LEN=20) :: modname = 'thermcell_cld'
905	CHARACTER (LEN=80) :: abort_message
906
907	LOGICAL vtest(klon), down
908	LOGICAL zsat(klon)
909
910	EXTERNAL scopy
911
912	INTEGER ncorrec, ll
913	SAVE ncorrec
914	DATA ncorrec/0/
915	!$OMP THREADPRIVATE(ncorrec)
916
917
918
919	! -----------------------------------------------------------------------
920	! initialisation:
921	! ---------------
922
923	IF (first) THEN
924	ALLOCATE (zmix0(klon))
925	ALLOCATE (zmax0(klon))
926	ALLOCATE (f0(klon))
927	first = .FALSE.
928	END IF
929
930	sorties = .FALSE.
931	! print*,'NOUVEAU DETR PLUIE '
932	IF (ngrid/=klon) THEN
933	PRINT *
934	PRINT *, 'STOP dans convadj'
935	PRINT *, 'ngrid =', ngrid
936	PRINT *, 'klon =', klon
937	END IF
938
939	! Initialisation
940	rlvcp = rlvtt/rcpd
941	reps = rd/rv
942	! initialisations de zqsat
943	DO ll = 1, nlay
944	DO ig = 1, ngrid
945	zqsat(ig, ll) = 0.
946	zqsatth(ig, ll) = 0.
947	END DO
948	END DO
949
950	! on met le first a true pour le premier passage de la journée
951	DO ig = 1, klon
952	test(ig) = 0
953	END DO
954	IF (debut) THEN
955	DO ig = 1, klon
956	test(ig) = 1
957	f0(ig) = 0.
958	zmax0(ig) = 0.
959	END DO
960	END IF
961	DO ig = 1, klon
962	IF ((.NOT. debut) .AND. (f0(ig)<1.E-10)) THEN
963	test(ig) = 1
964	END IF
965	END DO
966	! do ig=1,klon
967	! print*,'test(ig)',test(ig),zmax0(ig)
968	! enddo
969	nuage = .FALSE.
970	! -----------------------------------------------------------------------
971	! AM Calcul de T,q,ql a partir de Tl et qT
972	! ---------------------------------------------------
973
974	! Pr Tprec=Tl calcul de qsat
975	! Si qsat>qT T=Tl, q=qT
976	! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
977	! On cherche DDT < DDT0
978
979	! defaut
980	DO ll = 1, nlay
981	DO ig = 1, ngrid
982	zo(ig, ll) = po(ig, ll)
983	zl(ig, ll) = 0.
984	zh(ig, ll) = pt(ig, ll)
985	END DO
986	END DO
987	DO ig = 1, ngrid
988	zsat(ig) = .FALSE.
989	END DO
990
991
992	DO ll = 1, nlay
993	! les points insatures sont definitifs
994	DO ig = 1, ngrid
995	tbef(ig) = pt(ig, ll)
996	zdelta = max(0., sign(1.,rtt-tbef(ig)))
997	qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
998	qsatbef(ig) = min(0.5, qsatbef(ig))
999	zcor = 1./(1.-retv*qsatbef(ig))
1000	qsatbef(ig) = qsatbef(ig)*zcor
1001	zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig))>1.E-10)
1002	END DO
1003
1004	DO ig = 1, ngrid
1005	IF (zsat(ig) .AND. (1==1)) THEN
1006	qlbef = max(0., po(ig,ll)-qsatbef(ig))
1007	! si sature: ql est surestime, d'ou la sous-relax
1008	dt = 0.5rlvcpqlbef
1009	! write(18,*),'DT0=',DT
1010	! on pourra enchainer 2 ou 3 calculs sans Do while
1011	DO WHILE (abs(dt)>ddt0)
1012	! il faut verifier si c,a conserve quand on repasse en insature ...
1013	tbef(ig) = tbef(ig) + dt
1014	zdelta = max(0., sign(1.,rtt-tbef(ig)))
1015	qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
1016	qsatbef(ig) = min(0.5, qsatbef(ig))
1017	zcor = 1./(1.-retv*qsatbef(ig))
1018	qsatbef(ig) = qsatbef(ig)*zcor
1019	! on veut le signe de qlbef
1020	qlbef = po(ig, ll) - qsatbef(ig)
1021	zdelta = max(0., sign(1.,rtt-tbef(ig)))
1022	zcvm5 = r5les(1.-zdelta) + r5ieszdelta
1023	zcor = 1./(1.-retv*qsatbef(ig))
1024	dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
1025	num = -tbef(ig) + pt(ig, ll) + rlvcp*qlbef
1026	denom = 1. + rlvcp*dqsat_dt
1027	IF (denom<1.E-10) THEN
1028	PRINT *, 'pb denom'
1029	END IF
1030	dt = num/denom
1031	END DO
1032	! on ecrit de maniere conservative (sat ou non)
1033	zl(ig, ll) = max(0., qlbef)
1034	! T = Tl +Lv/Cp ql
1035	zh(ig, ll) = pt(ig, ll) + rlvcp*zl(ig, ll)
1036	zo(ig, ll) = po(ig, ll) - zl(ig, ll)
1037	END IF
1038	! on ecrit zqsat
1039	zqsat(ig, ll) = qsatbef(ig)
1040	END DO
1041	END DO
1042	! AM fin
1043
1044	! -----------------------------------------------------------------------
1045	! incrementation eventuelle de tendances precedentes:
1046	! ---------------------------------------------------
1047
1048	! print*,'0 OK convect8'
1049
1050	DO l = 1, nlay
1051	DO ig = 1, ngrid
1052	zpspsk(ig, l) = (pplay(ig,l)/100000.)**rkappa
1053	! zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
1054	! zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
1055	zu(ig, l) = pu(ig, l)
1056	zv(ig, l) = pv(ig, l)
1057	! zo(ig,l)=po(ig,l)
1058	! ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
1059	! AM attention zh est maintenant le profil de T et plus le profil de
1060	! theta !
1061
1062	! T-> Theta
1063	ztv(ig, l) = zh(ig, l)/zpspsk(ig, l)
1064	! AM Theta_v
1065	ztv(ig, l) = ztv(ig, l)(1.+retv(zo(ig,l))-zl(ig,l))
1066	! AM Thetal
1067	zthl(ig, l) = pt(ig, l)/zpspsk(ig, l)
1068
1069	END DO
1070	END DO
1071
1072	! print*,'1 OK convect8'
1073	! --------------------
1074
1075
1076	! + + + + + + + + + + +
1077
1078
1079	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
1080	! wh,wt,wo ...
1081
1082	! + + + + + + + + + + + zh,zu,zv,zo,rho
1083
1084
1085	! -------------------- zlev(1)
1086	! \\\\\\\\\\\\\\\\\\\\
1087
1088
1089
1090	! -----------------------------------------------------------------------
1091	! Calcul des altitudes des couches
1092	! -----------------------------------------------------------------------
1093
1094	DO l = 2, nlay
1095	DO ig = 1, ngrid
1096	zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
1097	END DO
1098	END DO
1099	DO ig = 1, ngrid
1100	zlev(ig, 1) = 0.
1101	zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
1102	END DO
1103	DO l = 1, nlay
1104	DO ig = 1, ngrid
1105	zlay(ig, l) = pphi(ig, l)/rg
1106	END DO
1107	END DO
1108	! calcul de deltaz
1109	DO l = 1, nlay
1110	DO ig = 1, ngrid
1111	deltaz(ig, l) = zlev(ig, l+1) - zlev(ig, l)
1112	END DO
1113	END DO
1114
1115	! print*,'2 OK convect8'
1116	! -----------------------------------------------------------------------
1117	! Calcul des densites
1118	! -----------------------------------------------------------------------
1119
1120	DO l = 1, nlay
1121	DO ig = 1, ngrid
1122	! rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
1123	rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)rdztv(ig,l))
1124	END DO
1125	END DO
1126
1127	DO l = 2, nlay
1128	DO ig = 1, ngrid
1129	rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
1130	END DO
1131	END DO
1132
1133	DO k = 1, nlay
1134	DO l = 1, nlay + 1
1135	DO ig = 1, ngrid
1136	wa(ig, k, l) = 0.
1137	END DO
1138	END DO
1139	END DO
1140	! Cr:ajout:calcul de la masse
1141	DO l = 1, nlay
1142	DO ig = 1, ngrid
1143	! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1144	masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
1145	END DO
1146	END DO
1147	! print*,'3 OK convect8'
1148	! ------------------------------------------------------------------
1149	! Calcul de w2, quarre de w a partir de la cape
1150	! a partir de w2, on calcule wa, vitesse de l'ascendance
1151
1152	! ATTENTION: Dans cette version, pour cause d'economie de memoire,
1153	! w2 est stoke dans wa
1154
1155	! ATTENTION: dans convect8, on n'utilise le calcule des wa
1156	! independants par couches que pour calculer l'entrainement
1157	! a la base et la hauteur max de l'ascendance.
1158
1159	! Indicages:
1160	! l'ascendance provenant du niveau k traverse l'interface l avec
1161	! une vitesse wa(k,l).
1162
1163	! --------------------
1164
1165	! + + + + + + + + + +
1166
1167	! wa(k,l) ---- -------------------- l
1168	! /\
1169	! /\|\|\ + + + + + + + + + +
1170	! \|\|
1171	! \|\| --------------------
1172	! \|\|
1173	! \|\| + + + + + + + + + +
1174	! \|\|
1175	! \|\| --------------------
1176	! \|\|__
1177	! \|___ + + + + + + + + + + k
1178
1179	! --------------------
1180
1181
1182
1183	! ------------------------------------------------------------------
1184
1185	! CR: ponderation entrainement des couches instables
1186	! def des alim_star tels que alim=f*alim_star
1187	DO l = 1, klev
1188	DO ig = 1, ngrid
1189	alim_star(ig, l) = 0.
1190	alim(ig, l) = 0.
1191	END DO
1192	END DO
1193	! determination de la longueur de la couche d entrainement
1194	DO ig = 1, ngrid
1195	lentr(ig) = 1
1196	END DO
1197
1198	! on ne considere que les premieres couches instables
1199	therm = .FALSE.
1200	DO k = nlay - 2, 1, -1
1201	DO ig = 1, ngrid
1202	IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
1203	lentr(ig) = k + 1
1204	therm = .TRUE.
1205	END IF
1206	END DO
1207	END DO
1208
1209	! determination du lmin: couche d ou provient le thermique
1210	DO ig = 1, ngrid
1211	lmin(ig) = 1
1212	END DO
1213	DO ig = 1, ngrid
1214	DO l = nlay, 2, -1
1215	IF (ztv(ig,l-1)>ztv(ig,l)) THEN
1216	lmin(ig) = l - 1
1217	END IF
1218	END DO
1219	END DO
1220
1221	! definition de l'entrainement des couches
1222	DO l = 1, klev - 1
1223	DO ig = 1, ngrid
1224	IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<lentr(ig)) THEN
1225	! def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
1226	alim_star(ig, l) = max((ztv(ig,l)-ztv(ig,l+1)), 0.) & ! s
1227	! *(zlev(ig,l+1)-zlev(ig,l))
1228	*sqrt(zlev(ig,l+1))
1229	! alim_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
1230	! s /zlev(ig,lentr(ig)+2)))**(3./2.)
1231	END IF
1232	END DO
1233	END DO
1234
1235	! pas de thermique si couche 1 stable
1236	DO ig = 1, ngrid
1237	! if (lmin(ig).gt.1) then
1238	! CRnouveau test
1239	IF (alim_star(ig,1)<1.E-10) THEN
1240	DO l = 1, klev
1241	alim_star(ig, l) = 0.
1242	END DO
1243	END IF
1244	END DO
1245	! calcul de l entrainement total
1246	DO ig = 1, ngrid
1247	alim_star_tot(ig) = 0.
1248	entr_star_tot(ig) = 0.
1249	detr_star_tot(ig) = 0.
1250	END DO
1251	DO ig = 1, ngrid
1252	DO k = 1, klev
1253	alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, k)
1254	END DO
1255	END DO
1256
1257	! Calcul entrainement normalise
1258	DO ig = 1, ngrid
1259	IF (alim_star_tot(ig)>1.E-10) THEN
1260	! do l=1,lentr(ig)
1261	DO l = 1, klev
1262	! def possibles pour entr_star: zdthetadz, dthetadz, zdtheta
1263	alim_star(ig, l) = alim_star(ig, l)/alim_star_tot(ig)
1264	END DO
1265	END IF
1266	END DO
1267
1268	! print*,'fin calcul alim_star'
1269
1270	! AM:initialisations
1271	DO k = 1, nlay
1272	DO ig = 1, ngrid
1273	ztva(ig, k) = ztv(ig, k)
1274	ztla(ig, k) = zthl(ig, k)
1275	zqla(ig, k) = 0.
1276	zqta(ig, k) = po(ig, k)
1277	zsat(ig) = .FALSE.
1278	END DO
1279	END DO
1280	DO k = 1, klev
1281	DO ig = 1, ngrid
1282	detr_star(ig, k) = 0.
1283	entr_star(ig, k) = 0.
1284	detr(ig, k) = 0.
1285	entr(ig, k) = 0.
1286	END DO
1287	END DO
1288	! print*,'7 OK convect8'
1289	DO k = 1, klev + 1
1290	DO ig = 1, ngrid
1291	zw2(ig, k) = 0.
1292	fmc(ig, k) = 0.
1293	! CR
1294	f_star(ig, k) = 0.
1295	! RC
1296	larg_cons(ig, k) = 0.
1297	larg_detr(ig, k) = 0.
1298	wa_moy(ig, k) = 0.
1299	END DO
1300	END DO
1301
1302	! n print*,'8 OK convect8'
1303	DO ig = 1, ngrid
1304	linter(ig) = 1.
1305	lmaxa(ig) = 1
1306	lmix(ig) = 1
1307	wmaxa(ig) = 0.
1308	END DO
1309
1310	nu_min = l_mix
1311	nu_max = 1000.
1312	! do ig=1,ngrid
1313	! nu_max=wmax_sec(ig)
1314	! enddo
1315	DO ig = 1, ngrid
1316	DO k = 1, klev
1317	nu(ig, k) = 0.
1318	nu_e(ig, k) = 0.
1319	END DO
1320	END DO
1321	! Calcul de l'excès de température du à la diffusion turbulente
1322	DO ig = 1, ngrid
1323	DO l = 1, klev
1324	dtheta(ig, l) = 0.
1325	END DO
1326	END DO
1327	DO ig = 1, ngrid
1328	DO l = 1, lentr(ig) - 1
1329	dtheta(ig, l) = sqrt(10.0.4zlev(ig,l+1)*21.*((ztv(ig,l+1)- &
1330	ztv(ig,l))/(zlev(ig,l+1)-zlev(ig,l)))**2)
1331	END DO
1332	END DO
1333	! do l=1,nlay-2
1334	DO l = 1, klev - 1
1335	DO ig = 1, ngrid
1336	IF (ztv(ig,l)>ztv(ig,l+1) .AND. alim_star(ig,l)>1.E-10 .AND. &
1337	zw2(ig,l)<1E-10) THEN
1338	! AM
1339	! test:on rajoute un excès de T dans couche alim
1340	! ztla(ig,l)=zthl(ig,l)+dtheta(ig,l)
1341	ztla(ig, l) = zthl(ig, l)
1342	! test: on rajoute un excès de q dans la couche alim
1343	! zqta(ig,l)=po(ig,l)+0.001
1344	zqta(ig, l) = po(ig, l)
1345	zqla(ig, l) = zl(ig, l)
1346	! AM
1347	f_star(ig, l+1) = alim_star(ig, l)
1348	! test:calcul de dteta
1349	zw2(ig, l+1) = 2.rg(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
1350	(zlev(ig,l+1)-zlev(ig,l))0.4pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
1351	w_est(ig, l+1) = zw2(ig, l+1)
1352	larg_detr(ig, l) = 0.
1353	! print*,'coucou boucle 1'
1354	ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+alim_star(ig, &
1355	l))>1.E-10) THEN
1356	! print*,'coucou boucle 2'
1357	! estimation du detrainement a partir de la geometrie du pas
1358	! precedent
1359	IF ((test(ig)==1) .OR. ((.NOT. debut) .AND. (f0(ig)<1.E-10))) THEN
1360	detr_star(ig, l) = 0.
1361	entr_star(ig, l) = 0.
1362	! print*,'coucou test(ig)',test(ig),f0(ig),zmax0(ig)
1363	ELSE
1364	! print*,'coucou debut detr'
1365	! tests sur la definition du detr
1366	IF (zqla(ig,l-1)>1.E-10) THEN
1367	nuage = .TRUE.
1368	END IF
1369
1370	w_est(ig, l+1) = zw2(ig, l)((f_star(ig,l))*2)/(f_star(ig,l)+ &
1371	alim_star(ig,l))*2 + 2.rg(ztva(ig,l-1)-ztv(ig,l))/ztv(ig, l)( &
1372	zlev(ig,l+1)-zlev(ig,l))
1373	IF (w_est(ig,l+1)<0.) THEN
1374	w_est(ig, l+1) = zw2(ig, l)
1375	END IF
1376	IF (l>2) THEN
1377	IF ((w_est(ig,l+1)>w_est(ig,l)) .AND. (zlev(ig, &
1378	l+1)<zmax_sec(ig)) .AND. (zqla(ig,l-1)<1.E-10)) THEN
1379	detr_star(ig, l) = max(0., (rhobarz(ig, &
1380	l+1)sqrt(w_est(ig,l+1))sqrt(nu(ig,l)* &
1381	zlev(ig,l+1))-rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(nu(ig,l)* &
1382	zlev(ig,l)))/(r_aspect*zmax_sec(ig)))
1383	ELSE IF ((zlev(ig,l+1)<zmax_sec(ig)) .AND. (zqla(ig, &
1384	l-1)<1.E-10)) THEN
1385	detr_star(ig, l) = -f0(ig)*f_star(ig, lmix(ig))/(rhobarz(ig, &
1386	lmix(ig))wmaxa(ig))(rhobarz(ig,l+1)*sqrt(w_est(ig, &
1387	l+1))*((zmax_sec(ig)-zlev(ig,l+1))/((zmax_sec(ig)-zlev(ig, &
1388	lmix(ig)))))*2.-rhobarz(ig,l)sqrt(w_est(ig, &
1389	l))*((zmax_sec(ig)-zlev(ig,l))/((zmax_sec(ig)-zlev(ig,lmix(ig &
1390	)))))**2.)
1391	ELSE
1392	detr_star(ig, l) = 0.002f0(ig)f_star(ig, l)* &
1393	(zlev(ig,l+1)-zlev(ig,l))
1394
1395	END IF
1396	ELSE
1397	detr_star(ig, l) = 0.
1398	END IF
1399
1400	detr_star(ig, l) = detr_star(ig, l)/f0(ig)
1401	IF (nuage) THEN
1402	entr_star(ig, l) = 0.4*detr_star(ig, l)
1403	ELSE
1404	entr_star(ig, l) = 0.4*detr_star(ig, l)
1405	END IF
1406
1407	IF ((detr_star(ig,l))>f_star(ig,l)) THEN
1408	detr_star(ig, l) = f_star(ig, l)
1409	! entr_star(ig,l)=0.
1410	END IF
1411
1412	IF ((l<lentr(ig))) THEN
1413	entr_star(ig, l) = 0.
1414	! detr_star(ig,l)=0.
1415	END IF
1416
1417	! print*,'ok detr_star'
1418	END IF
1419	! prise en compte du detrainement dans le calcul du flux
1420	f_star(ig, l+1) = f_star(ig, l) + alim_star(ig, l) + &
1421	entr_star(ig, l) - detr_star(ig, l)
1422	! test
1423	! if (f_star(ig,l+1).lt.0.) then
1424	! f_star(ig,l+1)=0.
1425	! entr_star(ig,l)=0.
1426	! detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)
1427	! endif
1428	! test sur le signe de f_star
1429	IF (f_star(ig,l+1)>1.E-10) THEN
1430	! then
1431	! test
1432	! if (((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10)) then
1433	! AM on melange Tl et qt du thermique
1434	! on rajoute un excès de T dans la couche alim
1435	! if (l.lt.lentr(ig)) then
1436	! ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1437	! s
1438	! (alim_star(ig,l)+entr_star(ig,l))*(zthl(ig,l)+dtheta(ig,l)))
1439	! s /(f_star(ig,l+1)+detr_star(ig,l))
1440	! else
1441	ztla(ig, l) = (f_star(ig,l)*ztla(ig,l-1)+(alim_star(ig, &
1442	l)+entr_star(ig,l))*zthl(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1443	! s /(f_star(ig,l+1))
1444	! endif
1445	! on rajoute un excès de q dans la couche alim
1446	! if (l.lt.lentr(ig)) then
1447	! zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1448	! s (alim_star(ig,l)+entr_star(ig,l))*(po(ig,l)+0.001))
1449	! s /(f_star(ig,l+1)+detr_star(ig,l))
1450	! else
1451	zqta(ig, l) = (f_star(ig,l)*zqta(ig,l-1)+(alim_star(ig, &
1452	l)+entr_star(ig,l))*po(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1453	! s /(f_star(ig,l+1))
1454	! endif
1455	! AM on en deduit thetav et ql du thermique
1456	! CR test
1457	! Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1458	tbef(ig) = ztla(ig, l)*zpspsk(ig, l)
1459	zdelta = max(0., sign(1.,rtt-tbef(ig)))
1460	qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
1461	qsatbef(ig) = min(0.5, qsatbef(ig))
1462	zcor = 1./(1.-retv*qsatbef(ig))
1463	qsatbef(ig) = qsatbef(ig)*zcor
1464	zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig))>1.E-10)
1465
1466	IF (zsat(ig) .AND. (1==1)) THEN
1467	qlbef = max(0., zqta(ig,l)-qsatbef(ig))
1468	dt = 0.5rlvcpqlbef
1469	! write(17,*)'DT0=',DT
1470	DO WHILE (abs(dt)>ddt0)
1471	! print*,'aie'
1472	tbef(ig) = tbef(ig) + dt
1473	zdelta = max(0., sign(1.,rtt-tbef(ig)))
1474	qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
1475	qsatbef(ig) = min(0.5, qsatbef(ig))
1476	zcor = 1./(1.-retv*qsatbef(ig))
1477	qsatbef(ig) = qsatbef(ig)*zcor
1478	qlbef = zqta(ig, l) - qsatbef(ig)
1479
1480	zdelta = max(0., sign(1.,rtt-tbef(ig)))
1481	zcvm5 = r5les(1.-zdelta) + r5ieszdelta
1482	zcor = 1./(1.-retv*qsatbef(ig))
1483	dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
1484	num = -tbef(ig) + ztla(ig, l)zpspsk(ig, l) + rlvcpqlbef
1485	denom = 1. + rlvcp*dqsat_dt
1486	IF (denom<1.E-10) THEN
1487	PRINT *, 'pb denom'
1488	END IF
1489	dt = num/denom
1490	! write(17,*)'DT=',DT
1491	END DO
1492	zqla(ig, l) = max(0., zqta(ig,l)-qsatbef(ig))
1493	zqla(ig, l) = max(0., qlbef)
1494	! zqla(ig,l)=0.
1495	END IF
1496	! zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1497
1498	! on ecrit de maniere conservative (sat ou non)
1499	! T = Tl +Lv/Cp ql
1500	! CR rq utilisation de humidite specifique ou rapport de melange?
1501	ztva(ig, l) = ztla(ig, l)zpspsk(ig, l) + rlvcpzqla(ig, l)
1502	ztva(ig, l) = ztva(ig, l)/zpspsk(ig, l)
1503	! on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1504	zha(ig, l) = ztva(ig, l)
1505	! if (l.lt.lentr(ig)) then
1506	! ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)
1507	! s -zqla(ig,l))-zqla(ig,l)) + 0.1
1508	! else
1509	ztva(ig, l) = ztva(ig, l)(1.+retv(zqta(ig,l)-zqla(ig, &
1510	l))-zqla(ig,l))
1511	! endif
1512	! ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
1513	! s /(1.-retv*zqla(ig,l))
1514	! ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1515	! ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l)
1516	! s /(1.-retv*zqta(ig,l))
1517	! s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1518	! s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1519	! write(13,)zqla(ig,l),zqla(ig,l)/(1.-retvzqla(ig,l))
1520	! on ecrit zqsat
1521	zqsatth(ig, l) = qsatbef(ig)
1522	! enddo
1523	! DO ig=1,ngrid
1524	! if (zw2(ig,l).ge.1.e-10.and.
1525	! s f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
1526	! mise a jour de la vitesse ascendante (l'air entraine de la couche
1527	! consideree commence avec une vitesse nulle).
1528
1529	! if (f_star(ig,l+1).gt.1.e-10) then
1530	zw2(ig, l+1) = zw2(ig, l)* & ! s
1531	! ((f_star(ig,l)-detr_star(ig,l))**2)
1532	! s /f_star(ig,l+1)**2+
1533	((f_star(ig,l))2)/(f_star(ig,l+1)+detr_star(ig,l))2 + & ! s
1534	! /(f_star(ig,l+1))**2+
1535	2.rg(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
1536	! s (f_star(ig,l)/f_star(ig,l+1))*2
1537
1538	END IF
1539	END IF
1540
1541	IF (zw2(ig,l+1)<0.) THEN
1542	linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
1543	ig,l))
1544	zw2(ig, l+1) = 0.
1545	! print*,'linter=',linter(ig)
1546	! else if ((zw2(ig,l+1).lt.1.e-10).and.(zw2(ig,l+1).ge.0.)) then
1547	! linter(ig)=l+1
1548	! print*,'linter=l',zw2(ig,l),zw2(ig,l+1)
1549	ELSE
1550	wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
1551	! wa_moy(ig,l+1)=zw2(ig,l+1)
1552	END IF
1553	IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
1554	! lmix est le niveau de la couche ou w (wa_moy) est maximum
1555	lmix(ig) = l + 1
1556	wmaxa(ig) = wa_moy(ig, l+1)
1557	END IF
1558	END DO
1559	END DO
1560	PRINT *, 'fin calcul zw2'
1561
1562	! Calcul de la couche correspondant a la hauteur du thermique
1563	DO ig = 1, ngrid
1564	lmax(ig) = lentr(ig)
1565	END DO
1566	DO ig = 1, ngrid
1567	DO l = nlay, lentr(ig) + 1, -1
1568	IF (zw2(ig,l)<=1.E-10) THEN
1569	lmax(ig) = l - 1
1570	END IF
1571	END DO
1572	END DO
1573	! pas de thermique si couche 1 stable
1574	DO ig = 1, ngrid
1575	IF (lmin(ig)>1) THEN
1576	lmax(ig) = 1
1577	lmin(ig) = 1
1578	lentr(ig) = 1
1579	END IF
1580	END DO
1581
1582	! Determination de zw2 max
1583	DO ig = 1, ngrid
1584	wmax(ig) = 0.
1585	END DO
1586
1587	DO l = 1, nlay
1588	DO ig = 1, ngrid
1589	IF (l<=lmax(ig)) THEN
1590	IF (zw2(ig,l)<0.) THEN
1591	PRINT *, 'pb2 zw2<0'
1592	END IF
1593	zw2(ig, l) = sqrt(zw2(ig,l))
1594	wmax(ig) = max(wmax(ig), zw2(ig,l))
1595	ELSE
1596	zw2(ig, l) = 0.
1597	END IF
1598	END DO
1599	END DO
1600
1601	! Longueur caracteristique correspondant a la hauteur des thermiques.
1602	DO ig = 1, ngrid
1603	zmax(ig) = 0.
1604	zlevinter(ig) = zlev(ig, 1)
1605	END DO
1606	DO ig = 1, ngrid
1607	! calcul de zlevinter
1608	zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
1609	zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
1610	! pour le cas ou on prend tjs lmin=1
1611	! zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1612	zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,1))
1613	zmax0(ig) = zmax(ig)
1614	WRITE (11, *) 'ig,lmax,linter', ig, lmax(ig), linter(ig)
1615	WRITE (12, *) 'ig,zlevinter,zmax', ig, zmax(ig), zlevinter(ig)
1616	END DO
1617
1618	! Calcul de zmax_sec et wmax_sec
1619	CALL fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, f0, &
1620	zpspsk, alim, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, zmax_sec2, &
1621	wmax_sec2)
1622
1623	PRINT *, 'avant fermeture'
1624	! Fermeture,determination de f
1625	! en lmax f=d-e
1626	DO ig = 1, ngrid
1627	! entr_star(ig,lmax(ig))=0.
1628	! f_star(ig,lmax(ig)+1)=0.
1629	! detr_star(ig,lmax(ig))=f_star(ig,lmax(ig))+entr_star(ig,lmax(ig))
1630	! s +alim_star(ig,lmax(ig))
1631	END DO
1632
1633	DO ig = 1, ngrid
1634	alim_star2(ig) = 0.
1635	END DO
1636	! calcul de entr_star_tot
1637	DO ig = 1, ngrid
1638	DO k = 1, lmix(ig)
1639	entr_star_tot(ig) = entr_star_tot(ig) & ! s
1640	! +entr_star(ig,k)
1641	+alim_star(ig, k)
1642	! s -detr_star(ig,k)
1643	detr_star_tot(ig) = detr_star_tot(ig) & ! s
1644	! +alim_star(ig,k)
1645	-detr_star(ig, k) + entr_star(ig, k)
1646	END DO
1647	END DO
1648
1649	DO ig = 1, ngrid
1650	IF (alim_star_tot(ig)<1.E-10) THEN
1651	f(ig) = 0.
1652	ELSE
1653	! do k=lmin(ig),lentr(ig)
1654	DO k = 1, lentr(ig)
1655	alim_star2(ig) = alim_star2(ig) + alim_star(ig, k)*2/(rho(ig,k)( &
1656	zlev(ig,k+1)-zlev(ig,k)))
1657	END DO
1658	IF ((zmax_sec(ig)>1.E-10) .AND. (1==1)) THEN
1659	f(ig) = wmax_sec(ig)/(max(500.,zmax_sec(ig))r_aspectalim_star2(ig))
1660	f(ig) = f(ig) + (f0(ig)-f(ig))exp((-ptimestep/zmax_sec(ig))wmax_sec &
1661	(ig))
1662	ELSE
1663	f(ig) = wmax(ig)/(max(500.,zmax(ig))r_aspectalim_star2(ig))
1664	f(ig) = f(ig) + (f0(ig)-f(ig))exp((-ptimestep/zmax(ig))wmax(ig))
1665	END IF
1666	END IF
1667	f0(ig) = f(ig)
1668	END DO
1669	PRINT *, 'apres fermeture'
1670	! Calcul de l'entrainement
1671	DO ig = 1, ngrid
1672	DO k = 1, klev
1673	alim(ig, k) = f(ig)*alim_star(ig, k)
1674	END DO
1675	END DO
1676	! CR:test pour entrainer moins que la masse
1677	! do ig=1,ngrid
1678	! do l=1,lentr(ig)
1679	! if ((alim(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
1680	! alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
1681	! s -0.9*masse(ig,l)/ptimestep
1682	! alim(ig,l)=0.9*masse(ig,l)/ptimestep
1683	! endif
1684	! enddo
1685	! enddo
1686	! calcul du détrainement
1687	DO ig = 1, klon
1688	DO k = 1, klev
1689	detr(ig, k) = f(ig)*detr_star(ig, k)
1690	IF (detr(ig,k)<0.) THEN
1691	! print*,'detr1<0!!!'
1692	END IF
1693	END DO
1694	DO k = 1, klev
1695	entr(ig, k) = f(ig)*entr_star(ig, k)
1696	IF (entr(ig,k)<0.) THEN
1697	! print*,'entr1<0!!!'
1698	END IF
1699	END DO
1700	END DO
1701
1702	! do ig=1,ngrid
1703	! do l=1,klev
1704	! if (((detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep).gt.
1705	! s (masse(ig,l))) then
1706	! print*,'d2+e2+a2>m2','ig=',ig,'l=',l,'lmax(ig)=',lmax(ig),'d+e+a='
1707	! s,(detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep,'m=',masse(ig,l)
1708	! endif
1709	! enddo
1710	! enddo
1711	! Calcul des flux
1712
1713	DO ig = 1, ngrid
1714	DO l = 1, lmax(ig)
1715	! do l=1,klev
1716	! fmc(ig,l+1)=f(ig)*f_star(ig,l+1)
1717	fmc(ig, l+1) = fmc(ig, l) + alim(ig, l) + entr(ig, l) - detr(ig, l)
1718	! print*,'??!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1719	! s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1720	! s 'f+1=',fmc(ig,l+1)
1721	IF (fmc(ig,l+1)<0.) THEN
1722	PRINT *, 'fmc1<0', l + 1, lmax(ig), fmc(ig, l+1)
1723	fmc(ig, l+1) = fmc(ig, l)
1724	detr(ig, l) = alim(ig, l) + entr(ig, l)
1725	! fmc(ig,l+1)=0.
1726	! print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1727	END IF
1728	! if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1729	! f_old=fmc(ig,l+1)
1730	! fmc(ig,l+1)=fmc(ig,l)
1731	! detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1732	! endif
1733
1734	! if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1735	! f_old=fmc(ig,l+1)
1736	! fmc(ig,l+1)=fmc(ig,l)
1737	! detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l)
1738	! endif
1739	! rajout du test sur alpha croissant
1740	! if test
1741	! if (1.eq.0) then
1742
1743	IF (l==klev) THEN
1744	PRINT *, 'THERMCELL PB ig=', ig, ' l=', l
1745	abort_message = 'THERMCELL PB'
1746	CALL abort_physic(modname, abort_message, 1)
1747	END IF
1748	! if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1749	! s (l.ge.lentr(ig)).and.
1750	IF ((zw2(ig,l+1)>1.E-10) .AND. (zw2(ig,l)>1.E-10) .AND. (l>=lentr(ig))) &
1751	THEN
1752	IF (((fmc(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1)))>(fmc(ig,l)/ &
1753	(rhobarz(ig,l)*zw2(ig,l))))) THEN
1754	f_old = fmc(ig, l+1)
1755	fmc(ig, l+1) = fmc(ig, l)rhobarz(ig, l+1)zw2(ig, l+1)/ &
1756	(rhobarz(ig,l)*zw2(ig,l))
1757	detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1758	! detr(ig,l)=(fmc(ig,l+1)-fmc(ig,l))/(0.4-1.)
1759	! entr(ig,l)=0.4*detr(ig,l)
1760	! entr(ig,l)=fmc(ig,l+1)-fmc(ig,l)+detr(ig,l)
1761	END IF
1762	END IF
1763	IF ((fmc(ig,l+1)>fmc(ig,l)) .AND. (l>lentr(ig))) THEN
1764	f_old = fmc(ig, l+1)
1765	fmc(ig, l+1) = fmc(ig, l)
1766	detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1767	END IF
1768	IF (detr(ig,l)>fmc(ig,l)) THEN
1769	detr(ig, l) = fmc(ig, l)
1770	entr(ig, l) = fmc(ig, l+1) - alim(ig, l)
1771	END IF
1772	IF (fmc(ig,l+1)<0.) THEN
1773	detr(ig, l) = detr(ig, l) + fmc(ig, l+1)
1774	fmc(ig, l+1) = 0.
1775	PRINT *, 'fmc2<0', l + 1, lmax(ig)
1776	END IF
1777
1778	! test pour ne pas avoir f=0 et d=e/=0
1779	! if (fmc(ig,l+1).lt.1.e-10) then
1780	! detr(ig,l+1)=0.
1781	! entr(ig,l+1)=0.
1782	! zqla(ig,l+1)=0.
1783	! zw2(ig,l+1)=0.
1784	! lmax(ig)=l+1
1785	! zmax(ig)=zlev(ig,lmax(ig))
1786	! endif
1787	IF (zw2(ig,l+1)>1.E-10) THEN
1788	IF ((((fmc(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1)))>1.)) THEN
1789	f_old = fmc(ig, l+1)
1790	fmc(ig, l+1) = rhobarz(ig, l+1)*zw2(ig, l+1)
1791	zw2(ig, l+1) = 0.
1792	zqla(ig, l+1) = 0.
1793	detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1794	lmax(ig) = l + 1
1795	zmax(ig) = zlev(ig, lmax(ig))
1796	PRINT *, 'alpha>1', l + 1, lmax(ig)
1797	END IF
1798	END IF
1799	! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
1800	! endif test
1801	! endif
1802	END DO
1803	END DO
1804	DO ig = 1, ngrid
1805	! if (fmc(ig,lmax(ig)+1).ne.0.) then
1806	fmc(ig, lmax(ig)+1) = 0.
1807	entr(ig, lmax(ig)) = 0.
1808	detr(ig, lmax(ig)) = fmc(ig, lmax(ig)) + entr(ig, lmax(ig)) + &
1809	alim(ig, lmax(ig))
1810	! endif
1811	END DO
1812	! test sur le signe de fmc
1813	DO ig = 1, ngrid
1814	DO l = 1, klev + 1
1815	IF (fmc(ig,l)<0.) THEN
1816	PRINT *, 'fm1<0!!!', 'ig=', ig, 'l=', l, 'a=', alim(ig, l-1), 'e=', &
1817	entr(ig, l-1), 'f=', fmc(ig, l-1), 'd=', detr(ig, l-1), 'f+1=', &
1818	fmc(ig, l)
1819	END IF
1820	END DO
1821	END DO
1822	! test de verification
1823	DO ig = 1, ngrid
1824	DO l = 1, lmax(ig)
1825	IF ((abs(fmc(ig,l+1)-fmc(ig,l)-alim(ig,l)-entr(ig,l)+ &
1826	detr(ig,l)))>1.E-4) THEN
1827	! print*,'pbcm!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1828	! s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1829	! s 'f+1=',fmc(ig,l+1)
1830	END IF
1831	IF (detr(ig,l)<0.) THEN
1832	PRINT *, 'detrdemi<0!!!'
1833	END IF
1834	END DO
1835	END DO
1836
1837	! RC
1838	! CR def de zmix continu (profil parabolique des vitesses)
1839	DO ig = 1, ngrid
1840	IF (lmix(ig)>1.) THEN
1841	! test
1842	IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
1843	(zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
1844	zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
1845	(zlev(ig,lmix(ig)))))>1E-10) THEN
1846
1847	zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
1848	)2-(zlev(ig,lmix(ig)+1))2)-(zw2(ig,lmix(ig))-zw2(ig, &
1849	lmix(ig)+1))((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))/ &
1850	(2.((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))((zlev(ig,lmix(ig)))- &
1851	(zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
1852	zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
1853	ELSE
1854	zmix(ig) = zlev(ig, lmix(ig))
1855	PRINT *, 'pb zmix'
1856	END IF
1857	ELSE
1858	zmix(ig) = 0.
1859	END IF
1860	! test
1861	IF ((zmax(ig)-zmix(ig))<=0.) THEN
1862	zmix(ig) = 0.9*zmax(ig)
1863	! print*,'pb zmix>zmax'
1864	END IF
1865	END DO
1866	DO ig = 1, klon
1867	zmix0(ig) = zmix(ig)
1868	END DO
1869
1870	! calcul du nouveau lmix correspondant
1871	DO ig = 1, ngrid
1872	DO l = 1, klev
1873	IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
1874	lmix(ig) = l
1875	END IF
1876	END DO
1877	END DO
1878
1879	! ne devrait pas arriver!!!!!
1880	DO ig = 1, ngrid
1881	DO l = 1, klev
1882	IF (detr(ig,l)>(fmc(ig,l)+alim(ig,l))+entr(ig,l)) THEN
1883	PRINT *, 'detr2>fmc2!!!', 'ig=', ig, 'l=', l, 'd=', detr(ig, l), &
1884	'f=', fmc(ig, l), 'lmax=', lmax(ig)
1885	! detr(ig,l)=fmc(ig,l)+alim(ig,l)+entr(ig,l)
1886	! entr(ig,l)=0.
1887	! fmc(ig,l+1)=0.
1888	! zw2(ig,l+1)=0.
1889	! zqla(ig,l+1)=0.
1890	PRINT *, 'pb!fm=0 et f_star>0', l, lmax(ig)
1891	! lmax(ig)=l
1892	END IF
1893	END DO
1894	END DO
1895	DO ig = 1, ngrid
1896	DO l = lmax(ig) + 1, klev + 1
1897	! fmc(ig,l)=0.
1898	! detr(ig,l)=0.
1899	! entr(ig,l)=0.
1900	! zw2(ig,l)=0.
1901	! zqla(ig,l)=0.
1902	END DO
1903	END DO
1904
1905	! Calcul du detrainement lors du premier passage
1906	! print*,'9 OK convect8'
1907	! print*,'WA1 ',wa_moy
1908
1909	! determination de l'indice du debut de la mixed layer ou w decroit
1910
1911	! calcul de la largeur de chaque ascendance dans le cas conservatif.
1912	! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
1913	! d'une couche est égale à la hauteur de la couche alimentante.
1914	! La vitesse maximale dans l'ascendance est aussi prise comme estimation
1915	! de la vitesse d'entrainement horizontal dans la couche alimentante.
1916
1917	DO l = 2, nlay
1918	DO ig = 1, ngrid
1919	IF (l<=lmax(ig) .AND. (test(ig)==1)) THEN
1920	zw = max(wa_moy(ig,l), 1.E-10)
1921	larg_cons(ig, l) = zmax(ig)r_aspectfmc(ig, l)/(rhobarz(ig,l)*zw)
1922	END IF
1923	END DO
1924	END DO
1925
1926	DO l = 2, nlay
1927	DO ig = 1, ngrid
1928	IF (l<=lmax(ig) .AND. (test(ig)==1)) THEN
1929	! if (idetr.eq.0) then
1930	! cette option est finalement en dur.
1931	IF ((l_mix*zlev(ig,l))<0.) THEN
1932	PRINT , 'pb l_mixzlev<0'
1933	END IF
1934	! CR: test: nouvelle def de lambda
1935	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1936	IF (zw2(ig,l)>1.E-10) THEN
1937	larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
1938	ELSE
1939	larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
1940	END IF
1941	! else if (idetr.eq.1) then
1942	! larg_detr(ig,l)=larg_cons(ig,l)
1943	! s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
1944	! else if (idetr.eq.2) then
1945	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1946	! s *sqrt(wa_moy(ig,l))
1947	! else if (idetr.eq.4) then
1948	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1949	! s *wa_moy(ig,l)
1950	! endif
1951	END IF
1952	END DO
1953	END DO
1954
1955	! print*,'10 OK convect8'
1956	! print*,'WA2 ',wa_moy
1957	! cal1cul de la fraction de la maille concernée par l'ascendance en tenant
1958	! compte de l'epluchage du thermique.
1959
1960
1961	DO l = 2, nlay
1962	DO ig = 1, ngrid
1963	IF (larg_cons(ig,l)>1. .AND. (test(ig)==1)) THEN
1964	! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
1965	fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
1966	! test
1967	fraca(ig, l) = max(fraca(ig,l), 0.)
1968	fraca(ig, l) = min(fraca(ig,l), 0.5)
1969	fracd(ig, l) = 1. - fraca(ig, l)
1970	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
1971	ELSE
1972	! wa_moy(ig,l)=0.
1973	fraca(ig, l) = 0.
1974	fracc(ig, l) = 0.
1975	fracd(ig, l) = 1.
1976	END IF
1977	END DO
1978	END DO
1979	! CR: calcul de fracazmix
1980	DO ig = 1, ngrid
1981	IF (test(ig)==1) THEN
1982	fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
1983	(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
1984	fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca( &
1985	ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
1986	END IF
1987	END DO
1988
1989	DO l = 2, nlay
1990	DO ig = 1, ngrid
1991	IF (larg_cons(ig,l)>1. .AND. (test(ig)==1)) THEN
1992	IF (l>lmix(ig)) THEN
1993	! test
1994	IF (zmax(ig)-zmix(ig)<1.E-10) THEN
1995	! print*,'pb xxx'
1996	xxx(ig, l) = (lmax(ig)+1.-l)/(lmax(ig)+1.-lmix(ig))
1997	ELSE
1998	xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
1999	END IF
2000	IF (idetr==0) THEN
2001	fraca(ig, l) = fracazmix(ig)
2002	ELSE IF (idetr==1) THEN
2003	fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
2004	ELSE IF (idetr==2) THEN
2005	fraca(ig, l) = fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
2006	ELSE
2007	fraca(ig, l) = fracazmix(ig)xxx(ig, l)*2
2008	END IF
2009	! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
2010	fraca(ig, l) = max(fraca(ig,l), 0.)
2011	fraca(ig, l) = min(fraca(ig,l), 0.5)
2012	fracd(ig, l) = 1. - fraca(ig, l)
2013	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
2014	END IF
2015	END IF
2016	END DO
2017	END DO
2018
2019	PRINT *, 'fin calcul fraca'
2020	! print*,'11 OK convect8'
2021	! print*,'Ea3 ',wa_moy
2022	! ------------------------------------------------------------------
2023	! Calcul de fracd, wd
2024	! somme wa - wd = 0
2025	! ------------------------------------------------------------------
2026
2027
2028	DO ig = 1, ngrid
2029	fm(ig, 1) = 0.
2030	fm(ig, nlay+1) = 0.
2031	END DO
2032
2033	DO l = 2, nlay
2034	DO ig = 1, ngrid
2035	IF (test(ig)==1) THEN
2036	fm(ig, l) = fraca(ig, l)wa_moy(ig, l)rhobarz(ig, l)
2037	! CR:test
2038	IF (alim(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) &
2039	THEN
2040	fm(ig, l) = fm(ig, l-1)
2041	! write(1,*)'ajustement fm, l',l
2042	END IF
2043	! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
2044	! RC
2045	END IF
2046	END DO
2047	DO ig = 1, ngrid
2048	IF (fracd(ig,l)<0.1 .AND. (test(ig)==1)) THEN
2049	abort_message = 'fracd trop petit'
2050	CALL abort_physic(modname, abort_message, 1)
2051	ELSE
2052	! vitesse descendante "diagnostique"
2053	wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
2054	END IF
2055	END DO
2056	END DO
2057
2058	DO l = 1, nlay + 1
2059	DO ig = 1, ngrid
2060	IF (test(ig)==0) THEN
2061	fm(ig, l) = fmc(ig, l)
2062	END IF
2063	END DO
2064	END DO
2065
2066	! fin du first
2067	DO l = 1, nlay
2068	DO ig = 1, ngrid
2069	! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
2070	masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
2071	END DO
2072	END DO
2073
2074	! print*,'12 OK convect8'
2075	! print*,'WA4 ',wa_moy
2076	! c------------------------------------------------------------------
2077	! calcul du transport vertical
2078	! ------------------------------------------------------------------
2079
2080	GO TO 4444
2081	! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
2082	DO l = 2, nlay - 1
2083	DO ig = 1, ngrid
2084	IF (fm(ig,l+1)ptimestep>masse(ig,l) .AND. fm(ig,l+1)ptimestep>masse( &
2085	ig,l+1)) THEN
2086	PRINT *, 'WARN!!! FM>M ig=', ig, ' l=', l, ' FM=', &
2087	fm(ig, l+1)*ptimestep, ' M=', masse(ig, l), masse(ig, l+1)
2088	END IF
2089	END DO
2090	END DO
2091
2092	DO l = 1, nlay
2093	DO ig = 1, ngrid
2094	IF ((alim(ig,l)+entr(ig,l))*ptimestep>masse(ig,l)) THEN
2095	PRINT *, 'WARN!!! E>M ig=', ig, ' l=', l, ' E==', &
2096	(entr(ig,l)+alim(ig,l))*ptimestep, ' M=', masse(ig, l)
2097	END IF
2098	END DO
2099	END DO
2100
2101	DO l = 1, nlay
2102	DO ig = 1, ngrid
2103	IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
2104	! print*,'WARN!!! fm exagere ig=',ig,' l=',l
2105	! s ,' FM=',fm(ig,l)
2106	END IF
2107	IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
2108	! print*,'WARN!!! masse exagere ig=',ig,' l=',l
2109	! s ,' M=',masse(ig,l)
2110	! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
2111	! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
2112	! print*,'zlev(ig,l+1),zlev(ig,l)'
2113	! s ,zlev(ig,l+1),zlev(ig,l)
2114	! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
2115	! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
2116	END IF
2117	IF (.NOT. alim(ig,l)>=0. .OR. .NOT. alim(ig,l)<=10.) THEN
2118	! print*,'WARN!!! entr exagere ig=',ig,' l=',l
2119	! s ,' E=',entr(ig,l)
2120	END IF
2121	END DO
2122	END DO
2123
2124	4444 CONTINUE
2125
2126	! CR:redefinition du entr
2127	! CR:test:on ne change pas la def du entr mais la def du fm
2128	DO l = 1, nlay
2129	DO ig = 1, ngrid
2130	IF (test(ig)==1) THEN
2131	detr(ig, l) = fm(ig, l) + alim(ig, l) - fm(ig, l+1)
2132	IF (detr(ig,l)<0.) THEN
2133	! entr(ig,l)=entr(ig,l)-detr(ig,l)
2134	fm(ig, l+1) = fm(ig, l) + alim(ig, l)
2135	detr(ig, l) = 0.
2136	! write(11,*)'l,ig,entr',l,ig,entr(ig,l)
2137	! print*,'WARNING !!! detrainement negatif ',ig,l
2138	END IF
2139	END IF
2140	END DO
2141	END DO
2142	! RC
2143
2144	IF (w2di==1) THEN
2145	fm0 = fm0 + ptimestep*(fm-fm0)/tho
2146	entr0 = entr0 + ptimestep*(alim+entr-entr0)/tho
2147	ELSE
2148	fm0 = fm
2149	entr0 = alim + entr
2150	detr0 = detr
2151	alim0 = alim
2152	! zoa=zqta
2153	! entr0=alim
2154	END IF
2155
2156	IF (1==1) THEN
2157	! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2158	! . ,zh,zdhadj,zha)
2159	! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2160	! . ,zo,pdoadj,zoa)
2161	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zthl, &
2162	zdthladj, zta)
2163	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, po, pdoadj, &
2164	zoa)
2165	ELSE
2166	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
2167	zdhadj, zha)
2168	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
2169	pdoadj, zoa)
2170	END IF
2171
2172	IF (1==0) THEN
2173	CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
2174	zu, zv, pduadj, pdvadj, zua, zva)
2175	ELSE
2176	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
2177	zua)
2178	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
2179	zva)
2180	END IF
2181
2182	! Calcul des moments
2183	! do l=1,nlay
2184	! do ig=1,ngrid
2185	! zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
2186	! zf2=zf/(1.-zf)
2187	! thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
2188	! wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
2189	! enddo
2190	! enddo
2191
2192
2193
2194
2195
2196
2197	! print*,'13 OK convect8'
2198	! print*,'WA5 ',wa_moy
2199	DO l = 1, nlay
2200	DO ig = 1, ngrid
2201	! pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
2202	pdtadj(ig, l) = zdthladj(ig, l)*zpspsk(ig, l)
2203	END DO
2204	END DO
2205
2206
2207	! do l=1,nlay
2208	! do ig=1,ngrid
2209	! if(abs(pdtadj(ig,l))*86400..gt.500.) then
2210	! print*,'WARN!!! ig=',ig,' l=',l
2211	! s ,' pdtadj=',pdtadj(ig,l)
2212	! endif
2213	! if(abs(pdoadj(ig,l))*86400..gt.1.) then
2214	! print*,'WARN!!! ig=',ig,' l=',l
2215	! s ,' pdoadj=',pdoadj(ig,l)
2216	! endif
2217	! enddo
2218	! enddo
2219
2220	! print*,'14 OK convect8'
2221	! ------------------------------------------------------------------
2222	! Calculs pour les sorties
2223	! ------------------------------------------------------------------
2224	! calcul de fraca pour les sorties
2225	DO l = 2, klev
2226	DO ig = 1, klon
2227	IF (zw2(ig,l)>1.E-10) THEN
2228	fraca(ig, l) = fm(ig, l)/(rhobarz(ig,l)*zw2(ig,l))
2229	ELSE
2230	fraca(ig, l) = 0.
2231	END IF
2232	END DO
2233	END DO
2234	IF (sorties) THEN
2235	DO l = 1, nlay
2236	DO ig = 1, ngrid
2237	zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
2238	zld(ig, l) = fracd(ig, l)*zmax(ig)
2239	IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
2240	(1.-fracd(ig,l))
2241	END DO
2242	END DO
2243	! CR calcul du niveau de condensation
2244	! initialisation
2245	DO ig = 1, ngrid
2246	nivcon(ig) = 0.
2247	zcon(ig) = 0.
2248	END DO
2249	DO k = nlay, 1, -1
2250	DO ig = 1, ngrid
2251	IF (zqla(ig,k)>1E-10) THEN
2252	nivcon(ig) = k
2253	zcon(ig) = zlev(ig, k)
2254	END IF
2255	! if (zcon(ig).gt.1.e-10) then
2256	! nuage=.true.
2257	! else
2258	! nuage=.false.
2259	! endif
2260	END DO
2261	END DO
2262
2263	DO l = 1, nlay
2264	DO ig = 1, ngrid
2265	zf = fraca(ig, l)
2266	zf2 = zf/(1.-zf)
2267	thetath2(ig, l) = zf2(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))*2
2268	wth2(ig, l) = zf2(zw2(ig,l))*2
2269	! print*,'wth2=',wth2(ig,l)
2270	wth3(ig, l) = zf2(1-2.fraca(ig,l))/(1-fraca(ig,l))zw2(ig, l) &
2271	zw2(ig, l)*zw2(ig, l)
2272	q2(ig, l) = zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
2273	! test: on calcul q2/po=ratqsc
2274	! if (nuage) then
2275	ratqscth(ig, l) = sqrt(q2(ig,l))/(po(ig,l)*1000.)
2276	! else
2277	! ratqscth(ig,l)=0.
2278	! endif
2279	END DO
2280	END DO
2281	! calcul du ratqscdiff
2282	sum = 0.
2283	sumdiff = 0.
2284	ratqsdiff(:, :) = 0.
2285	DO ig = 1, ngrid
2286	DO l = 1, lentr(ig)
2287	sum = sum + alim_star(ig, l)zqta(ig, l)1000.
2288	END DO
2289	END DO
2290	DO ig = 1, ngrid
2291	DO l = 1, lentr(ig)
2292	zf = fraca(ig, l)
2293	zf2 = zf/(1.-zf)
2294	sumdiff = sumdiff + alim_star(ig, l)(zqta(ig,l)1000.-sum)**2
2295	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
2296	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
2297	END DO
2298	END DO
2299	DO l = 1, klev
2300	DO ig = 1, ngrid
2301	ratqsdiff(ig, l) = sqrt(sumdiff)/(po(ig,l)*1000.)
2302	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
2303	END DO
2304	END DO
2305
2306	END IF
2307
2308	! print*,'19 OK convect8'
2309	RETURN
2310	END SUBROUTINE thermcell_cld
2311
2312	SUBROUTINE thermcell_eau(ngrid, nlay, ptimestep, pplay, pplev, pphi, pu, pv, &
2313	pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
2314	! ,pu_therm,pv_therm
2315	, r_aspect, l_mix, w2di, tho)
2316
2317	USE dimphy
2318	IMPLICIT NONE
2319
2320	! =======================================================================
2321
2322	! Calcul du transport verticale dans la couche limite en presence
2323	! de "thermiques" explicitement representes
2324
2325	! Réécriture à partir d'un listing papier à Habas, le 14/02/00
2326
2327	! le thermique est supposé homogène et dissipé par mélange avec
2328	! son environnement. la longueur l_mix contrôle l'efficacité du
2329	! mélange
2330
2331	! Le calcul du transport des différentes espèces se fait en prenant
2332	! en compte:
2333	! 1. un flux de masse montant
2334	! 2. un flux de masse descendant
2335	! 3. un entrainement
2336	! 4. un detrainement
2337
2338	! =======================================================================
2339
2340	! -----------------------------------------------------------------------
2341	! declarations:
2342	! -------------
2343
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	include "YOMCST.h"
3276
3277	! arguments:
3278	! ----------
3279
3280	INTEGER ngrid, nlay, w2di
3281	REAL tho
3282	REAL ptimestep, l_mix, r_aspect
3283	REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
3284	REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
3285	REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
3286	REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
3287	REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
3288	REAL pphi(ngrid, nlay)
3289
3290	INTEGER idetr
3291	SAVE idetr
3292	DATA idetr/3/
3293	!$OMP THREADPRIVATE(idetr)
3294
3295	! local:
3296	! ------
3297
3298	INTEGER ig, k, l, lmaxa(klon), lmix(klon)
3299	REAL zsortie1d(klon)
3300	! CR: on remplace lmax(klon,klev+1)
3301	INTEGER lmax(klon), lmin(klon), lentr(klon)
3302	REAL linter(klon)
3303	REAL zmix(klon), fracazmix(klon)
3304	! RC
3305	REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
3306
3307	REAL zlev(klon, klev+1), zlay(klon, klev)
3308	REAL zh(klon, klev), zdhadj(klon, klev)
3309	REAL ztv(klon, klev)
3310	REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
3311	REAL wh(klon, klev+1)
3312	REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
3313	REAL zla(klon, klev+1)
3314	REAL zwa(klon, klev+1)
3315	REAL zld(klon, klev+1)
3316	REAL zwd(klon, klev+1)
3317	REAL zsortie(klon, klev)
3318	REAL zva(klon, klev)
3319	REAL zua(klon, klev)
3320	REAL zoa(klon, klev)
3321
3322	REAL zha(klon, klev)
3323	REAL wa_moy(klon, klev+1)
3324	REAL fraca(klon, klev+1)
3325	REAL fracc(klon, klev+1)
3326	REAL zf, zf2
3327	REAL thetath2(klon, klev), wth2(klon, klev)
3328	! common/comtherm/thetath2,wth2
3329
3330	REAL count_time
3331	INTEGER ialt
3332
3333	LOGICAL sorties
3334	REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
3335	REAL zpspsk(klon, klev)
3336
3337	! real wmax(klon,klev),wmaxa(klon)
3338	REAL wmax(klon), wmaxa(klon)
3339	REAL wa(klon, klev, klev+1)
3340	REAL wd(klon, klev+1)
3341	REAL larg_part(klon, klev, klev+1)
3342	REAL fracd(klon, klev+1)
3343	REAL xxx(klon, klev+1)
3344	REAL larg_cons(klon, klev+1)
3345	REAL larg_detr(klon, klev+1)
3346	REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
3347	REAL pu_therm(klon, klev), pv_therm(klon, klev)
3348	REAL fm(klon, klev+1), entr(klon, klev)
3349	REAL fmc(klon, klev+1)
3350
3351	! CR:nouvelles variables
3352	REAL f_star(klon, klev+1), entr_star(klon, klev)
3353	REAL entr_star_tot(klon), entr_star2(klon)
3354	REAL f(klon), f0(klon)
3355	REAL zlevinter(klon)
3356	LOGICAL first
3357	DATA first/.FALSE./
3358	SAVE first
3359	!$OMP THREADPRIVATE(first)
3360	! RC
3361
3362	CHARACTER *2 str2
3363	CHARACTER *10 str10
3364
3365	CHARACTER (LEN=20) :: modname = 'thermcell'
3366	CHARACTER (LEN=80) :: abort_message
3367
3368	LOGICAL vtest(klon), down
3369
3370	EXTERNAL scopy
3371
3372	INTEGER ncorrec, ll
3373	SAVE ncorrec
3374	DATA ncorrec/0/
3375	!$OMP THREADPRIVATE(ncorrec)
3376
3377
3378	! -----------------------------------------------------------------------
3379	! initialisation:
3380	! ---------------
3381
3382	sorties = .TRUE.
3383	IF (ngrid/=klon) THEN
3384	PRINT *
3385	PRINT *, 'STOP dans convadj'
3386	PRINT *, 'ngrid =', ngrid
3387	PRINT *, 'klon =', klon
3388	END IF
3389
3390	! -----------------------------------------------------------------------
3391	! incrementation eventuelle de tendances precedentes:
3392	! ---------------------------------------------------
3393
3394	! print*,'0 OK convect8'
3395
3396	DO l = 1, nlay
3397	DO ig = 1, ngrid
3398	zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
3399	zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
3400	zu(ig, l) = pu(ig, l)
3401	zv(ig, l) = pv(ig, l)
3402	zo(ig, l) = po(ig, l)
3403	ztv(ig, l) = zh(ig, l)(1.+0.61zo(ig,l))
3404	END DO
3405	END DO
3406
3407	! print*,'1 OK convect8'
3408	! --------------------
3409
3410
3411	! + + + + + + + + + + +
3412
3413
3414	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
3415	! wh,wt,wo ...
3416
3417	! + + + + + + + + + + + zh,zu,zv,zo,rho
3418
3419
3420	! -------------------- zlev(1)
3421	! \\\\\\\\\\\\\\\\\\\\
3422
3423
3424
3425	! -----------------------------------------------------------------------
3426	! Calcul des altitudes des couches
3427	! -----------------------------------------------------------------------
3428
3429	DO l = 2, nlay
3430	DO ig = 1, ngrid
3431	zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
3432	END DO
3433	END DO
3434	DO ig = 1, ngrid
3435	zlev(ig, 1) = 0.
3436	zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
3437	END DO
3438	DO l = 1, nlay
3439	DO ig = 1, ngrid
3440	zlay(ig, l) = pphi(ig, l)/rg
3441	END DO
3442	END DO
3443
3444	! print*,'2 OK convect8'
3445	! -----------------------------------------------------------------------
3446	! Calcul des densites
3447	! -----------------------------------------------------------------------
3448
3449	DO l = 1, nlay
3450	DO ig = 1, ngrid
3451	rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)rdzh(ig,l))
3452	END DO
3453	END DO
3454
3455	DO l = 2, nlay
3456	DO ig = 1, ngrid
3457	rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
3458	END DO
3459	END DO
3460
3461	DO k = 1, nlay
3462	DO l = 1, nlay + 1
3463	DO ig = 1, ngrid
3464	wa(ig, k, l) = 0.
3465	END DO
3466	END DO
3467	END DO
3468
3469	! print*,'3 OK convect8'
3470	! ------------------------------------------------------------------
3471	! Calcul de w2, quarre de w a partir de la cape
3472	! a partir de w2, on calcule wa, vitesse de l'ascendance
3473
3474	! ATTENTION: Dans cette version, pour cause d'economie de memoire,
3475	! w2 est stoke dans wa
3476
3477	! ATTENTION: dans convect8, on n'utilise le calcule des wa
3478	! independants par couches que pour calculer l'entrainement
3479	! a la base et la hauteur max de l'ascendance.
3480
3481	! Indicages:
3482	! l'ascendance provenant du niveau k traverse l'interface l avec
3483	! une vitesse wa(k,l).
3484
3485	! --------------------
3486
3487	! + + + + + + + + + +
3488
3489	! wa(k,l) ---- -------------------- l
3490	! /\
3491	! /\|\|\ + + + + + + + + + +
3492	! \|\|
3493	! \|\| --------------------
3494	! \|\|
3495	! \|\| + + + + + + + + + +
3496	! \|\|
3497	! \|\| --------------------
3498	! \|\|__
3499	! \|___ + + + + + + + + + + k
3500
3501	! --------------------
3502
3503
3504
3505	! ------------------------------------------------------------------
3506
3507	! CR: ponderation entrainement des couches instables
3508	! def des entr_star tels que entr=f*entr_star
3509	DO l = 1, klev
3510	DO ig = 1, ngrid
3511	entr_star(ig, l) = 0.
3512	END DO
3513	END DO
3514	! determination de la longueur de la couche d entrainement
3515	DO ig = 1, ngrid
3516	lentr(ig) = 1
3517	END DO
3518
3519	! on ne considere que les premieres couches instables
3520	DO k = nlay - 2, 1, -1
3521	DO ig = 1, ngrid
3522	IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
3523	lentr(ig) = k
3524	END IF
3525	END DO
3526	END DO
3527
3528	! determination du lmin: couche d ou provient le thermique
3529	DO ig = 1, ngrid
3530	lmin(ig) = 1
3531	END DO
3532	DO ig = 1, ngrid
3533	DO l = nlay, 2, -1
3534	IF (ztv(ig,l-1)>ztv(ig,l)) THEN
3535	lmin(ig) = l - 1
3536	END IF
3537	END DO
3538	END DO
3539
3540	! definition de l'entrainement des couches
3541	DO l = 1, klev - 1
3542	DO ig = 1, ngrid
3543	IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
3544	entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))*(zlev(ig,l+1)-zlev(ig,l))
3545	END IF
3546	END DO
3547	END DO
3548	! pas de thermique si couches 1->5 stables
3549	DO ig = 1, ngrid
3550	IF (lmin(ig)>5) THEN
3551	DO l = 1, klev
3552	entr_star(ig, l) = 0.
3553	END DO
3554	END IF
3555	END DO
3556	! calcul de l entrainement total
3557	DO ig = 1, ngrid
3558	entr_star_tot(ig) = 0.
3559	END DO
3560	DO ig = 1, ngrid
3561	DO k = 1, klev
3562	entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
3563	END DO
3564	END DO
3565
3566	PRINT *, 'fin calcul entr_star'
3567	DO k = 1, klev
3568	DO ig = 1, ngrid
3569	ztva(ig, k) = ztv(ig, k)
3570	END DO
3571	END DO
3572	! RC
3573	! print*,'7 OK convect8'
3574	DO k = 1, klev + 1
3575	DO ig = 1, ngrid
3576	zw2(ig, k) = 0.
3577	fmc(ig, k) = 0.
3578	! CR
3579	f_star(ig, k) = 0.
3580	! RC
3581	larg_cons(ig, k) = 0.
3582	larg_detr(ig, k) = 0.
3583	wa_moy(ig, k) = 0.
3584	END DO
3585	END DO
3586
3587	! print*,'8 OK convect8'
3588	DO ig = 1, ngrid
3589	linter(ig) = 1.
3590	lmaxa(ig) = 1
3591	lmix(ig) = 1
3592	wmaxa(ig) = 0.
3593	END DO
3594
3595	! CR:
3596	DO l = 1, nlay - 2
3597	DO ig = 1, ngrid
3598	IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
3599	zw2(ig,l)<1E-10) THEN
3600	f_star(ig, l+1) = entr_star(ig, l)
3601	! test:calcul de dteta
3602	zw2(ig, l+1) = 2.rg(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
3603	(zlev(ig,l+1)-zlev(ig,l))0.4pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
3604	larg_detr(ig, l) = 0.
3605	ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
3606	l)>1.E-10)) THEN
3607	f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
3608	ztva(ig, l) = (f_star(ig,l)ztva(ig,l-1)+entr_star(ig,l)ztv(ig,l))/ &
3609	f_star(ig, l+1)
3610	zw2(ig, l+1) = zw2(ig, l)(f_star(ig,l)/f_star(ig,l+1))*2 + &
3611	2.rg(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
3612	END IF
3613	! determination de zmax continu par interpolation lineaire
3614	IF (zw2(ig,l+1)<0.) THEN
3615	! test
3616	IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
3617	PRINT *, 'pb linter'
3618	END IF
3619	linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
3620	ig,l))
3621	zw2(ig, l+1) = 0.
3622	lmaxa(ig) = l
3623	ELSE
3624	IF (zw2(ig,l+1)<0.) THEN
3625	PRINT *, 'pb1 zw2<0'
3626	END IF
3627	wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
3628	END IF
3629	IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
3630	! lmix est le niveau de la couche ou w (wa_moy) est maximum
3631	lmix(ig) = l + 1
3632	wmaxa(ig) = wa_moy(ig, l+1)
3633	END IF
3634	END DO
3635	END DO
3636	PRINT *, 'fin calcul zw2'
3637
3638	! Calcul de la couche correspondant a la hauteur du thermique
3639	DO ig = 1, ngrid
3640	lmax(ig) = lentr(ig)
3641	END DO
3642	DO ig = 1, ngrid
3643	DO l = nlay, lentr(ig) + 1, -1
3644	IF (zw2(ig,l)<=1.E-10) THEN
3645	lmax(ig) = l - 1
3646	END IF
3647	END DO
3648	END DO
3649	! pas de thermique si couches 1->5 stables
3650	DO ig = 1, ngrid
3651	IF (lmin(ig)>5) THEN
3652	lmax(ig) = 1
3653	lmin(ig) = 1
3654	END IF
3655	END DO
3656
3657	! Determination de zw2 max
3658	DO ig = 1, ngrid
3659	wmax(ig) = 0.
3660	END DO
3661
3662	DO l = 1, nlay
3663	DO ig = 1, ngrid
3664	IF (l<=lmax(ig)) THEN
3665	IF (zw2(ig,l)<0.) THEN
3666	PRINT *, 'pb2 zw2<0'
3667	END IF
3668	zw2(ig, l) = sqrt(zw2(ig,l))
3669	wmax(ig) = max(wmax(ig), zw2(ig,l))
3670	ELSE
3671	zw2(ig, l) = 0.
3672	END IF
3673	END DO
3674	END DO
3675
3676	! Longueur caracteristique correspondant a la hauteur des thermiques.
3677	DO ig = 1, ngrid
3678	zmax(ig) = 0.
3679	zlevinter(ig) = zlev(ig, 1)
3680	END DO
3681	DO ig = 1, ngrid
3682	! calcul de zlevinter
3683	zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
3684	zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
3685	zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
3686	END DO
3687
3688	PRINT *, 'avant fermeture'
3689	! Fermeture,determination de f
3690	DO ig = 1, ngrid
3691	entr_star2(ig) = 0.
3692	END DO
3693	DO ig = 1, ngrid
3694	IF (entr_star_tot(ig)<1.E-10) THEN
3695	f(ig) = 0.
3696	ELSE
3697	DO k = lmin(ig), lentr(ig)
3698	entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)*2/(rho(ig,k)( &
3699	zlev(ig,k+1)-zlev(ig,k)))
3700	END DO
3701	! Nouvelle fermeture
3702	f(ig) = wmax(ig)/(max(500.,zmax(ig))r_aspectentr_star2(ig))* &
3703	entr_star_tot(ig)
3704	! test
3705	! if (first) then
3706	! f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
3707	! s *wmax(ig))
3708	! endif
3709	END IF
3710	! f0(ig)=f(ig)
3711	! first=.true.
3712	END DO
3713	PRINT *, 'apres fermeture'
3714
3715	! Calcul de l'entrainement
3716	DO k = 1, klev
3717	DO ig = 1, ngrid
3718	entr(ig, k) = f(ig)*entr_star(ig, k)
3719	END DO
3720	END DO
3721	! Calcul des flux
3722	DO ig = 1, ngrid
3723	DO l = 1, lmax(ig) - 1
3724	fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
3725	END DO
3726	END DO
3727
3728	! RC
3729
3730
3731	! print*,'9 OK convect8'
3732	! print*,'WA1 ',wa_moy
3733
3734	! determination de l'indice du debut de la mixed layer ou w decroit
3735
3736	! calcul de la largeur de chaque ascendance dans le cas conservatif.
3737	! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
3738	! d'une couche est égale à la hauteur de la couche alimentante.
3739	! La vitesse maximale dans l'ascendance est aussi prise comme estimation
3740	! de la vitesse d'entrainement horizontal dans la couche alimentante.
3741
3742	DO l = 2, nlay
3743	DO ig = 1, ngrid
3744	IF (l<=lmaxa(ig)) THEN
3745	zw = max(wa_moy(ig,l), 1.E-10)
3746	larg_cons(ig, l) = zmax(ig)r_aspectfmc(ig, l)/(rhobarz(ig,l)*zw)
3747	END IF
3748	END DO
3749	END DO
3750
3751	DO l = 2, nlay
3752	DO ig = 1, ngrid
3753	IF (l<=lmaxa(ig)) THEN
3754	! if (idetr.eq.0) then
3755	! cette option est finalement en dur.
3756	IF ((l_mix*zlev(ig,l))<0.) THEN
3757	PRINT , 'pb l_mixzlev<0'
3758	END IF
3759	larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
3760	! else if (idetr.eq.1) then
3761	! larg_detr(ig,l)=larg_cons(ig,l)
3762	! s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
3763	! else if (idetr.eq.2) then
3764	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3765	! s *sqrt(wa_moy(ig,l))
3766	! else if (idetr.eq.4) then
3767	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3768	! s *wa_moy(ig,l)
3769	! endif
3770	END IF
3771	END DO
3772	END DO
3773
3774	! print*,'10 OK convect8'
3775	! print*,'WA2 ',wa_moy
3776	! calcul de la fraction de la maille concernée par l'ascendance en tenant
3777	! compte de l'epluchage du thermique.
3778
3779	! CR def de zmix continu (profil parabolique des vitesses)
3780	DO ig = 1, ngrid
3781	IF (lmix(ig)>1.) THEN
3782	! test
3783	IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
3784	(zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
3785	zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
3786	(zlev(ig,lmix(ig)))))>1E-10) THEN
3787
3788	zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
3789	)2-(zlev(ig,lmix(ig)+1))2)-(zw2(ig,lmix(ig))-zw2(ig, &
3790	lmix(ig)+1))((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))/ &
3791	(2.((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))((zlev(ig,lmix(ig)))- &
3792	(zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
3793	zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3794	ELSE
3795	zmix(ig) = zlev(ig, lmix(ig))
3796	PRINT *, 'pb zmix'
3797	END IF
3798	ELSE
3799	zmix(ig) = 0.
3800	END IF
3801	! test
3802	IF ((zmax(ig)-zmix(ig))<0.) THEN
3803	zmix(ig) = 0.99*zmax(ig)
3804	! print*,'pb zmix>zmax'
3805	END IF
3806	END DO
3807
3808	! calcul du nouveau lmix correspondant
3809	DO ig = 1, ngrid
3810	DO l = 1, klev
3811	IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
3812	lmix(ig) = l
3813	END IF
3814	END DO
3815	END DO
3816
3817	DO l = 2, nlay
3818	DO ig = 1, ngrid
3819	IF (larg_cons(ig,l)>1.) THEN
3820	! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
3821	fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
3822	! test
3823	fraca(ig, l) = max(fraca(ig,l), 0.)
3824	fraca(ig, l) = min(fraca(ig,l), 0.5)
3825	fracd(ig, l) = 1. - fraca(ig, l)
3826	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3827	ELSE
3828	! wa_moy(ig,l)=0.
3829	fraca(ig, l) = 0.
3830	fracc(ig, l) = 0.
3831	fracd(ig, l) = 1.
3832	END IF
3833	END DO
3834	END DO
3835	! CR: calcul de fracazmix
3836	DO ig = 1, ngrid
3837	fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
3838	(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
3839	fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
3840	,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3841	END DO
3842
3843	DO l = 2, nlay
3844	DO ig = 1, ngrid
3845	IF (larg_cons(ig,l)>1.) THEN
3846	IF (l>lmix(ig)) THEN
3847	! test
3848	IF (zmax(ig)-zmix(ig)<1.E-10) THEN
3849	! print*,'pb xxx'
3850	xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
3851	ELSE
3852	xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3853	END IF
3854	IF (idetr==0) THEN
3855	fraca(ig, l) = fracazmix(ig)
3856	ELSE IF (idetr==1) THEN
3857	fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
3858	ELSE IF (idetr==2) THEN
3859	fraca(ig, l) = fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
3860	ELSE
3861	fraca(ig, l) = fracazmix(ig)xxx(ig, l)*2
3862	END IF
3863	! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3864	fraca(ig, l) = max(fraca(ig,l), 0.)
3865	fraca(ig, l) = min(fraca(ig,l), 0.5)
3866	fracd(ig, l) = 1. - fraca(ig, l)
3867	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3868	END IF
3869	END IF
3870	END DO
3871	END DO
3872
3873	PRINT *, 'fin calcul fraca'
3874	! print*,'11 OK convect8'
3875	! print*,'Ea3 ',wa_moy
3876	! ------------------------------------------------------------------
3877	! Calcul de fracd, wd
3878	! somme wa - wd = 0
3879	! ------------------------------------------------------------------
3880
3881
3882	DO ig = 1, ngrid
3883	fm(ig, 1) = 0.
3884	fm(ig, nlay+1) = 0.
3885	END DO
3886
3887	DO l = 2, nlay
3888	DO ig = 1, ngrid
3889	fm(ig, l) = fraca(ig, l)wa_moy(ig, l)rhobarz(ig, l)
3890	! CR:test
3891	IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
3892	fm(ig, l) = fm(ig, l-1)
3893	! write(1,*)'ajustement fm, l',l
3894	END IF
3895	! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3896	! RC
3897	END DO
3898	DO ig = 1, ngrid
3899	IF (fracd(ig,l)<0.1) THEN
3900	abort_message = 'fracd trop petit'
3901	CALL abort_physic(modname, abort_message, 1)
3902	ELSE
3903	! vitesse descendante "diagnostique"
3904	wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
3905	END IF
3906	END DO
3907	END DO
3908
3909	DO l = 1, nlay
3910	DO ig = 1, ngrid
3911	! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3912	masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
3913	END DO
3914	END DO
3915
3916	! print*,'12 OK convect8'
3917	! print*,'WA4 ',wa_moy
3918	! c------------------------------------------------------------------
3919	! calcul du transport vertical
3920	! ------------------------------------------------------------------
3921
3922	GO TO 4444
3923	! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3924	DO l = 2, nlay - 1
3925	DO ig = 1, ngrid
3926	IF (fm(ig,l+1)ptimestep>masse(ig,l) .AND. fm(ig,l+1)ptimestep>masse( &
3927	ig,l+1)) THEN
3928	! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
3929	! s ,fm(ig,l+1)*ptimestep
3930	! s ,' M=',masse(ig,l),masse(ig,l+1)
3931	END IF
3932	END DO
3933	END DO
3934
3935	DO l = 1, nlay
3936	DO ig = 1, ngrid
3937	IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
3938	! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
3939	! s ,entr(ig,l)*ptimestep
3940	! s ,' M=',masse(ig,l)
3941	END IF
3942	END DO
3943	END DO
3944
3945	DO l = 1, nlay
3946	DO ig = 1, ngrid
3947	IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
3948	! print*,'WARN!!! fm exagere ig=',ig,' l=',l
3949	! s ,' FM=',fm(ig,l)
3950	END IF
3951	IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
3952	! print*,'WARN!!! masse exagere ig=',ig,' l=',l
3953	! s ,' M=',masse(ig,l)
3954	! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3955	! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3956	! print*,'zlev(ig,l+1),zlev(ig,l)'
3957	! s ,zlev(ig,l+1),zlev(ig,l)
3958	! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3959	! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3960	END IF
3961	IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
3962	! print*,'WARN!!! entr exagere ig=',ig,' l=',l
3963	! s ,' E=',entr(ig,l)
3964	END IF
3965	END DO
3966	END DO
3967
3968	4444 CONTINUE
3969
3970	! CR:redefinition du entr
3971	DO l = 1, nlay
3972	DO ig = 1, ngrid
3973	detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
3974	IF (detr(ig,l)<0.) THEN
3975	entr(ig, l) = entr(ig, l) - detr(ig, l)
3976	detr(ig, l) = 0.
3977	! print*,'WARNING !!! detrainement negatif ',ig,l
3978	END IF
3979	END DO
3980	END DO
3981	! RC
3982	IF (w2di==1) THEN
3983	fm0 = fm0 + ptimestep*(fm-fm0)/tho
3984	entr0 = entr0 + ptimestep*(entr-entr0)/tho
3985	ELSE
3986	fm0 = fm
3987	entr0 = entr
3988	END IF
3989
3990	IF (1==1) THEN
3991	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
3992	zha)
3993	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
3994	zoa)
3995	ELSE
3996	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
3997	zdhadj, zha)
3998	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
3999	pdoadj, zoa)
4000	END IF
4001
4002	IF (1==0) THEN
4003	CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
4004	zu, zv, pduadj, pdvadj, zua, zva)
4005	ELSE
4006	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
4007	zua)
4008	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
4009	zva)
4010	END IF
4011
4012	DO l = 1, nlay
4013	DO ig = 1, ngrid
4014	zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
4015	zf2 = zf/(1.-zf)
4016	thetath2(ig, l) = zf2(zha(ig,l)-zh(ig,l))*2
4017	wth2(ig, l) = zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
4018	END DO
4019	END DO
4020
4021
4022
4023	! print*,'13 OK convect8'
4024	! print*,'WA5 ',wa_moy
4025	DO l = 1, nlay
4026	DO ig = 1, ngrid
4027	pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
4028	END DO
4029	END DO
4030
4031
4032	! do l=1,nlay
4033	! do ig=1,ngrid
4034	! if(abs(pdtadj(ig,l))*86400..gt.500.) then
4035	! print*,'WARN!!! ig=',ig,' l=',l
4036	! s ,' pdtadj=',pdtadj(ig,l)
4037	! endif
4038	! if(abs(pdoadj(ig,l))*86400..gt.1.) then
4039	! print*,'WARN!!! ig=',ig,' l=',l
4040	! s ,' pdoadj=',pdoadj(ig,l)
4041	! endif
4042	! enddo
4043	! enddo
4044
4045	! print*,'14 OK convect8'
4046	! ------------------------------------------------------------------
4047	! Calculs pour les sorties
4048	! ------------------------------------------------------------------
4049
4050	IF (sorties) THEN
4051	DO l = 1, nlay
4052	DO ig = 1, ngrid
4053	zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
4054	zld(ig, l) = fracd(ig, l)*zmax(ig)
4055	IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
4056	(1.-fracd(ig,l))
4057	END DO
4058	END DO
4059
4060	! deja fait
4061	! do l=1,nlay
4062	! do ig=1,ngrid
4063	! detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4064	! if (detr(ig,l).lt.0.) then
4065	! entr(ig,l)=entr(ig,l)-detr(ig,l)
4066	! detr(ig,l)=0.
4067	! print*,'WARNING !!! detrainement negatif ',ig,l
4068	! endif
4069	! enddo
4070	! enddo
4071
4072	! print*,'15 OK convect8'
4073
4074
4075	! #define und
4076	GO TO 123
4077	#ifdef und
4078	CALL writeg1d(1, nlay, wd, 'wd ', 'wd ')
4079	CALL writeg1d(1, nlay, zwa, 'wa ', 'wa ')
4080	CALL writeg1d(1, nlay, fracd, 'fracd ', 'fracd ')
4081	CALL writeg1d(1, nlay, fraca, 'fraca ', 'fraca ')
4082	CALL writeg1d(1, nlay, wa_moy, 'wam ', 'wam ')
4083	CALL writeg1d(1, nlay, zla, 'la ', 'la ')
4084	CALL writeg1d(1, nlay, zld, 'ld ', 'ld ')
4085	CALL writeg1d(1, nlay, pt, 'pt ', 'pt ')
4086	CALL writeg1d(1, nlay, zh, 'zh ', 'zh ')
4087	CALL writeg1d(1, nlay, zha, 'zha ', 'zha ')
4088	CALL writeg1d(1, nlay, zu, 'zu ', 'zu ')
4089	CALL writeg1d(1, nlay, zv, 'zv ', 'zv ')
4090	CALL writeg1d(1, nlay, zo, 'zo ', 'zo ')
4091	CALL writeg1d(1, nlay, wh, 'wh ', 'wh ')
4092	CALL writeg1d(1, nlay, wu, 'wu ', 'wu ')
4093	CALL writeg1d(1, nlay, wv, 'wv ', 'wv ')
4094	CALL writeg1d(1, nlay, wo, 'w15uo ', 'wXo ')
4095	CALL writeg1d(1, nlay, zdhadj, 'zdhadj ', 'zdhadj ')
4096	CALL writeg1d(1, nlay, pduadj, 'pduadj ', 'pduadj ')
4097	CALL writeg1d(1, nlay, pdvadj, 'pdvadj ', 'pdvadj ')
4098	CALL writeg1d(1, nlay, pdoadj, 'pdoadj ', 'pdoadj ')
4099	CALL writeg1d(1, nlay, entr, 'entr ', 'entr ')
4100	CALL writeg1d(1, nlay, detr, 'detr ', 'detr ')
4101	CALL writeg1d(1, nlay, fm, 'fm ', 'fm ')
4102
4103	CALL writeg1d(1, nlay, pdtadj, 'pdtadj ', 'pdtadj ')
4104	CALL writeg1d(1, nlay, pplay, 'pplay ', 'pplay ')
4105	CALL writeg1d(1, nlay, pplev, 'pplev ', 'pplev ')
4106
4107	! recalcul des flux en diagnostique...
4108	! print*,'PAS DE TEMPS ',ptimestep
4109	CALL dt2f(pplev, pplay, pt, pdtadj, wh)
4110	CALL writeg1d(1, nlay, wh, 'wh2 ', 'wh2 ')
4111	#endif
4112	123 CONTINUE
4113
4114	END IF
4115
4116	! if(wa_moy(1,4).gt.1.e-10) stop
4117
4118	! print*,'19 OK convect8'
4119	RETURN
4120	END SUBROUTINE thermcell
4121
4122	SUBROUTINE dqthermcell(ngrid, nlay, ptimestep, fm, entr, masse, q, dq, qa)
4123	USE dimphy
4124	IMPLICIT NONE
4125
4126	! =======================================================================
4127
4128	! Calcul du transport verticale dans la couche limite en presence
4129	! de "thermiques" explicitement representes
4130	! calcul du dq/dt une fois qu'on connait les ascendances
4131
4132	! =======================================================================
4133
4134	INTEGER ngrid, nlay
4135
4136	REAL ptimestep
4137	REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4138	REAL entr(ngrid, nlay)
4139	REAL q(ngrid, nlay)
4140	REAL dq(ngrid, nlay)
4141
4142	REAL qa(klon, klev), detr(klon, klev), wqd(klon, klev+1)
4143
4144	INTEGER ig, k
4145
4146	! calcul du detrainement
4147
4148	DO k = 1, nlay
4149	DO ig = 1, ngrid
4150	detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4151	! test
4152	IF (detr(ig,k)<0.) THEN
4153	entr(ig, k) = entr(ig, k) - detr(ig, k)
4154	detr(ig, k) = 0.
4155	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
4156	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
4157	END IF
4158	IF (fm(ig,k+1)<0.) THEN
4159	! print*,'fm2<0!!!'
4160	END IF
4161	IF (entr(ig,k)<0.) THEN
4162	! print*,'entr2<0!!!'
4163	END IF
4164	END DO
4165	END DO
4166
4167	! calcul de la valeur dans les ascendances
4168	DO ig = 1, ngrid
4169	qa(ig, 1) = q(ig, 1)
4170	END DO
4171
4172	DO k = 2, nlay
4173	DO ig = 1, ngrid
4174	IF ((fm(ig,k+1)+detr(ig,k))ptimestep>1.E-5masse(ig,k)) THEN
4175	qa(ig, k) = (fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k))/ &
4176	(fm(ig,k+1)+detr(ig,k))
4177	ELSE
4178	qa(ig, k) = q(ig, k)
4179	END IF
4180	IF (qa(ig,k)<0.) THEN
4181	! print*,'qa<0!!!'
4182	END IF
4183	IF (q(ig,k)<0.) THEN
4184	! print*,'q<0!!!'
4185	END IF
4186	END DO
4187	END DO
4188
4189	DO k = 2, nlay
4190	DO ig = 1, ngrid
4191	! wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
4192	wqd(ig, k) = fm(ig, k)*q(ig, k)
4193	IF (wqd(ig,k)<0.) THEN
4194	! print*,'wqd<0!!!'
4195	END IF
4196	END DO
4197	END DO
4198	DO ig = 1, ngrid
4199	wqd(ig, 1) = 0.
4200	wqd(ig, nlay+1) = 0.
4201	END DO
4202
4203	DO k = 1, nlay
4204	DO ig = 1, ngrid
4205	dq(ig, k) = (detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k)-wqd(ig,k)+wqd(ig,k+ &
4206	1))/masse(ig, k)
4207	! if (dq(ig,k).lt.0.) then
4208	! print*,'dq<0!!!'
4209	! endif
4210	END DO
4211	END DO
4212
4213	RETURN
4214	END SUBROUTINE dqthermcell
4215	SUBROUTINE dvthermcell(ngrid, nlay, ptimestep, fm, entr, masse, fraca, larga, &
4216	u, v, du, dv, ua, va)
4217	USE dimphy
4218	IMPLICIT NONE
4219
4220	! =======================================================================
4221
4222	! Calcul du transport verticale dans la couche limite en presence
4223	! de "thermiques" explicitement representes
4224	! calcul du dq/dt une fois qu'on connait les ascendances
4225
4226	! =======================================================================
4227
4228	INTEGER ngrid, nlay
4229
4230	REAL ptimestep
4231	REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4232	REAL fraca(ngrid, nlay+1)
4233	REAL larga(ngrid)
4234	REAL entr(ngrid, nlay)
4235	REAL u(ngrid, nlay)
4236	REAL ua(ngrid, nlay)
4237	REAL du(ngrid, nlay)
4238	REAL v(ngrid, nlay)
4239	REAL va(ngrid, nlay)
4240	REAL dv(ngrid, nlay)
4241
4242	REAL qa(klon, klev), detr(klon, klev)
4243	REAL wvd(klon, klev+1), wud(klon, klev+1)
4244	REAL gamma0, gamma(klon, klev+1)
4245	REAL dua, dva
4246	INTEGER iter
4247
4248	INTEGER ig, k
4249
4250	! calcul du detrainement
4251
4252	DO k = 1, nlay
4253	DO ig = 1, ngrid
4254	detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4255	END DO
4256	END DO
4257
4258	! calcul de la valeur dans les ascendances
4259	DO ig = 1, ngrid
4260	ua(ig, 1) = u(ig, 1)
4261	va(ig, 1) = v(ig, 1)
4262	END DO
4263
4264	DO k = 2, nlay
4265	DO ig = 1, ngrid
4266	IF ((fm(ig,k+1)+detr(ig,k))ptimestep>1.E-5masse(ig,k)) THEN
4267	! On itère sur la valeur du coeff de freinage.
4268	! gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4269	gamma0 = masse(ig, k)sqrt(0.5(fraca(ig,k+1)+fraca(ig, &
4270	k)))*0.5/larga(ig)
4271	! gamma0=0.
4272	! la première fois on multiplie le coefficient de freinage
4273	! par le module du vent dans la couche en dessous.
4274	dua = ua(ig, k-1) - u(ig, k-1)
4275	dva = va(ig, k-1) - v(ig, k-1)
4276	DO iter = 1, 5
4277	gamma(ig, k) = gamma0sqrt(dua2+dva*2)
4278	ua(ig, k) = (fm(ig,k)*ua(ig,k-1)+(entr(ig,k)+gamma(ig, &
4279	k))*u(ig,k))/(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4280	va(ig, k) = (fm(ig,k)*va(ig,k-1)+(entr(ig,k)+gamma(ig, &
4281	k))*v(ig,k))/(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4282	! print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4283	dua = ua(ig, k) - u(ig, k)
4284	dva = va(ig, k) - v(ig, k)
4285	END DO
4286	ELSE
4287	ua(ig, k) = u(ig, k)
4288	va(ig, k) = v(ig, k)
4289	gamma(ig, k) = 0.
4290	END IF
4291	END DO
4292	END DO
4293
4294	DO k = 2, nlay
4295	DO ig = 1, ngrid
4296	wud(ig, k) = fm(ig, k)*u(ig, k)
4297	wvd(ig, k) = fm(ig, k)*v(ig, k)
4298	END DO
4299	END DO
4300	DO ig = 1, ngrid
4301	wud(ig, 1) = 0.
4302	wud(ig, nlay+1) = 0.
4303	wvd(ig, 1) = 0.
4304	wvd(ig, nlay+1) = 0.
4305	END DO
4306
4307	DO k = 1, nlay
4308	DO ig = 1, ngrid
4309	du(ig, k) = ((detr(ig,k)+gamma(ig,k))*ua(ig,k)-(entr(ig,k)+gamma(ig, &
4310	k))*u(ig,k)-wud(ig,k)+wud(ig,k+1))/masse(ig, k)
4311	dv(ig, k) = ((detr(ig,k)+gamma(ig,k))*va(ig,k)-(entr(ig,k)+gamma(ig, &
4312	k))*v(ig,k)-wvd(ig,k)+wvd(ig,k+1))/masse(ig, k)
4313	END DO
4314	END DO
4315
4316	RETURN
4317	END SUBROUTINE dvthermcell
4318	SUBROUTINE dqthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, frac, q, dq, &
4319	qa)
4320	USE dimphy
4321	IMPLICIT NONE
4322
4323	! =======================================================================
4324
4325	! Calcul du transport verticale dans la couche limite en presence
4326	! de "thermiques" explicitement representes
4327	! calcul du dq/dt une fois qu'on connait les ascendances
4328
4329	! =======================================================================
4330
4331	INTEGER ngrid, nlay
4332
4333	REAL ptimestep
4334	REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4335	REAL entr(ngrid, nlay), frac(ngrid, nlay)
4336	REAL q(ngrid, nlay)
4337	REAL dq(ngrid, nlay)
4338
4339	REAL qa(klon, klev), detr(klon, klev), wqd(klon, klev+1)
4340	REAL qe(klon, klev), zf, zf2
4341
4342	INTEGER ig, k
4343
4344	! calcul du detrainement
4345
4346	DO k = 1, nlay
4347	DO ig = 1, ngrid
4348	detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4349	END DO
4350	END DO
4351
4352	! calcul de la valeur dans les ascendances
4353	DO ig = 1, ngrid
4354	qa(ig, 1) = q(ig, 1)
4355	qe(ig, 1) = q(ig, 1)
4356	END DO
4357
4358	DO k = 2, nlay
4359	DO ig = 1, ngrid
4360	IF ((fm(ig,k+1)+detr(ig,k))ptimestep>1.E-5masse(ig,k)) THEN
4361	zf = 0.5*(frac(ig,k)+frac(ig,k+1))
4362	zf2 = 1./(1.-zf)
4363	qa(ig, k) = (fm(ig,k)qa(ig,k-1)+zf2entr(ig,k)*q(ig,k))/ &
4364	(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2)
4365	qe(ig, k) = (q(ig,k)-zfqa(ig,k))zf2
4366	ELSE
4367	qa(ig, k) = q(ig, k)
4368	qe(ig, k) = q(ig, k)
4369	END IF
4370	END DO
4371	END DO
4372
4373	DO k = 2, nlay
4374	DO ig = 1, ngrid
4375	! wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
4376	wqd(ig, k) = fm(ig, k)*qe(ig, k)
4377	END DO
4378	END DO
4379	DO ig = 1, ngrid
4380	wqd(ig, 1) = 0.
4381	wqd(ig, nlay+1) = 0.
4382	END DO
4383
4384	DO k = 1, nlay
4385	DO ig = 1, ngrid
4386	dq(ig, k) = (detr(ig,k)qa(ig,k)-entr(ig,k)qe(ig,k)-wqd(ig,k)+wqd(ig,k &
4387	+1))/masse(ig, k)
4388	END DO
4389	END DO
4390
4391	RETURN
4392	END SUBROUTINE dqthermcell2
4393	SUBROUTINE dvthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, fraca, &
4394	larga, u, v, du, dv, ua, va)
4395	USE dimphy
4396	IMPLICIT NONE
4397
4398	! =======================================================================
4399
4400	! Calcul du transport verticale dans la couche limite en presence
4401	! de "thermiques" explicitement representes
4402	! calcul du dq/dt une fois qu'on connait les ascendances
4403
4404	! =======================================================================
4405
4406	INTEGER ngrid, nlay
4407
4408	REAL ptimestep
4409	REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4410	REAL fraca(ngrid, nlay+1)
4411	REAL larga(ngrid)
4412	REAL entr(ngrid, nlay)
4413	REAL u(ngrid, nlay)
4414	REAL ua(ngrid, nlay)
4415	REAL du(ngrid, nlay)
4416	REAL v(ngrid, nlay)
4417	REAL va(ngrid, nlay)
4418	REAL dv(ngrid, nlay)
4419
4420	REAL qa(klon, klev), detr(klon, klev), zf, zf2
4421	REAL wvd(klon, klev+1), wud(klon, klev+1)
4422	REAL gamma0, gamma(klon, klev+1)
4423	REAL ue(klon, klev), ve(klon, klev)
4424	REAL dua, dva
4425	INTEGER iter
4426
4427	INTEGER ig, k
4428
4429	! calcul du detrainement
4430
4431	DO k = 1, nlay
4432	DO ig = 1, ngrid
4433	detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4434	END DO
4435	END DO
4436
4437	! calcul de la valeur dans les ascendances
4438	DO ig = 1, ngrid
4439	ua(ig, 1) = u(ig, 1)
4440	va(ig, 1) = v(ig, 1)
4441	ue(ig, 1) = u(ig, 1)
4442	ve(ig, 1) = v(ig, 1)
4443	END DO
4444
4445	DO k = 2, nlay
4446	DO ig = 1, ngrid
4447	IF ((fm(ig,k+1)+detr(ig,k))ptimestep>1.E-5masse(ig,k)) THEN
4448	! On itère sur la valeur du coeff de freinage.
4449	! gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4450	gamma0 = masse(ig, k)sqrt(0.5(fraca(ig,k+1)+fraca(ig, &
4451	k)))0.5/larga(ig)1.
4452	! s *0.5
4453	! gamma0=0.
4454	zf = 0.5*(fraca(ig,k)+fraca(ig,k+1))
4455	zf = 0.
4456	zf2 = 1./(1.-zf)
4457	! la première fois on multiplie le coefficient de freinage
4458	! par le module du vent dans la couche en dessous.
4459	dua = ua(ig, k-1) - u(ig, k-1)
4460	dva = va(ig, k-1) - v(ig, k-1)
4461	DO iter = 1, 5
4462	! On choisit une relaxation lineaire.
4463	gamma(ig, k) = gamma0
4464	! On choisit une relaxation quadratique.
4465	gamma(ig, k) = gamma0sqrt(dua2+dva*2)
4466	ua(ig, k) = (fm(ig,k)ua(ig,k-1)+(zf2entr(ig,k)+gamma(ig, &
4467	k))u(ig,k))/(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zf*zf2+gamma(ig,k) &
4468	)
4469	va(ig, k) = (fm(ig,k)va(ig,k-1)+(zf2entr(ig,k)+gamma(ig, &
4470	k))v(ig,k))/(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zf*zf2+gamma(ig,k) &
4471	)
4472	! print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4473	dua = ua(ig, k) - u(ig, k)
4474	dva = va(ig, k) - v(ig, k)
4475	ue(ig, k) = (u(ig,k)-zfua(ig,k))zf2
4476	ve(ig, k) = (v(ig,k)-zfva(ig,k))zf2
4477	END DO
4478	ELSE
4479	ua(ig, k) = u(ig, k)
4480	va(ig, k) = v(ig, k)
4481	ue(ig, k) = u(ig, k)
4482	ve(ig, k) = v(ig, k)
4483	gamma(ig, k) = 0.
4484	END IF
4485	END DO
4486	END DO
4487
4488	DO k = 2, nlay
4489	DO ig = 1, ngrid
4490	wud(ig, k) = fm(ig, k)*ue(ig, k)
4491	wvd(ig, k) = fm(ig, k)*ve(ig, k)
4492	END DO
4493	END DO
4494	DO ig = 1, ngrid
4495	wud(ig, 1) = 0.
4496	wud(ig, nlay+1) = 0.
4497	wvd(ig, 1) = 0.
4498	wvd(ig, nlay+1) = 0.
4499	END DO
4500
4501	DO k = 1, nlay
4502	DO ig = 1, ngrid
4503	du(ig, k) = ((detr(ig,k)+gamma(ig,k))*ua(ig,k)-(entr(ig,k)+gamma(ig, &
4504	k))*ue(ig,k)-wud(ig,k)+wud(ig,k+1))/masse(ig, k)
4505	dv(ig, k) = ((detr(ig,k)+gamma(ig,k))*va(ig,k)-(entr(ig,k)+gamma(ig, &
4506	k))*ve(ig,k)-wvd(ig,k)+wvd(ig,k+1))/masse(ig, k)
4507	END DO
4508	END DO
4509
4510	RETURN
4511	END SUBROUTINE dvthermcell2
4512	SUBROUTINE thermcell_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, &
4513	pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
4514	! ,pu_therm,pv_therm
4515	, r_aspect, l_mix, w2di, tho)
4516
4517	USE dimphy
4518	IMPLICIT NONE
4519
4520	! =======================================================================
4521
4522	! Calcul du transport verticale dans la couche limite en presence
4523	! de "thermiques" explicitement representes
4524
4525	! Réécriture à partir d'un listing papier à Habas, le 14/02/00
4526
4527	! le thermique est supposé homogène et dissipé par mélange avec
4528	! son environnement. la longueur l_mix contrôle l'efficacité du
4529	! mélange
4530
4531	! Le calcul du transport des différentes espèces se fait en prenant
4532	! en compte:
4533	! 1. un flux de masse montant
4534	! 2. un flux de masse descendant
4535	! 3. un entrainement
4536	! 4. un detrainement
4537
4538	! =======================================================================
4539
4540	! -----------------------------------------------------------------------
4541	! declarations:
4542	! -------------
4543
4544	include "YOMCST.h"
4545
4546	! arguments:
4547	! ----------
4548
4549	INTEGER ngrid, nlay, w2di
4550	REAL tho
4551	REAL ptimestep, l_mix, r_aspect
4552	REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
4553	REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
4554	REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
4555	REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
4556	REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
4557	REAL pphi(ngrid, nlay)
4558
4559	INTEGER idetr
4560	SAVE idetr
4561	DATA idetr/3/
4562	!$OMP THREADPRIVATE(idetr)
4563
4564	! local:
4565	! ------
4566
4567	INTEGER ig, k, l, lmaxa(klon), lmix(klon)
4568	REAL zsortie1d(klon)
4569	! CR: on remplace lmax(klon,klev+1)
4570	INTEGER lmax(klon), lmin(klon), lentr(klon)
4571	REAL linter(klon)
4572	REAL zmix(klon), fracazmix(klon)
4573	! RC
4574	REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
4575
4576	REAL zlev(klon, klev+1), zlay(klon, klev)
4577	REAL zh(klon, klev), zdhadj(klon, klev)
4578	REAL ztv(klon, klev)
4579	REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
4580	REAL wh(klon, klev+1)
4581	REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
4582	REAL zla(klon, klev+1)
4583	REAL zwa(klon, klev+1)
4584	REAL zld(klon, klev+1)
4585	REAL zwd(klon, klev+1)
4586	REAL zsortie(klon, klev)
4587	REAL zva(klon, klev)
4588	REAL zua(klon, klev)
4589	REAL zoa(klon, klev)
4590
4591	REAL zha(klon, klev)
4592	REAL wa_moy(klon, klev+1)
4593	REAL fraca(klon, klev+1)
4594	REAL fracc(klon, klev+1)
4595	REAL zf, zf2
4596	REAL thetath2(klon, klev), wth2(klon, klev)
4597	! common/comtherm/thetath2,wth2
4598
4599	REAL count_time
4600	INTEGER ialt
4601
4602	LOGICAL sorties
4603	REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
4604	REAL zpspsk(klon, klev)
4605
4606	! real wmax(klon,klev),wmaxa(klon)
4607	REAL wmax(klon), wmaxa(klon)
4608	REAL wa(klon, klev, klev+1)
4609	REAL wd(klon, klev+1)
4610	REAL larg_part(klon, klev, klev+1)
4611	REAL fracd(klon, klev+1)
4612	REAL xxx(klon, klev+1)
4613	REAL larg_cons(klon, klev+1)
4614	REAL larg_detr(klon, klev+1)
4615	REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
4616	REAL pu_therm(klon, klev), pv_therm(klon, klev)
4617	REAL fm(klon, klev+1), entr(klon, klev)
4618	REAL fmc(klon, klev+1)
4619
4620	! CR:nouvelles variables
4621	REAL f_star(klon, klev+1), entr_star(klon, klev)
4622	REAL entr_star_tot(klon), entr_star2(klon)
4623	REAL f(klon), f0(klon)
4624	REAL zlevinter(klon)
4625	LOGICAL first
4626	DATA first/.FALSE./
4627	SAVE first
4628	!$OMP THREADPRIVATE(first)
4629	! RC
4630
4631	CHARACTER *2 str2
4632	CHARACTER *10 str10
4633
4634	CHARACTER (LEN=20) :: modname = 'thermcell_sec'
4635	CHARACTER (LEN=80) :: abort_message
4636
4637	LOGICAL vtest(klon), down
4638
4639	EXTERNAL scopy
4640
4641	INTEGER ncorrec, ll
4642	SAVE ncorrec
4643	DATA ncorrec/0/
4644	!$OMP THREADPRIVATE(ncorrec)
4645
4646
4647	! -----------------------------------------------------------------------
4648	! initialisation:
4649	! ---------------
4650
4651	sorties = .TRUE.
4652	IF (ngrid/=klon) THEN
4653	PRINT *
4654	PRINT *, 'STOP dans convadj'
4655	PRINT *, 'ngrid =', ngrid
4656	PRINT *, 'klon =', klon
4657	END IF
4658
4659	! -----------------------------------------------------------------------
4660	! incrementation eventuelle de tendances precedentes:
4661	! ---------------------------------------------------
4662
4663	! print*,'0 OK convect8'
4664
4665	DO l = 1, nlay
4666	DO ig = 1, ngrid
4667	zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
4668	zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
4669	zu(ig, l) = pu(ig, l)
4670	zv(ig, l) = pv(ig, l)
4671	zo(ig, l) = po(ig, l)
4672	ztv(ig, l) = zh(ig, l)(1.+0.61zo(ig,l))
4673	END DO
4674	END DO
4675
4676	! print*,'1 OK convect8'
4677	! --------------------
4678
4679
4680	! + + + + + + + + + + +
4681
4682
4683	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
4684	! wh,wt,wo ...
4685
4686	! + + + + + + + + + + + zh,zu,zv,zo,rho
4687
4688
4689	! -------------------- zlev(1)
4690	! \\\\\\\\\\\\\\\\\\\\
4691
4692
4693
4694	! -----------------------------------------------------------------------
4695	! Calcul des altitudes des couches
4696	! -----------------------------------------------------------------------
4697
4698	DO l = 2, nlay
4699	DO ig = 1, ngrid
4700	zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
4701	END DO
4702	END DO
4703	DO ig = 1, ngrid
4704	zlev(ig, 1) = 0.
4705	zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
4706	END DO
4707	DO l = 1, nlay
4708	DO ig = 1, ngrid
4709	zlay(ig, l) = pphi(ig, l)/rg
4710	END DO
4711	END DO
4712
4713	! print*,'2 OK convect8'
4714	! -----------------------------------------------------------------------
4715	! Calcul des densites
4716	! -----------------------------------------------------------------------
4717
4718	DO l = 1, nlay
4719	DO ig = 1, ngrid
4720	rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)rdzh(ig,l))
4721	END DO
4722	END DO
4723
4724	DO l = 2, nlay
4725	DO ig = 1, ngrid
4726	rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
4727	END DO
4728	END DO
4729
4730	DO k = 1, nlay
4731	DO l = 1, nlay + 1
4732	DO ig = 1, ngrid
4733	wa(ig, k, l) = 0.
4734	END DO
4735	END DO
4736	END DO
4737
4738	! print*,'3 OK convect8'
4739	! ------------------------------------------------------------------
4740	! Calcul de w2, quarre de w a partir de la cape
4741	! a partir de w2, on calcule wa, vitesse de l'ascendance
4742
4743	! ATTENTION: Dans cette version, pour cause d'economie de memoire,
4744	! w2 est stoke dans wa
4745
4746	! ATTENTION: dans convect8, on n'utilise le calcule des wa
4747	! independants par couches que pour calculer l'entrainement
4748	! a la base et la hauteur max de l'ascendance.
4749
4750	! Indicages:
4751	! l'ascendance provenant du niveau k traverse l'interface l avec
4752	! une vitesse wa(k,l).
4753
4754	! --------------------
4755
4756	! + + + + + + + + + +
4757
4758	! wa(k,l) ---- -------------------- l
4759	! /\
4760	! /\|\|\ + + + + + + + + + +
4761	! \|\|
4762	! \|\| --------------------
4763	! \|\|
4764	! \|\| + + + + + + + + + +
4765	! \|\|
4766	! \|\| --------------------
4767	! \|\|__
4768	! \|___ + + + + + + + + + + k
4769
4770	! --------------------
4771
4772
4773
4774	! ------------------------------------------------------------------
4775
4776	! CR: ponderation entrainement des couches instables
4777	! def des entr_star tels que entr=f*entr_star
4778	DO l = 1, klev
4779	DO ig = 1, ngrid
4780	entr_star(ig, l) = 0.
4781	END DO
4782	END DO
4783	! determination de la longueur de la couche d entrainement
4784	DO ig = 1, ngrid
4785	lentr(ig) = 1
4786	END DO
4787
4788	! on ne considere que les premieres couches instables
4789	DO k = nlay - 2, 1, -1
4790	DO ig = 1, ngrid
4791	IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
4792	lentr(ig) = k
4793	END IF
4794	END DO
4795	END DO
4796
4797	! determination du lmin: couche d ou provient le thermique
4798	DO ig = 1, ngrid
4799	lmin(ig) = 1
4800	END DO
4801	DO ig = 1, ngrid
4802	DO l = nlay, 2, -1
4803	IF (ztv(ig,l-1)>ztv(ig,l)) THEN
4804	lmin(ig) = l - 1
4805	END IF
4806	END DO
4807	END DO
4808
4809	! definition de l'entrainement des couches
4810	DO l = 1, klev - 1
4811	DO ig = 1, ngrid
4812	IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
4813	entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))** & ! s
4814	! (zlev(ig,l+1)-zlev(ig,l))
4815	sqrt(zlev(ig,l+1))
4816	END IF
4817	END DO
4818	END DO
4819	! pas de thermique si couche 1 stable
4820	DO ig = 1, ngrid
4821	IF (lmin(ig)>1) THEN
4822	DO l = 1, klev
4823	entr_star(ig, l) = 0.
4824	END DO
4825	END IF
4826	END DO
4827	! calcul de l entrainement total
4828	DO ig = 1, ngrid
4829	entr_star_tot(ig) = 0.
4830	END DO
4831	DO ig = 1, ngrid
4832	DO k = 1, klev
4833	entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
4834	END DO
4835	END DO
4836
4837	! print*,'fin calcul entr_star'
4838	DO k = 1, klev
4839	DO ig = 1, ngrid
4840	ztva(ig, k) = ztv(ig, k)
4841	END DO
4842	END DO
4843	! RC
4844	! print*,'7 OK convect8'
4845	DO k = 1, klev + 1
4846	DO ig = 1, ngrid
4847	zw2(ig, k) = 0.
4848	fmc(ig, k) = 0.
4849	! CR
4850	f_star(ig, k) = 0.
4851	! RC
4852	larg_cons(ig, k) = 0.
4853	larg_detr(ig, k) = 0.
4854	wa_moy(ig, k) = 0.
4855	END DO
4856	END DO
4857
4858	! print*,'8 OK convect8'
4859	DO ig = 1, ngrid
4860	linter(ig) = 1.
4861	lmaxa(ig) = 1
4862	lmix(ig) = 1
4863	wmaxa(ig) = 0.
4864	END DO
4865
4866	! CR:
4867	DO l = 1, nlay - 2
4868	DO ig = 1, ngrid
4869	IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
4870	zw2(ig,l)<1E-10) THEN
4871	f_star(ig, l+1) = entr_star(ig, l)
4872	! test:calcul de dteta
4873	zw2(ig, l+1) = 2.rg(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
4874	(zlev(ig,l+1)-zlev(ig,l))0.4pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
4875	larg_detr(ig, l) = 0.
4876	ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
4877	l)>1.E-10)) THEN
4878	f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
4879	ztva(ig, l) = (f_star(ig,l)ztva(ig,l-1)+entr_star(ig,l)ztv(ig,l))/ &
4880	f_star(ig, l+1)
4881	zw2(ig, l+1) = zw2(ig, l)(f_star(ig,l)/f_star(ig,l+1))*2 + &
4882	2.rg(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
4883	END IF
4884	! determination de zmax continu par interpolation lineaire
4885	IF (zw2(ig,l+1)<0.) THEN
4886	! test
4887	IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
4888	! print*,'pb linter'
4889	END IF
4890	linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
4891	ig,l))
4892	zw2(ig, l+1) = 0.
4893	lmaxa(ig) = l
4894	ELSE
4895	IF (zw2(ig,l+1)<0.) THEN
4896	! print*,'pb1 zw2<0'
4897	END IF
4898	wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
4899	END IF
4900	IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
4901	! lmix est le niveau de la couche ou w (wa_moy) est maximum
4902	lmix(ig) = l + 1
4903	wmaxa(ig) = wa_moy(ig, l+1)
4904	END IF
4905	END DO
4906	END DO
4907	! print*,'fin calcul zw2'
4908
4909	! Calcul de la couche correspondant a la hauteur du thermique
4910	DO ig = 1, ngrid
4911	lmax(ig) = lentr(ig)
4912	END DO
4913	DO ig = 1, ngrid
4914	DO l = nlay, lentr(ig) + 1, -1
4915	IF (zw2(ig,l)<=1.E-10) THEN
4916	lmax(ig) = l - 1
4917	END IF
4918	END DO
4919	END DO
4920	! pas de thermique si couche 1 stable
4921	DO ig = 1, ngrid
4922	IF (lmin(ig)>1) THEN
4923	lmax(ig) = 1
4924	lmin(ig) = 1
4925	END IF
4926	END DO
4927
4928	! Determination de zw2 max
4929	DO ig = 1, ngrid
4930	wmax(ig) = 0.
4931	END DO
4932
4933	DO l = 1, nlay
4934	DO ig = 1, ngrid
4935	IF (l<=lmax(ig)) THEN
4936	IF (zw2(ig,l)<0.) THEN
4937	! print*,'pb2 zw2<0'
4938	END IF
4939	zw2(ig, l) = sqrt(zw2(ig,l))
4940	wmax(ig) = max(wmax(ig), zw2(ig,l))
4941	ELSE
4942	zw2(ig, l) = 0.
4943	END IF
4944	END DO
4945	END DO
4946
4947	! Longueur caracteristique correspondant a la hauteur des thermiques.
4948	DO ig = 1, ngrid
4949	zmax(ig) = 0.
4950	zlevinter(ig) = zlev(ig, 1)
4951	END DO
4952	DO ig = 1, ngrid
4953	! calcul de zlevinter
4954	zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
4955	zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
4956	zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
4957	END DO
4958
4959	! print*,'avant fermeture'
4960	! Fermeture,determination de f
4961	DO ig = 1, ngrid
4962	entr_star2(ig) = 0.
4963	END DO
4964	DO ig = 1, ngrid
4965	IF (entr_star_tot(ig)<1.E-10) THEN
4966	f(ig) = 0.
4967	ELSE
4968	DO k = lmin(ig), lentr(ig)
4969	entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)*2/(rho(ig,k)( &
4970	zlev(ig,k+1)-zlev(ig,k)))
4971	END DO
4972	! Nouvelle fermeture
4973	f(ig) = wmax(ig)/(max(500.,zmax(ig))r_aspectentr_star2(ig))* &
4974	entr_star_tot(ig)
4975	! test
4976	! if (first) then
4977	! f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
4978	! s *wmax(ig))
4979	! endif
4980	END IF
4981	! f0(ig)=f(ig)
4982	! first=.true.
4983	END DO
4984	! print*,'apres fermeture'
4985
4986	! Calcul de l'entrainement
4987	DO k = 1, klev
4988	DO ig = 1, ngrid
4989	entr(ig, k) = f(ig)*entr_star(ig, k)
4990	END DO
4991	END DO
4992	! CR:test pour entrainer moins que la masse
4993	DO ig = 1, ngrid
4994	DO l = 1, lentr(ig)
4995	IF ((entr(ig,l)ptimestep)>(0.9masse(ig,l))) THEN
4996	entr(ig, l+1) = entr(ig, l+1) + entr(ig, l) - &
4997	0.9*masse(ig, l)/ptimestep
4998	entr(ig, l) = 0.9*masse(ig, l)/ptimestep
4999	END IF
5000	END DO
5001	END DO
5002	! CR: fin test
5003	! Calcul des flux
5004	DO ig = 1, ngrid
5005	DO l = 1, lmax(ig) - 1
5006	fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
5007	END DO
5008	END DO
5009
5010	! RC
5011
5012
5013	! print*,'9 OK convect8'
5014	! print*,'WA1 ',wa_moy
5015
5016	! determination de l'indice du debut de la mixed layer ou w decroit
5017
5018	! calcul de la largeur de chaque ascendance dans le cas conservatif.
5019	! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
5020	! d'une couche est égale à la hauteur de la couche alimentante.
5021	! La vitesse maximale dans l'ascendance est aussi prise comme estimation
5022	! de la vitesse d'entrainement horizontal dans la couche alimentante.
5023
5024	DO l = 2, nlay
5025	DO ig = 1, ngrid
5026	IF (l<=lmaxa(ig)) THEN
5027	zw = max(wa_moy(ig,l), 1.E-10)
5028	larg_cons(ig, l) = zmax(ig)r_aspectfmc(ig, l)/(rhobarz(ig,l)*zw)
5029	END IF
5030	END DO
5031	END DO
5032
5033	DO l = 2, nlay
5034	DO ig = 1, ngrid
5035	IF (l<=lmaxa(ig)) THEN
5036	! if (idetr.eq.0) then
5037	! cette option est finalement en dur.
5038	IF ((l_mix*zlev(ig,l))<0.) THEN
5039	! print,'pb l_mixzlev<0'
5040	END IF
5041	! CR: test: nouvelle def de lambda
5042	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5043	IF (zw2(ig,l)>1.E-10) THEN
5044	larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
5045	ELSE
5046	larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
5047	END IF
5048	! RC
5049	! else if (idetr.eq.1) then
5050	! larg_detr(ig,l)=larg_cons(ig,l)
5051	! s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
5052	! else if (idetr.eq.2) then
5053	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5054	! s *sqrt(wa_moy(ig,l))
5055	! else if (idetr.eq.4) then
5056	! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5057	! s *wa_moy(ig,l)
5058	! endif
5059	END IF
5060	END DO
5061	END DO
5062
5063	! print*,'10 OK convect8'
5064	! print*,'WA2 ',wa_moy
5065	! calcul de la fraction de la maille concernée par l'ascendance en tenant
5066	! compte de l'epluchage du thermique.
5067
5068	! CR def de zmix continu (profil parabolique des vitesses)
5069	DO ig = 1, ngrid
5070	IF (lmix(ig)>1.) THEN
5071	! test
5072	IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
5073	(zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
5074	zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
5075	(zlev(ig,lmix(ig)))))>1E-10) THEN
5076
5077	zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
5078	)2-(zlev(ig,lmix(ig)+1))2)-(zw2(ig,lmix(ig))-zw2(ig, &
5079	lmix(ig)+1))((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))/ &
5080	(2.((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))((zlev(ig,lmix(ig)))- &
5081	(zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
5082	zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
5083	ELSE
5084	zmix(ig) = zlev(ig, lmix(ig))
5085	! print*,'pb zmix'
5086	END IF
5087	ELSE
5088	zmix(ig) = 0.
5089	END IF
5090	! test
5091	IF ((zmax(ig)-zmix(ig))<0.) THEN
5092	zmix(ig) = 0.99*zmax(ig)
5093	! print*,'pb zmix>zmax'
5094	END IF
5095	END DO
5096
5097	! calcul du nouveau lmix correspondant
5098	DO ig = 1, ngrid
5099	DO l = 1, klev
5100	IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
5101	lmix(ig) = l
5102	END IF
5103	END DO
5104	END DO
5105
5106	DO l = 2, nlay
5107	DO ig = 1, ngrid
5108	IF (larg_cons(ig,l)>1.) THEN
5109	! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
5110	fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
5111	! test
5112	fraca(ig, l) = max(fraca(ig,l), 0.)
5113	fraca(ig, l) = min(fraca(ig,l), 0.5)
5114	fracd(ig, l) = 1. - fraca(ig, l)
5115	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
5116	ELSE
5117	! wa_moy(ig,l)=0.
5118	fraca(ig, l) = 0.
5119	fracc(ig, l) = 0.
5120	fracd(ig, l) = 1.
5121	END IF
5122	END DO
5123	END DO
5124	! CR: calcul de fracazmix
5125	DO ig = 1, ngrid
5126	fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
5127	(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
5128	fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
5129	,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
5130	END DO
5131
5132	DO l = 2, nlay
5133	DO ig = 1, ngrid
5134	IF (larg_cons(ig,l)>1.) THEN
5135	IF (l>lmix(ig)) THEN
5136	! test
5137	IF (zmax(ig)-zmix(ig)<1.E-10) THEN
5138	! print*,'pb xxx'
5139	xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
5140	ELSE
5141	xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
5142	END IF
5143	IF (idetr==0) THEN
5144	fraca(ig, l) = fracazmix(ig)
5145	ELSE IF (idetr==1) THEN
5146	fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
5147	ELSE IF (idetr==2) THEN
5148	fraca(ig, l) = fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
5149	ELSE
5150	fraca(ig, l) = fracazmix(ig)xxx(ig, l)*2
5151	END IF
5152	! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
5153	fraca(ig, l) = max(fraca(ig,l), 0.)
5154	fraca(ig, l) = min(fraca(ig,l), 0.5)
5155	fracd(ig, l) = 1. - fraca(ig, l)
5156	fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
5157	END IF
5158	END IF
5159	END DO
5160	END DO
5161
5162	! print*,'fin calcul fraca'
5163	! print*,'11 OK convect8'
5164	! print*,'Ea3 ',wa_moy
5165	! ------------------------------------------------------------------
5166	! Calcul de fracd, wd
5167	! somme wa - wd = 0
5168	! ------------------------------------------------------------------
5169
5170
5171	DO ig = 1, ngrid
5172	fm(ig, 1) = 0.
5173	fm(ig, nlay+1) = 0.
5174	END DO
5175
5176	DO l = 2, nlay
5177	DO ig = 1, ngrid
5178	fm(ig, l) = fraca(ig, l)wa_moy(ig, l)rhobarz(ig, l)
5179	! CR:test
5180	IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
5181	fm(ig, l) = fm(ig, l-1)
5182	! write(1,*)'ajustement fm, l',l
5183	END IF
5184	! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
5185	! RC
5186	END DO
5187	DO ig = 1, ngrid
5188	IF (fracd(ig,l)<0.1) THEN
5189	abort_message = 'fracd trop petit'
5190	CALL abort_physic(modname, abort_message, 1)
5191	ELSE
5192	! vitesse descendante "diagnostique"
5193	wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
5194	END IF
5195	END DO
5196	END DO
5197
5198	DO l = 1, nlay
5199	DO ig = 1, ngrid
5200	! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
5201	masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
5202	END DO
5203	END DO
5204
5205	! print*,'12 OK convect8'
5206	! print*,'WA4 ',wa_moy
5207	! c------------------------------------------------------------------
5208	! calcul du transport vertical
5209	! ------------------------------------------------------------------
5210
5211	GO TO 4444
5212	! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
5213	DO l = 2, nlay - 1
5214	DO ig = 1, ngrid
5215	IF (fm(ig,l+1)ptimestep>masse(ig,l) .AND. fm(ig,l+1)ptimestep>masse( &
5216	ig,l+1)) THEN
5217	! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
5218	! s ,fm(ig,l+1)*ptimestep
5219	! s ,' M=',masse(ig,l),masse(ig,l+1)
5220	END IF
5221	END DO
5222	END DO
5223
5224	DO l = 1, nlay
5225	DO ig = 1, ngrid
5226	IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
5227	! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
5228	! s ,entr(ig,l)*ptimestep
5229	! s ,' M=',masse(ig,l)
5230	END IF
5231	END DO
5232	END DO
5233
5234	DO l = 1, nlay
5235	DO ig = 1, ngrid
5236	IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
5237	! print*,'WARN!!! fm exagere ig=',ig,' l=',l
5238	! s ,' FM=',fm(ig,l)
5239	END IF
5240	IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
5241	! print*,'WARN!!! masse exagere ig=',ig,' l=',l
5242	! s ,' M=',masse(ig,l)
5243	! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
5244	! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
5245	! print*,'zlev(ig,l+1),zlev(ig,l)'
5246	! s ,zlev(ig,l+1),zlev(ig,l)
5247	! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
5248	! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
5249	END IF
5250	IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
5251	! print*,'WARN!!! entr exagere ig=',ig,' l=',l
5252	! s ,' E=',entr(ig,l)
5253	END IF
5254	END DO
5255	END DO
5256
5257	4444 CONTINUE
5258
5259	! CR:redefinition du entr
5260	DO l = 1, nlay
5261	DO ig = 1, ngrid
5262	detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
5263	IF (detr(ig,l)<0.) THEN
5264	entr(ig, l) = entr(ig, l) - detr(ig, l)
5265	detr(ig, l) = 0.
5266	! print*,'WARNING !!! detrainement negatif ',ig,l
5267	END IF
5268	END DO
5269	END DO
5270	! RC
5271	IF (w2di==1) THEN
5272	fm0 = fm0 + ptimestep*(fm-fm0)/tho
5273	entr0 = entr0 + ptimestep*(entr-entr0)/tho
5274	ELSE
5275	fm0 = fm
5276	entr0 = entr
5277	END IF
5278
5279	IF (1==1) THEN
5280	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
5281	zha)
5282	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
5283	zoa)
5284	ELSE
5285	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
5286	zdhadj, zha)
5287	CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
5288	pdoadj, zoa)
5289	END IF
5290
5291	IF (1==0) THEN
5292	CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
5293	zu, zv, pduadj, pdvadj, zua, zva)
5294	ELSE
5295	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
5296	zua)
5297	CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
5298	zva)
5299	END IF
5300
5301	DO l = 1, nlay
5302	DO ig = 1, ngrid
5303	zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
5304	zf2 = zf/(1.-zf)
5305	thetath2(ig, l) = zf2(zha(ig,l)-zh(ig,l))*2
5306	wth2(ig, l) = zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
5307	END DO
5308	END DO
5309
5310
5311
5312	! print*,'13 OK convect8'
5313	! print*,'WA5 ',wa_moy
5314	DO l = 1, nlay
5315	DO ig = 1, ngrid
5316	pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
5317	END DO
5318	END DO
5319
5320
5321	! do l=1,nlay
5322	! do ig=1,ngrid
5323	! if(abs(pdtadj(ig,l))*86400..gt.500.) then
5324	! print*,'WARN!!! ig=',ig,' l=',l
5325	! s ,' pdtadj=',pdtadj(ig,l)
5326	! endif
5327	! if(abs(pdoadj(ig,l))*86400..gt.1.) then
5328	! print*,'WARN!!! ig=',ig,' l=',l
5329	! s ,' pdoadj=',pdoadj(ig,l)
5330	! endif
5331	! enddo
5332	! enddo
5333
5334	! print*,'14 OK convect8'
5335	! ------------------------------------------------------------------
5336	! Calculs pour les sorties
5337	! ------------------------------------------------------------------
5338
5339	RETURN
5340	END SUBROUTINE thermcell_sec
5341

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

Download in other formats:

Original Format