Context Navigation

thermcell.F90 @ 2283

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

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

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

Download in other formats:

Original Format