Context Navigation

thermcell.F90 @ 4794

Last change on this file since 4794 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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.5 KB

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

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

Download in other formats:

Original Format