Context Navigation

thermcell.F @ 1227

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

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

Download in other formats:

Original Format