Context Navigation

thermcell.F @ 898

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

Rev	Line
[878]	1	SUBROUTINE calcul_sec(ngrid,nlay,ptimestep
	2	s ,pplay,pplev,pphi,zlev
[542]	3	s ,pu,pv,pt,po
[878]	4	s ,zmax,wmax,zw2,lmix
[542]	5	c s ,pu_therm,pv_therm
	6	s ,r_aspect,l_mix,w2di,tho)
[878]	7
[542]	8	IMPLICIT NONE
	9
	10	c=======================================================================
	11	c
	12	c Calcul du transport verticale dans la couche limite en presence
	13	c de "thermiques" explicitement representes
	14	c
[878]	15	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
[542]	16	c
[878]	17	c le thermique est supposé homogène et dissipé par mélange avec
	18	c son environnement. la longueur l_mix contrôle l'efficacité du
	19	c mélange
[542]	20	c
[878]	21	c Le calcul du transport des différentes espèces se fait en prenant
[542]	22	c en compte:
	23	c 1. un flux de masse montant
	24	c 2. un flux de masse descendant
	25	c 3. un entrainement
	26	c 4. un detrainement
	27	c
	28	c=======================================================================
	29
	30	c-----------------------------------------------------------------------
	31	c declarations:
	32	c -------------
	33
[878]	34	#include "dimensions.h"
	35	#include "dimphy.h"
[542]	36	#include "YOMCST.h"
	37
	38	c arguments:
	39	c ----------
	40
	41	INTEGER ngrid,nlay,w2di,tho
	42	real ptimestep,l_mix,r_aspect
	43	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
	44	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
	45	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
	46	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
	47	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	48	real pphi(ngrid,nlay)
	49
	50	integer idetr
	51	save idetr
	52	data idetr/3/
[878]	53
[542]	54	c local:
	55	c ------
	56
	57	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
	58	real zsortie1d(klon)
	59	c CR: on remplace lmax(klon,klev+1)
	60	INTEGER lmax(klon),lmin(klon),lentr(klon)
	61	real linter(klon)
	62	real zmix(klon), fracazmix(klon)
	63	c RC
	64	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
	65
	66	real zlev(klon,klev+1),zlay(klon,klev)
	67	REAL zh(klon,klev),zdhadj(klon,klev)
	68	REAL ztv(klon,klev)
	69	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
	70	REAL wh(klon,klev+1)
	71	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
	72	real zla(klon,klev+1)
	73	real zwa(klon,klev+1)
	74	real zld(klon,klev+1)
	75	real zwd(klon,klev+1)
	76	real zsortie(klon,klev)
	77	real zva(klon,klev)
	78	real zua(klon,klev)
	79	real zoa(klon,klev)
	80
	81	real zha(klon,klev)
	82	real wa_moy(klon,klev+1)
	83	real fraca(klon,klev+1)
	84	real fracc(klon,klev+1)
	85	real zf,zf2
[878]	86	real thetath2(klon,klev),wth2(klon,klev)
	87	common/comtherm/thetath2,wth2
[542]	88
	89	real count_time
	90	integer isplit,nsplit,ialt
	91	parameter (nsplit=10)
	92	data isplit/0/
	93	save isplit
[878]	94
[542]	95	logical sorties
	96	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
	97	real zpspsk(klon,klev)
	98
	99	c real wmax(klon,klev),wmaxa(klon)
	100	real wmax(klon),wmaxa(klon)
	101	real wa(klon,klev,klev+1)
	102	real wd(klon,klev+1)
	103	real larg_part(klon,klev,klev+1)
	104	real fracd(klon,klev+1)
	105	real xxx(klon,klev+1)
	106	real larg_cons(klon,klev+1)
	107	real larg_detr(klon,klev+1)
	108	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
	109	real pu_therm(klon,klev),pv_therm(klon,klev)
	110	real fm(klon,klev+1),entr(klon,klev)
	111	real fmc(klon,klev+1)
[878]	112
[542]	113	cCR:nouvelles variables
	114	real f_star(klon,klev+1),entr_star(klon,klev)
	115	real entr_star_tot(klon),entr_star2(klon)
[878]	116	real zalim(klon)
	117	integer lalim(klon)
	118	real norme(klon)
[542]	119	real f(klon), f0(klon)
	120	real zlevinter(klon)
[878]	121	logical therm
[542]	122	logical first
	123	data first /.false./
	124	save first
	125	cRC
	126
	127	character*2 str2
	128	character*10 str10
	129
	130	LOGICAL vtest(klon),down
	131
	132	EXTERNAL SCOPY
	133
	134	integer ncorrec,ll
	135	save ncorrec
	136	data ncorrec/0/
[878]	137
[542]	138	c
	139	c-----------------------------------------------------------------------
	140	c initialisation:
	141	c ---------------
	142	c
	143	sorties=.true.
	144	IF(ngrid.NE.klon) THEN
	145	PRINT*
	146	PRINT*,'STOP dans convadj'
	147	PRINT*,'ngrid =',ngrid
	148	PRINT*,'klon =',klon
	149	ENDIF
	150	c
	151	c-----------------------------------------------------------------------
	152	c incrementation eventuelle de tendances precedentes:
	153	c ---------------------------------------------------
	154
