Context Navigation

thermcell_old.F90 @ 5446

Last change on this file since 5446 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 149.6 KB

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

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

Download in other formats:

Original Format