Context Navigation

thermcell_old.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 190.2 KB

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

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

Download in other formats:

Original Format