[878]	155	c print*,'0 OK convect8'
[542]	156
	157	DO 1010 l=1,nlay
	158	DO 1015 ig=1,ngrid
	159	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
	160	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
	161	zu(ig,l)=pu(ig,l)
	162	zv(ig,l)=pv(ig,l)
	163	zo(ig,l)=po(ig,l)
	164	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
	165	1015 CONTINUE
	166	1010 CONTINUE
	167
[878]	168	c print*,'1 OK convect8'
[542]	169	c --------------------
	170	c
	171	c
	172	c + + + + + + + + + + +
	173	c
	174	c
	175	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	176	c wh,wt,wo ...
	177	c
	178	c + + + + + + + + + + + zh,zu,zv,zo,rho
	179	c
	180	c
	181	c -------------------- zlev(1)
	182	c \\\\\\\\\\\\\\\\\\\\
	183	c
	184	c
	185
	186	c-----------------------------------------------------------------------
	187	c Calcul des altitudes des couches
	188	c-----------------------------------------------------------------------
	189
	190	do l=2,nlay
	191	do ig=1,ngrid
	192	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
	193	enddo
	194	enddo
	195	do ig=1,ngrid
	196	zlev(ig,1)=0.
	197	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
	198	enddo
	199	do l=1,nlay
	200	do ig=1,ngrid
	201	zlay(ig,l)=pphi(ig,l)/RG
	202	enddo
	203	enddo
	204
	205	c print*,'2 OK convect8'
	206	c-----------------------------------------------------------------------
	207	c Calcul des densites
	208	c-----------------------------------------------------------------------
	209
	210	do l=1,nlay
	211	do ig=1,ngrid
	212	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
	213	enddo
	214	enddo
	215
	216	do l=2,nlay
	217	do ig=1,ngrid
	218	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
	219	enddo
	220	enddo
	221
	222	do k=1,nlay
	223	do l=1,nlay+1
	224	do ig=1,ngrid
	225	wa(ig,k,l)=0.
	226	enddo
	227	enddo
	228	enddo
	229
	230	c print*,'3 OK convect8'
	231	c------------------------------------------------------------------
	232	c Calcul de w2, quarre de w a partir de la cape
	233	c a partir de w2, on calcule wa, vitesse de l'ascendance
	234	c
	235	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
	236	c w2 est stoke dans wa
	237	c
	238	c ATTENTION: dans convect8, on n'utilise le calcule des wa
	239	c independants par couches que pour calculer l'entrainement
	240	c a la base et la hauteur max de l'ascendance.
	241	c
	242	c Indicages:
	243	c l'ascendance provenant du niveau k traverse l'interface l avec
	244	c une vitesse wa(k,l).
	245	c
	246	c --------------------
	247	c
	248	c + + + + + + + + + +
	249	c
	250	c wa(k,l) ---- -------------------- l
	251	c /\
	252	c /\|\|\ + + + + + + + + + +
	253	c \|\|
	254	c \|\| --------------------
	255	c \|\|
	256	c \|\| + + + + + + + + + +
	257	c \|\|
	258	c \|\| --------------------
	259	c \|\|__
	260	c \|___ + + + + + + + + + + k
	261	c
	262	c --------------------
	263	c
	264	c
	265	c
	266	c------------------------------------------------------------------
	267
	268	cCR: ponderation entrainement des couches instables
	269	cdef des entr_star tels que entr=f*entr_star
	270	do l=1,klev
	271	do ig=1,ngrid
	272	entr_star(ig,l)=0.
	273	enddo
	274	enddo
	275	c determination de la longueur de la couche d entrainement
	276	do ig=1,ngrid
	277	lentr(ig)=1
	278	enddo
	279
	280	con ne considere que les premieres couches instables
[878]	281	therm=.false.
[542]	282	do k=nlay-2,1,-1
	283	do ig=1,ngrid
	284	if (ztv(ig,k).gt.ztv(ig,k+1).and.
	285	s ztv(ig,k+1).le.ztv(ig,k+2)) then
[878]	286	lentr(ig)=k+1
	287	therm=.true.
[542]	288	endif
	289	enddo
	290	enddo
[878]	291	climitation de la valeur du lentr
	292	c do ig=1,ngrid
	293	c lentr(ig)=min(5,lentr(ig))
	294	c enddo
[542]	295	c determination du lmin: couche d ou provient le thermique
	296	do ig=1,ngrid
	297	lmin(ig)=1
	298	enddo
	299	do ig=1,ngrid
	300	do l=nlay,2,-1
	301	if (ztv(ig,l-1).gt.ztv(ig,l)) then
	302	lmin(ig)=l-1
	303	endif
	304	enddo
	305	enddo
[878]	306	cinitialisations
	307	do ig=1,ngrid
	308	zalim(ig)=0.
	309	norme(ig)=0.
	310	lalim(ig)=1
	311	enddo
	312	do k=1,klev-1
	313	do ig=1,ngrid
	314	zalim(ig)=zalim(ig)+zlev(ig,k)*MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
	315	s /(zlev(ig,k+1)-zlev(ig,k)))
	316	c s *(zlev(ig,k+1)-zlev(ig,k))
	317	norme(ig)=norme(ig)+MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
	318	s /(zlev(ig,k+1)-zlev(ig,k)))
	319	c s *(zlev(ig,k+1)-zlev(ig,k))
	320	enddo
	321	enddo
	322	do ig=1,ngrid
	323	if (norme(ig).gt.1.e-10) then
	324	zalim(ig)=max(10.*zalim(ig)/norme(ig),zlev(ig,2))
	325	c zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig)))
	326	endif
	327	enddo
	328	cdétermination du lalim correspondant
	329	do k=1,klev-1
	330	do ig=1,ngrid
	331	if ((zalim(ig).gt.zlev(ig,k)).and.(zalim(ig).le.zlev(ig,k+1)))
	332	s then
	333	lalim(ig)=k
	334	endif
	335	enddo
	336	enddo
