Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

thermcell_old.F90

Last change on this file was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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

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

Download in other formats:

Original Format