Context Navigation

thermcell.F @ 944

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

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

Download in other formats:

Original Format