[542]	337	c
	338	c definition de l'entrainement des couches
	339	do l=1,klev-1
	340	do ig=1,ngrid
	341	if (ztv(ig,l).gt.ztv(ig,l+1).and.
[878]	342	s l.ge.lmin(ig).and.l.lt.lentr(ig)) then
	343	entr_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)
	344	c s *(zlev(ig,l+1)-zlev(ig,l))
	345	s *sqrt(zlev(ig,l+1))
	346	cautre def
	347	c entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
	348	c s /zlev(ig,lentr(ig)+2)))**(3./2.)
[542]	349	endif
	350	enddo
	351	enddo
[878]	352	cnouveau test
	353	c if (therm) then
	354	do l=1,klev-1
	355	do ig=1,ngrid
	356	if (ztv(ig,l).gt.ztv(ig,l+1).and.
	357	s l.ge.lmin(ig).and.l.le.lalim(ig)
	358	s .and.zalim(ig).gt.1.e-10) then
	359	c if (l.le.lentr(ig)) then
	360	c entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
	361	c s /zalim(ig)))**(3./2.)
	362	c write(10,*)zlev(ig,l),entr_star(ig,l)
	363	endif
	364	enddo
	365	enddo
	366	c endif
	367	c pas de thermique si couche 1 stable
[542]	368	do ig=1,ngrid
	369	if (lmin(ig).gt.5) then
	370	do l=1,klev
	371	entr_star(ig,l)=0.
	372	enddo
	373	endif
	374	enddo
	375	c calcul de l entrainement total
	376	do ig=1,ngrid
	377	entr_star_tot(ig)=0.
	378	enddo
	379	do ig=1,ngrid
	380	do k=1,klev
	381	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
	382	enddo
	383	enddo
[878]	384	c Calcul entrainement normalise
	385	do ig=1,ngrid
	386	if (entr_star_tot(ig).gt.1.e-10) then
	387	c do l=1,lentr(ig)
	388	do l=1,klev
	389	cdef possibles pour entr_star: zdthetadz, dthetadz, zdtheta
	390	entr_star(ig,l)=entr_star(ig,l)/entr_star_tot(ig)
	391	enddo
	392	endif
	393	enddo
[542]	394	c
[878]	395	c print*,'fin calcul entr_star'
[542]	396	do k=1,klev
	397	do ig=1,ngrid
	398	ztva(ig,k)=ztv(ig,k)
	399	enddo
	400	enddo
	401	cRC
	402	c print*,'7 OK convect8'
	403	do k=1,klev+1
	404	do ig=1,ngrid
	405	zw2(ig,k)=0.
	406	fmc(ig,k)=0.
	407	cCR
	408	f_star(ig,k)=0.
	409	cRC
	410	larg_cons(ig,k)=0.
	411	larg_detr(ig,k)=0.
	412	wa_moy(ig,k)=0.
	413	enddo
	414	enddo
	415
	416	c print*,'8 OK convect8'
	417	do ig=1,ngrid
	418	linter(ig)=1.
	419	lmaxa(ig)=1
	420	lmix(ig)=1
	421	wmaxa(ig)=0.
	422	enddo
	423
	424	cCR:
	425	do l=1,nlay-2
	426	do ig=1,ngrid
	427	if (ztv(ig,l).gt.ztv(ig,l+1)
	428	s .and.entr_star(ig,l).gt.1.e-10
	429	s .and.zw2(ig,l).lt.1e-10) then
	430	f_star(ig,l+1)=entr_star(ig,l)
	431	ctest:calcul de dteta
	432	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
	433	s *(zlev(ig,l+1)-zlev(ig,l))
	434	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	435	larg_detr(ig,l)=0.
	436	else if ((zw2(ig,l).ge.1e-10).and.
	437	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
	438	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
	439	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
	440	s *ztv(ig,l))/f_star(ig,l+1)
	441	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
	442	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	443	s *(zlev(ig,l+1)-zlev(ig,l))
	444	endif
	445	c determination de zmax continu par interpolation lineaire
	446	if (zw2(ig,l+1).lt.0.) then
	447	ctest
	448	if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
[878]	449	c print*,'pb linter'
[542]	450	endif
	451	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
	452	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	453	zw2(ig,l+1)=0.
	454	lmaxa(ig)=l
	455	else
	456	if (zw2(ig,l+1).lt.0.) then
[878]	457	c print*,'pb1 zw2<0'
[542]	458	endif
	459	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	460	endif
	461	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	462	c lmix est le niveau de la couche ou w (wa_moy) est maximum
	463	lmix(ig)=l+1
	464	wmaxa(ig)=wa_moy(ig,l+1)
	465	endif
	466	enddo
	467	enddo
[878]	468	c print*,'fin calcul zw2'
[542]	469	c
	470	c Calcul de la couche correspondant a la hauteur du thermique
	471	do ig=1,ngrid
	472	lmax(ig)=lentr(ig)
