Context Navigation

thermcell_old.F @ 1950

Last change on this file since 1950 was 1943, checked in by idelkadi, 10 years ago

Changes concerinng the Hourdin et al. 2002 version of the thermals and
Ayotte cases

calltherm.F90

call to thermcell_2002 modified :
1) iflag_thermals changed from 1 to >= 1000 ; iflag_thermals-1000
controls sub options.
2) thermals w and fraction in output files.

hbtm.F

Singularity in the 1/heat dependency of the Monin Obukov length L when
heat=0. Since 1/L is used rather than L, it is preferable to compute
directly L. There is a dependency in 1/u* then which is treated with a
threshold value.
(+ some cleaning in the syntax).

lmdz1d.F

If nday<0, -nday is used as the total number of time steps of the
simulation.
The option with imposed wtsurf and wqsurf read in lmdz1d.def was not
active anymore.
< IF(.NOT.ok_flux_surf) THEN
changed to

IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN

before

fsens=-wtsurf*rcpd*rho(1)
flat=-wqsurf*rlvtt*rho(1)

phys_output_write.F90 et phys_local_var_mod.F90

Removing the d_u_ajsb contribution to duthe (same for dv).
Those tendencies are not computed by the model ...
< zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
< d_u_ajsb(1:klon,1:klev)/pdtphys
---

zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys

! d_u_ajsb(1:klon,1:klev)/pdtphys

thermcell_dq.F90, thermcell_main.F90

Some cleaning

thermcell_old.F

Control by iflag_thermals >= 1000
wa_moy and fraca in outputs
+ cleaning

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

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

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

Download in other formats:

Original Format