Context Navigation

thermcell_old.F @ 901

Last change on this file since 901 was 878, checked in by Laurent Fairhead, 17 years ago
Bascule de la physique du LMD vers la physique avec thermiques LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 188.3 KB

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

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

Download in other formats:

Original Format