[878]	473	c lmax(ig)=lalim(ig)
[542]	474	enddo
	475	do ig=1,ngrid
	476	do l=nlay,lentr(ig)+1,-1
[878]	477	c do l=nlay,lalim(ig)+1,-1
[542]	478	if (zw2(ig,l).le.1.e-10) then
	479	lmax(ig)=l-1
	480	endif
	481	enddo
	482	enddo
[878]	483	c pas de thermique si couche 1 stable
[542]	484	do ig=1,ngrid
	485	if (lmin(ig).gt.5) then
	486	lmax(ig)=1
	487	lmin(ig)=1
[878]	488	lentr(ig)=1
	489	lalim(ig)=1
[542]	490	endif
	491	enddo
	492	c
	493	c Determination de zw2 max
	494	do ig=1,ngrid
	495	wmax(ig)=0.
	496	enddo
	497
	498	do l=1,nlay
	499	do ig=1,ngrid
	500	if (l.le.lmax(ig)) then
	501	if (zw2(ig,l).lt.0.)then
[878]	502	c print*,'pb2 zw2<0'
[542]	503	endif
	504	zw2(ig,l)=sqrt(zw2(ig,l))
	505	wmax(ig)=max(wmax(ig),zw2(ig,l))
	506	else
	507	zw2(ig,l)=0.
	508	endif
	509	enddo
	510	enddo
	511
	512	c Longueur caracteristique correspondant a la hauteur des thermiques.
	513	do ig=1,ngrid
	514	zmax(ig)=0.
	515	zlevinter(ig)=zlev(ig,1)
	516	enddo
	517	do ig=1,ngrid
	518	c calcul de zlevinter
	519	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
	520	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
	521	s -zlev(ig,lmax(ig)))
	522	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
	523	enddo
[878]	524	do ig=1,ngrid
	525	c write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig)
	526	enddo
	527	con stoppe après les calculs de zmax et wmax
	528	RETURN
[542]	529
[878]	530	c print*,'avant fermeture'
[542]	531	c Fermeture,determination de f
[878]	532	cAttention! entrainement normalisé ou pas?
[542]	533	do ig=1,ngrid
	534	entr_star2(ig)=0.
	535	enddo
	536	do ig=1,ngrid
	537	if (entr_star_tot(ig).LT.1.e-10) then
	538	f(ig)=0.
	539	else
	540	do k=lmin(ig),lentr(ig)
[878]	541	c do k=lmin(ig),lalim(ig)
[542]	542	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
	543	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
	544	enddo
	545	c Nouvelle fermeture
	546	f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
[878]	547	s *entr_star2(ig))
	548	c s *entr_star_tot(ig)
[542]	549	ctest
	550	c if (first) then
[878]	551	f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
	552	s *wmax(ig))
[542]	553	c endif
	554	endif
[878]	555	f0(ig)=f(ig)
[542]	556	c first=.true.
	557	enddo
[878]	558	c print*,'apres fermeture'
	559	con stoppe après la fermeture
	560	RETURN
[542]	561	c Calcul de l'entrainement
	562	do k=1,klev
	563	do ig=1,ngrid
	564	entr(ig,k)=f(ig)*entr_star(ig,k)
	565	enddo
	566	enddo
[878]	567	con stoppe après le calcul de entr
	568	c RETURN
	569	cCR:test pour entrainer moins que la masse
	570	c do ig=1,ngrid
	571	c do l=1,lentr(ig)
	572	c if ((entr(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
	573	c entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
	574	c s -0.9*masse(ig,l)/ptimestep
	575	c entr(ig,l)=0.9*masse(ig,l)/ptimestep
	576	c endif
	577	c enddo
	578	c enddo
	579	cCR: fin test
[542]	580	c Calcul des flux
	581	do ig=1,ngrid
	582	do l=1,lmax(ig)-1
	583	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
	584	enddo
	585	enddo
	586
	587	cRC
	588
	589
	590	c print*,'9 OK convect8'
	591	c print*,'WA1 ',wa_moy
	592
	593	c determination de l'indice du debut de la mixed layer ou w decroit
	594
	595	c calcul de la largeur de chaque ascendance dans le cas conservatif.
	596	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
[878]	597	c d'une couche est égale à la hauteur de la couche alimentante.
[542]	598	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
	599	c de la vitesse d'entrainement horizontal dans la couche alimentante.
	600
	601	do l=2,nlay
	602	do ig=1,ngrid
	603	if (l.le.lmaxa(ig)) then
	604	zw=max(wa_moy(ig,l),1.e-10)
	605	larg_cons(ig,l)=zmax(ig)*r_aspect
	606	s fmc(ig,l)/(rhobarz(ig,l)zw)
	607	endif
	608	enddo
	609	enddo
	610
	611	do l=2,nlay
	612	do ig=1,ngrid
	613	if (l.le.lmaxa(ig)) then
	614	c if (idetr.eq.0) then
	615	c cette option est finalement en dur.
	616	if ((l_mix*zlev(ig,l)).lt.0.)then
[878]	617	c print,'pb l_mixzlev<0'
[542]	618	endif
[878]	619	cCR: test: nouvelle def de lambda
	620	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	621	if (zw2(ig,l).gt.1.e-10) then
	622	larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
	623	else
[542]	624	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
[878]	625	endif
	626	cRC
[542]	627	c else if (idetr.eq.1) then
	628	c larg_detr(ig,l)=larg_cons(ig,l)
	629	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
	630	c else if (idetr.eq.2) then
	631	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	632	c s *sqrt(wa_moy(ig,l))
	633	c else if (idetr.eq.4) then
	634	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	635	c s *wa_moy(ig,l)
	636	c endif
	637	endif
	638	enddo
	639	enddo
	640
	641	c print*,'10 OK convect8'
	642	c print*,'WA2 ',wa_moy
[878]	643	c calcul de la fraction de la maille concernée par l'ascendance en tenant
[542]	644	c compte de l'epluchage du thermique.
	645	c
	646	cCR def de zmix continu (profil parabolique des vitesses)
	647	do ig=1,ngrid
	648	if (lmix(ig).gt.1.) then
	649	c test
	650	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	651	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
	652	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	653	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
	654	s then
	655	c
	656	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	657	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
	658	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	659	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
	660	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	661	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
	662	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	663	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
	664	else
	665	zmix(ig)=zlev(ig,lmix(ig))
[878]	666	c print*,'pb zmix'
[542]	667	endif
	668	else
	669	zmix(ig)=0.
	670	endif
	671	ctest
	672	if ((zmax(ig)-zmix(ig)).lt.0.) then
	673	zmix(ig)=0.99*zmax(ig)
	674	c print*,'pb zmix>zmax'
	675	endif
	676	enddo
	677	c
	678	c calcul du nouveau lmix correspondant
	679	do ig=1,ngrid
	680	do l=1,klev
	681	if (zmix(ig).ge.zlev(ig,l).and.
	682	s zmix(ig).lt.zlev(ig,l+1)) then
	683	lmix(ig)=l
	684	endif
	685	enddo
	686	enddo
	687	c
	688	do l=2,nlay
	689	do ig=1,ngrid
	690	if(larg_cons(ig,l).gt.1.) then
	691	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
	692	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
	693	s /(r_aspect*zmax(ig))
	694	c test
	695	fraca(ig,l)=max(fraca(ig,l),0.)
	696	fraca(ig,l)=min(fraca(ig,l),0.5)
	697	fracd(ig,l)=1.-fraca(ig,l)
	698	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
	699	else
	700	c wa_moy(ig,l)=0.
	701	fraca(ig,l)=0.
	702	fracc(ig,l)=0.
	703	fracd(ig,l)=1.
	704	endif
	705	enddo
	706	enddo
	707	cCR: calcul de fracazmix
	708	do ig=1,ngrid
	709	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
	710	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
	711	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
	712	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
	713	enddo
	714	c
	715	do l=2,nlay
	716	do ig=1,ngrid
	717	if(larg_cons(ig,l).gt.1.) then
	718	if (l.gt.lmix(ig)) then
	719	ctest
	720	if (zmax(ig)-zmix(ig).lt.1.e-10) then
	721	c print*,'pb xxx'
	722	xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
	723	else
	724	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
	725	endif
	726	if (idetr.eq.0) then
	727	fraca(ig,l)=fracazmix(ig)
	728	else if (idetr.eq.1) then
	729	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
	730	else if (idetr.eq.2) then
	731	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
	732	else
	733	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
	734	endif
	735	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
	736	fraca(ig,l)=max(fraca(ig,l),0.)
	737	fraca(ig,l)=min(fraca(ig,l),0.5)
	738	fracd(ig,l)=1.-fraca(ig,l)
	739	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
	740	endif
	741	endif
	742	enddo
	743	enddo
	744
[878]	745	c print*,'fin calcul fraca'
[542]	746	c print*,'11 OK convect8'
	747	c print*,'Ea3 ',wa_moy
	748	c------------------------------------------------------------------
	749	c Calcul de fracd, wd
	750	c somme wa - wd = 0
	751	c------------------------------------------------------------------
	752
	753
	754	do ig=1,ngrid
	755	fm(ig,1)=0.
	756	fm(ig,nlay+1)=0.
	757	enddo
	758
	759	do l=2,nlay
	760	do ig=1,ngrid
	761	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
	762	cCR:test
	763	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
	764	s .and.l.gt.lmix(ig)) then
	765	fm(ig,l)=fm(ig,l-1)
	766	c write(1,*)'ajustement fm, l',l
	767	endif
	768	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
	769	cRC
	770	enddo
	771	do ig=1,ngrid
	772	if(fracd(ig,l).lt.0.1) then
	773	stop'fracd trop petit'
	774	else
	775	c vitesse descendante "diagnostique"
	776	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
	777	endif
	778	enddo
	779	enddo
	780
	781	do l=1,nlay
	782	do ig=1,ngrid
	783	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
	784	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
	785	enddo
	786	enddo
	787
[878]	788	c print*,'12 OK convect8'
[542]	789	c print*,'WA4 ',wa_moy
	790	cc------------------------------------------------------------------
	791	c calcul du transport vertical
	792	c------------------------------------------------------------------
	793
	794	go to 4444
	795	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
	796	do l=2,nlay-1
	797	do ig=1,ngrid
	798	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
	799	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
	800	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
	801	c s ,fm(ig,l+1)*ptimestep
	802	c s ,' M=',masse(ig,l),masse(ig,l+1)
	803	endif
	804	enddo
	805	enddo
	806
	807	do l=1,nlay
	808	do ig=1,ngrid
	809	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
	810	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
	811	c s ,entr(ig,l)*ptimestep
	812	c s ,' M=',masse(ig,l)
	813	endif
	814	enddo
	815	enddo
	816
	817	do l=1,nlay
	818	do ig=1,ngrid
	819	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
	820	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
	821	c s ,' FM=',fm(ig,l)
	822	endif
	823	if(.not.masse(ig,l).ge.1.e-10
	824	s .or..not.masse(ig,l).le.1.e4) then
	825	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
	826	c s ,' M=',masse(ig,l)
	827	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
	828	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
	829	c print*,'zlev(ig,l+1),zlev(ig,l)'
	830	c s ,zlev(ig,l+1),zlev(ig,l)
	831	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
	832	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
	833	endif
	834	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
	835	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
	836	c s ,' E=',entr(ig,l)
	837	endif
	838	enddo
	839	enddo
	840
	841	4444 continue
	842
	843	cCR:redefinition du entr
	844	do l=1,nlay
	845	do ig=1,ngrid
	846	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
	847	if (detr(ig,l).lt.0.) then
[878]	848	c entr(ig,l)=entr(ig,l)-detr(ig,l)
	849	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
[542]	850	detr(ig,l)=0.
	851	c print*,'WARNING !!! detrainement negatif ',ig,l
	852	endif
	853	enddo
	854	enddo
	855	cRC
	856	if (w2di.eq.1) then
	857	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
	858	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
	859	else
	860	fm0=fm
	861	entr0=entr
	862	endif
	863
	864	if (1.eq.1) then
	865	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	866	. ,zh,zdhadj,zha)
	867	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	868	. ,zo,pdoadj,zoa)
	869	else
	870	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
	871	. ,zh,zdhadj,zha)
	872	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
	873	. ,zo,pdoadj,zoa)
	874	endif
	875
	876	if (1.eq.0) then
	877	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
	878	. ,fraca,zmax
	879	. ,zu,zv,pduadj,pdvadj,zua,zva)
	880	else
	881	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	882	. ,zu,pduadj,zua)
	883	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	884	. ,zv,pdvadj,zva)
	885	endif
	886
	887	do l=1,nlay
	888	do ig=1,ngrid
	889	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
	890	zf2=zf/(1.-zf)
	891	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
	892	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
	893	enddo
	894	enddo
	895
	896
	897
	898	c print*,'13 OK convect8'
	899	c print*,'WA5 ',wa_moy
	900	do l=1,nlay
	901	do ig=1,ngrid
	902	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
	903	enddo
	904	enddo
	905
	906
	907	c do l=1,nlay
	908	c do ig=1,ngrid
	909	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
	910	c print*,'WARN!!! ig=',ig,' l=',l
	911	c s ,' pdtadj=',pdtadj(ig,l)
	912	c endif
	913	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
	914	c print*,'WARN!!! ig=',ig,' l=',l
	915	c s ,' pdoadj=',pdoadj(ig,l)
	916	c endif
	917	c enddo
	918	c enddo
	919
[878]	920	c print*,'14 OK convect8'
[542]	921	c------------------------------------------------------------------
	922	c Calculs pour les sorties
	923	c------------------------------------------------------------------
	924
	925	if(sorties) then
	926	do l=1,nlay
	927	do ig=1,ngrid
	928	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
	929	zld(ig,l)=fracd(ig,l)*zmax(ig)
	930	if(1.-fracd(ig,l).gt.1.e-10)
	931	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
	932	enddo
	933	enddo
	934
	935	cdeja fait
	936	c do l=1,nlay
	937	c do ig=1,ngrid
	938	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
	939	c if (detr(ig,l).lt.0.) then
	940	c entr(ig,l)=entr(ig,l)-detr(ig,l)
	941	c detr(ig,l)=0.
	942	c print*,'WARNING !!! detrainement negatif ',ig,l
	943	c endif
	944	c enddo
	945	c enddo
	946
	947	c print*,'15 OK convect8'
	948
	949	isplit=isplit+1
	950
	951
	952	c #define und
	953	goto 123
	954	#ifdef und
	955	CALL writeg1d(1,nlay,wd,'wd ','wd ')
	956	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
	957	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
	958	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
	959	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
	960	CALL writeg1d(1,nlay,zla,'la ','la ')
	961	CALL writeg1d(1,nlay,zld,'ld ','ld ')
	962	CALL writeg1d(1,nlay,pt,'pt ','pt ')
	963	CALL writeg1d(1,nlay,zh,'zh ','zh ')
	964	CALL writeg1d(1,nlay,zha,'zha ','zha ')
	965	CALL writeg1d(1,nlay,zu,'zu ','zu ')
	966	CALL writeg1d(1,nlay,zv,'zv ','zv ')
	967	CALL writeg1d(1,nlay,zo,'zo ','zo ')
	968	CALL writeg1d(1,nlay,wh,'wh ','wh ')
	969	CALL writeg1d(1,nlay,wu,'wu ','wu ')
	970	CALL writeg1d(1,nlay,wv,'wv ','wv ')
	971	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
	972	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
	973	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
	974	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
	975	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
	976	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
	977	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
	978	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
	979
	980	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
	981	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
	982	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
	983
	984	c recalcul des flux en diagnostique...
	985	c print*,'PAS DE TEMPS ',ptimestep
	986	call dt2F(pplev,pplay,pt,pdtadj,wh)
	987	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
	988	#endif
	989	123 continue
[878]	990	#define troisD
[542]	991	#ifdef troisD
	992	c if (sorties) then
	993	print*,'Debut des wrgradsfi'
	994
	995	c print*,'16 OK convect8'
	996	call wrgradsfi(1,nlay,wd,'wd ','wd ')
	997	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
	998	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
	999	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
	1000	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
	1001	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
	1002	c print*,'WA6 ',wa_moy
	1003	call wrgradsfi(1,nlay,zla,'la ','la ')
	1004	call wrgradsfi(1,nlay,zld,'ld ','ld ')
	1005	call wrgradsfi(1,nlay,pt,'pt ','pt ')
	1006	call wrgradsfi(1,nlay,zh,'zh ','zh ')
	1007	call wrgradsfi(1,nlay,zha,'zha ','zha ')
	1008	call wrgradsfi(1,nlay,zua,'zua ','zua ')
	1009	call wrgradsfi(1,nlay,zva,'zva ','zva ')
	1010	call wrgradsfi(1,nlay,zu,'zu ','zu ')
	1011	call wrgradsfi(1,nlay,zv,'zv ','zv ')
	1012	call wrgradsfi(1,nlay,zo,'zo ','zo ')
	1013	call wrgradsfi(1,nlay,wh,'wh ','wh ')
	1014	call wrgradsfi(1,nlay,wu,'wu ','wu ')
	1015	call wrgradsfi(1,nlay,wv,'wv ','wv ')
	1016	call wrgradsfi(1,nlay,wo,'wo ','wo ')
	1017	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
	1018	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
	1019	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
	1020	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
	1021	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
	1022	call wrgradsfi(1,nlay,entr,'entr ','entr ')
	1023	call wrgradsfi(1,nlay,detr,'detr ','detr ')
	1024	call wrgradsfi(1,nlay,fm,'fm ','fm ')
	1025	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
	1026	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
	1027	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
	1028	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
	1029
	1030	call wrgradsfi(1,nlay,zo,'zo ','zo ')
	1031	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
	1032	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
	1033
	1034	cCR:nouveaux diagnostiques
	1035	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
	1036	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
	1037	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
	1038	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
	1039	zsortie1d(:)=lmax(:)
	1040	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
	1041	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
	1042	zsortie1d(:)=lmix(:)
	1043	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
	1044	zsortie1d(:)=lentr(:)
	1045	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
	1046
	1047	c print*,'17 OK convect8'
	1048
	1049	do k=1,klev/10
	1050	write(str2,'(i2.2)') k
	1051	str10='wa'//str2
	1052	do l=1,nlay
	1053	do ig=1,ngrid
	1054	zsortie(ig,l)=wa(ig,k,l)
	1055	enddo
	1056	enddo
	1057	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
	1058	do l=1,nlay
	1059	do ig=1,ngrid
	1060	zsortie(ig,l)=larg_part(ig,k,l)
	1061	enddo
	1062	enddo
	1063	str10='la'//str2
	1064	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
	1065	enddo
	1066
	1067
	1068	c print*,'18 OK convect8'
	1069	c endif
	1070	print*,'Fin des wrgradsfi'
	1071	#endif
	1072
	1073	endif
	1074
	1075	c if(wa_moy(1,4).gt.1.e-10) stop
	1076
[878]	1077	c print*,'19 OK convect8'
[542]	1078	return
	1079	end
	1080
[878]	1081	SUBROUTINE fermeture_seche(ngrid,nlay
	1082	s ,pplay,pplev,pphi,zlev,rhobarz,f0,zpspsk
	1083	s ,alim_star,zh,zo,lentr,lmin,nu_min,nu_max,r_aspect
	1084	s ,zmax,wmax)
[542]	1085
[878]	1086	IMPLICIT NONE
[542]	1087
[878]	1088	#include "dimensions.h"
	1089	#include "dimphy.h"
	1090	#include "YOMCST.h"
[542]	1091
[878]	1092	INTEGER ngrid,nlay
	1093	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	1094	real pphi(ngrid,nlay)
	1095	real zlev(klon,klev+1)
	1096	real alim_star(klon,klev)
	1097	real f0(klon)
	1098	integer lentr(klon)
	1099	integer lmin(klon)
	1100	real zmax(klon)
	1101	real wmax(klon)
	1102	real nu_min
	1103	real nu_max
	1104	real r_aspect
	1105	real rhobarz(klon,klev+1)
	1106	REAL zh(klon,klev)
	1107	real zo(klon,klev)
	1108	real zpspsk(klon,klev)
[542]	1109
[878]	1110	integer ig,l
[542]	1111
[878]	1112	real f_star(klon,klev+1)
	1113	real detr_star(klon,klev)
	1114	real entr_star(klon,klev)
	1115	real zw2(klon,klev+1)
	1116	real linter(klon)
	1117	integer lmix(klon)
	1118	integer lmax(klon)
	1119	real zlevinter(klon)
	1120	real wa_moy(klon,klev+1)
	1121	real wmaxa(klon)
	1122	REAL ztv(klon,klev)
	1123	REAL ztva(klon,klev)
	1124	real nu(klon,klev)
	1125	real zmax0_sec(klon)
	1126	save zmax0_sec
[542]	1127
[878]	1128	do l=1,nlay
[542]	1129	do ig=1,ngrid
[878]	1130	ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
	1131	ztv(ig,l)=ztv(ig,l)(1.+RETVzo(ig,l))
[542]	1132	enddo
	1133	enddo
[878]	1134	do l=1,nlay-2
[542]	1135	do ig=1,ngrid
[878]	1136	if (ztv(ig,l).gt.ztv(ig,l+1)
	1137	s .and.alim_star(ig,l).gt.1.e-10
	1138	s .and.zw2(ig,l).lt.1e-10) then
	1139	f_star(ig,l+1)=alim_star(ig,l)
	1140	ctest:calcul de dteta
	1141	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
	1142	s *(zlev(ig,l+1)-zlev(ig,l))
	1143	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	1144	else if ((zw2(ig,l).ge.1e-10).and.
	1145	s (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
	1146	cestimation du detrainement a partir de la geometrie du pas precedent
	1147	ctests sur la definition du detr
	1148	nu(ig,l)=(nu_min+nu_max)/2.
	1149	s *(1.-(nu_max-nu_min)/(nu_max+nu_min)
	1150	s *tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005)))
	1151
	1152	detr_star(ig,l)=rhobarz(ig,l)
	1153	s *sqrt(zw2(ig,l))
	1154	s /(r_aspectzmax0_sec(ig))
	1155	c s /(r_aspectzmax0(ig))
	1156	s (sqrt(nu(ig,l)*zlev(ig,l+1)
	1157	s /sqrt(zw2(ig,l)))
	1158	s -sqrt(nu(ig,l)*zlev(ig,l)
	1159	s /sqrt(zw2(ig,l))))
	1160	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
	1161	if ((detr_star(ig,l)).gt.f_star(ig,l)) then
	1162	detr_star(ig,l)=f_star(ig,l)
	1163	endif
	1164	entr_star(ig,l)=0.9*detr_star(ig,l)
	1165	if ((l.lt.lentr(ig))) then
	1166	entr_star(ig,l)=0.
	1167	c detr_star(ig,l)=0.
	1168	endif
	1169	print*,'ok detr_star'
	1170	cprise en compte du detrainement dans le calcul du flux
	1171	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
	1172	s -detr_star(ig,l)
	1173	ctest sur le signe de f_star
	1174	if ((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10) then
	1175	cAM on melange Tl et qt du thermique
	1176	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig,l)
	1177	s +alim_star(ig,l))
	1178	s *ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
	1179	zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)
	1180	s /(f_star(ig,l+1)+detr_star(ig,l)))**2+
	1181	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	1182	s *(zlev(ig,l+1)-zlev(ig,l))
[542]	1183	endif
[878]	1184	endif
[542]	1185	c
[878]	1186	if (zw2(ig,l+1).lt.0.) then
	1187	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
	1188	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	1189	zw2(ig,l+1)=0.
	1190	print*,'linter=',linter(ig)
[542]	1191	else
[878]	1192	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
[542]	1193	endif
[878]	1194	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	1195	c lmix est le niveau de la couche ou w (wa_moy) est maximum
	1196	lmix(ig)=l+1
	1197	wmaxa(ig)=wa_moy(ig,l+1)
	1198	endif
[542]	1199	enddo
	1200	enddo
[878]	1201	print*,'fin calcul zw2'
	1202	c
	1203	c Calcul de la couche correspondant a la hauteur du thermique
[542]	1204	do ig=1,ngrid
[878]	1205	lmax(ig)=lentr(ig)
[542]	1206	enddo
	1207	do ig=1,ngrid
[878]	1208	do l=nlay,lentr(ig)+1,-1
	1209	if (zw2(ig,l).le.1.e-10) then
	1210	lmax(ig)=l-1
[542]	1211	endif
	1212	enddo
	1213	enddo
[878]	1214	c pas de thermique si couche 1 stable
[542]	1215	do ig=1,ngrid
[878]	1216	if (lmin(ig).gt.1) then
	1217	lmax(ig)=1
	1218	lmin(ig)=1
	1219	lentr(ig)=1
	1220	endif
	1221	enddo
	1222	c
	1223	c Determination de zw2 max
[542]	1224	do ig=1,ngrid
[878]	1225	wmax(ig)=0.
[542]	1226	enddo
	1227
[878]	1228	do l=1,nlay
[542]	1229	do ig=1,ngrid
[878]	1230	if (l.le.lmax(ig)) then
	1231	if (zw2(ig,l).lt.0.)then
	1232	print*,'pb2 zw2<0'
	1233	endif
	1234	zw2(ig,l)=sqrt(zw2(ig,l))
	1235	wmax(ig)=max(wmax(ig),zw2(ig,l))
[542]	1236	else
[878]	1237	zw2(ig,l)=0.
[542]	1238	endif
[878]	1239	enddo
[542]	1240	enddo
	1241
[878]	1242	c Longueur caracteristique correspondant a la hauteur des thermiques.
	1243	do ig=1,ngrid
	1244	zmax(ig)=0.
	1245	zlevinter(ig)=zlev(ig,1)
[542]	1246	enddo
[878]	1247	do ig=1,ngrid
	1248	c calcul de zlevinter
	1249	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
	1250	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
	1251	s -zlev(ig,lmax(ig)))
	1252	cpour le cas ou on prend tjs lmin=1
	1253	c zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
	1254	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
	1255	zmax0_sec(ig)=zmax(ig)
[542]	1256	enddo
[878]	1257
	1258	RETURN
	1259	END

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

Download in other formats:

Original Format