Context Navigation

clmain.F @ 5444

Last change on this file since 5444 was 149, checked in by lmdz, 24 years ago
Continuation bug coefkz2 LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 57.1 KB

Rev	Line
[2]	1	SUBROUTINE clmain(dtime,pctsrf,t,q,u,v,
	2	. soil_model,ts,soilcap,soilflux,
	3	. paprs,pplay,radsol,snow,qsol,
	4	. xlat, rugos,
	5	. d_t,d_q,d_u,d_v,d_ts,
	6	. flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
	7	. rugmer, dflux_t,dflux_q,
	8	. zcoefh,zu1,zv1)
	9	cAA . itr, tr, flux_surf, d_tr)
	10	cAA REM:
	11	cAA-----
	12	cAA Tout ce qui a trait au traceurs est dans phytrac maintenant
	13	cAA pour l'instant le calcul de la couche limite pour les traceurs
	14	cAA se fait avec cltrac et ne tient pas compte de la differentiation
	15	cAA des sous-fraction de sol.
	16	cAA REM bis :
	17	cAA----------
	18	cAA Pour pouvoir extraire les coefficient d'echanges et le vent
	19	cAA dans la premiere couche, 3 champs supplementaires ont ete crees
	20	cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
	21	cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
	22	cAA si les informations des subsurfaces doivent etre prises en compte
	23	cAA il faudra sortir ces memes champs en leur ajoutant une dimension,
	24	cAA c'est a dire nbsrf (nbre de subsurface).
	25
	26	IMPLICIT none
	27	c======================================================================
	28	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	29	c Objet: interface de "couche limite" (diffusion verticale)
	30	c Arguments:
	31	c dtime----input-R- interval du temps (secondes)
	32	c t--------input-R- temperature (K)
	33	c q--------input-R- vapeur d'eau (kg/kg)
	34	c u--------input-R- vitesse u
	35	c v--------input-R- vitesse v
	36	c ts-------input-R- temperature du sol (en Kelvin)
	37	c paprs----input-R- pression a intercouche (Pa)
	38	c pplay----input-R- pression au milieu de couche (Pa)
	39	c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
	40	c capsol---input-R- inversion de l'effective capacite du sol (J/m2/K)
	41	c beta-----input-R- coefficient de l'evaporation reelle (0 a 1)
	42	c dif_grnd-input-R- coeff. de diffusion (chaleur) vers le sol profond
	43	c xlat-----input-R- latitude en degree
	44	c rugos----input-R- longeur de rugosite (en m)
	45	c
	46	c d_t------output-R- le changement pour "t"
	47	c d_q------output-R- le changement pour "q"
	48	c d_u------output-R- le changement pour "u"
	49	c d_v------output-R- le changement pour "v"
	50	c d_ts-----output-R- le changement pour "ts"
	51	c flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
	52	c (orientation positive vers le bas)
	53	c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
	54	c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
	55	c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
	56	c rugmer---output-R- longeur de rugosite sur mer (m)
	57	c dflux_t derive du flux sensible
	58	c dflux_q derive du flux latent
	59	cAA on rajoute en output yu1 et yv1 qui sont les vents dans
	60	cAA la premiere couche
	61	cAA ces 4 variables sont maintenant traites dans phytrac
	62	c itr--------input-I- nombre de traceurs
	63	c tr---------input-R- q. de traceurs
	64	c flux_surf--input-R- flux de traceurs a la surface
	65	c d_tr-------output-R tendance de traceurs
	66	c======================================================================
	67	#include "dimensions.h"
	68	#include "dimphy.h"
	69	#include "indicesol.h"
	70	c
	71	LOGICAL soil_model
	72	c
	73	REAL dtime
	74	REAL t(klon,klev), q(klon,klev)
	75	REAL u(klon,klev), v(klon,klev)
	76	REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
	77	REAL xlat(klon)
	78	REAL d_t(klon, klev), d_q(klon, klev)
	79	REAL d_u(klon, klev), d_v(klon, klev)
	80	REAL flux_t(klon,klev), flux_q(klon,klev)
	81	REAL dflux_t(klon), dflux_q(klon)
	82	REAL flux_u(klon,klev), flux_v(klon,klev)
	83	REAL rugmer(klon)
	84	REAL cdragh(klon), cdragm(klon)
	85	cAA INTEGER itr
	86	cAA REAL tr(klon,klev,nbtr)
	87	cAA REAL d_tr(klon,klev,nbtr)
	88	cAA REAL flux_surf(klon,nbtr)
	89	c
	90	REAL pctsrf(klon,nbsrf)
	91	REAL ts(klon,nbsrf)
	92	REAL d_ts(klon,nbsrf)
	93	REAL snow(klon,nbsrf)
	94	REAL qsol(klon,nbsrf)
	95	REAL rugos(klon,nbsrf)
	96	cAA
	97	REAL zcoefh(klon,klev)
	98	REAL zu1(klon)
	99	REAL zv1(klon)
	100	cAA
	101	c======================================================================
	102	EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
	103	c======================================================================
	104	REAL yts(klon), yrugos(klon), ypct(klon)
	105	REAL ycal(klon), ybeta(klon), ydif(klon)
	106	REAL yu1(klon), yv1(klon)
	107	REAL yrugm(klon), yrads(klon)
	108	REAL y_d_ts(klon)
	109	REAL y_d_t(klon, klev), y_d_q(klon, klev)
	110	REAL y_d_u(klon, klev), y_d_v(klon, klev)
	111	REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
	112	REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
	113	REAL y_dflux_t(klon), y_dflux_q(klon)
	114	REAL ycoefh(klon,klev), ycoefm(klon,klev)
	115	REAL ygamt(klon,2:klev) ! contre-gradient pour temperature
	116	REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
	117	REAL yu(klon,klev), yv(klon,klev)
	118	REAL yt(klon,klev), yq(klon,klev)
	119	REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
	120	cAA REAL ytr(klon,klev,nbtr)
	121	cAA REAL y_d_tr(klon,klev,nbtr)
	122	cAA REAL yflxsrf(klon,nbtr)
	123	c
[39]	124	LOGICAL contreg
	125	PARAMETER (contreg=.TRUE.)
	126	c
[2]	127	LOGICAL ok_nonloc
	128	PARAMETER (ok_nonloc=.FALSE.)
	129	REAL y_cd_h(klon), y_cd_m(klon)
	130	REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
	131	c
	132	#include "YOMCST.h"
	133	REAL u1lay(klon), v1lay(klon)
	134	REAL delp(klon,klev)
	135	REAL capsol(klon), beta(klon), dif_grnd(klon)
	136	REAL cal(klon)
	137	REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
	138	REAL totalflu(klon)
	139	INTEGER i, k, nsrf
	140	cAA INTEGER it
	141	INTEGER ni(klon), knon, j
	142	c======================================================================
	143	REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
	144	c======================================================================
	145	DO k = 1, klev ! epaisseur de couche
	146	DO i = 1, klon
	147	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	148	ENDDO
	149	ENDDO
	150	DO i = 1, klon ! vent de la premiere couche
	151	ccc zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
	152	zx_alf1 = 1.0
	153	zx_alf2 = 1.0 - zx_alf1
	154	u1lay(i) = u(i,1)zx_alf1 + u(i,2)zx_alf2
	155	v1lay(i) = v(i,1)zx_alf1 + v(i,2)zx_alf2
	156	ENDDO
	157	c
	158	c initialisation:
	159	c
	160	DO i = 1, klon
	161	rugmer(i) = 0.0
	162	cdragh(i) = 0.0
	163	cdragm(i) = 0.0
	164	dflux_t(i) = 0.0
	165	dflux_q(i) = 0.0
	166	zu1(i) = 0.0
	167	zv1(i) = 0.0
	168	ENDDO
	169	DO nsrf = 1, nbsrf
	170	DO i = 1, klon
	171	d_ts(i,nsrf) = 0.0
	172	ENDDO
	173	ENDDO
	174	DO k = 1, klev
	175	DO i = 1, klon
	176	d_t(i,k) = 0.0
	177	d_q(i,k) = 0.0
	178	flux_t(i,k) = 0.0
	179	flux_q(i,k) = 0.0
	180	d_u(i,k) = 0.0
	181	d_v(i,k) = 0.0
	182	flux_u(i,k) = 0.0
	183	flux_v(i,k) = 0.0
	184	zcoefh(i,k) = 0.0
	185	ENDDO
	186	ENDDO
	187	cAA IF (itr.GE.1) THEN
	188	cAA DO it = 1, itr
	189	cAA DO k = 1, klev
	190	cAA DO i = 1, klon
	191	cAA d_tr(i,k,it) = 0.0
	192	cAA ENDDO
	193	cAA ENDDO
	194	cAA ENDDO
	195	cAA ENDIF
	196	c
	197	c Boucler sur toutes les sous-fractions du sol:
	198	c
	199	DO 99999 nsrf = 1, nbsrf
	200	c
	201	c prescrire les proprietes du sol:
	202	CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
	203	IF (.NOT.soil_model) THEN
	204	DO i = 1, klon
	205	cal(i) = RCPD * capsol(i)
	206	totalflu(i) = radsol(i)
	207	ENDDO
	208	ELSE
	209	DO i = 1, klon
	210	totalflu(i) = soilflux(i,nsrf) + radsol(i)
	211	IF (nsrf.EQ.is_oce) THEN
	212	cal(i) = 0.0
	213	ELSE
	214	cal(i) = RCPD / soilcap(i,nsrf)
	215	ENDIF
	216	ENDDO
	217	ENDIF
	218	c
	219	c chercher les indices:
	220	DO j = 1, klon
	221	ni(j) = 0
	222	ENDDO
	223	knon = 0
	224	DO i = 1, klon
	225	IF (pctsrf(i,nsrf).GT.epsfra) THEN
	226	knon = knon + 1
	227	ni(knon) = i
	228	ENDIF
	229	ENDDO
	230	c
	231	IF (knon.EQ.0) GOTO 99999
	232	DO j = 1, knon
	233	i = ni(j)
	234	ypct(j) = pctsrf(i,nsrf)
	235	yts(j) = ts(i,nsrf)
	236	yrugos(j) = rugos(i,nsrf)
	237	yu1(j) = u1lay(i)
	238	yv1(j) = v1lay(i)
	239	yrads(j) = totalflu(i)
	240	ycal(j) = cal(i)
	241	ybeta(j) = beta(i)
	242	ydif(j) = dif_grnd(i)
	243	ypaprs(j,klev+1) = paprs(i,klev+1)
	244	ENDDO
	245	DO k = 1, klev
	246	DO j = 1, knon
	247	i = ni(j)
	248	ypaprs(j,k) = paprs(i,k)
	249	ypplay(j,k) = pplay(i,k)
	250	ydelp(j,k) = delp(i,k)
	251	yu(j,k) = u(i,k)
	252	yv(j,k) = v(i,k)
	253	yt(j,k) = t(i,k)
	254	yq(j,k) = q(i,k)
	255	ENDDO
	256	ENDDO
	257	c
	258	cAA IF (itr.GE.1) THEN
	259	cAA DO it = 1, itr
	260	cAA DO k = 1, klev
	261	cAA DO j = 1, knon
	262	cAA i = ni(j)
	263	cAA ytr(j,k,it) = tr(i,k,it)
	264	cAA ENDDO
	265	cAA ENDDO
	266	cAA DO j = 1, knon
	267	cAA i = ni(j)
	268	cAA yflxsrf(j,it) = flux_surf(i,it)
	269	cAA ENDDO
	270	cAA ENDDO
	271	cAA ENDIF
	272	c
	273	c calculer Cdrag et les coefficients d'echange
	274	CALL coefkz(nsrf, knon, ypaprs, ypplay,
	275	. yts, yrugos, yu, yv, yt, yq,
	276	. ycoefm, ycoefh)
[144]	277	CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
[39]	278	. ycoefm0, ycoefh0)
	279	DO k = 1, klev
	280	DO i = 1, knon
	281	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
	282	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
	283	ENDDO
	284	ENDDO
[2]	285	c
	286	c parametrisation non-locale:
	287	IF (ok_nonloc) THEN
	288	DO i = 1, knon
	289	y_cd_h(i) = ycoefh(i,1)
	290	y_cd_m(i) = ycoefm(i,1)
	291	ENDDO
	292	CALL nonlocal(knon, ypaprs, ypplay,
	293	. yts,ybeta,yu,yv,yt,yq,
	294	. y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
	295	DO k = 1, klev
	296	DO i = 1, knon
	297	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
	298	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
	299	ENDDO
	300	ENDDO
	301	ELSE
[39]	302	IF (.NOT.contreg) THEN
	303	DO k = 2, klev
	304	DO i = 1, knon
	305	ygamq(i,k) = 0.0
	306	ygamt(i,k) = 0.0
	307	ENDDO
	308	ENDDO
	309	ELSE
[2]	310	DO k = 3, klev
	311	DO i = 1, knon
	312	ygamq(i,k) = 0.0
	313	ygamt(i,k) = -1.0E-03
	314	ENDDO
	315	ENDDO
	316	DO i = 1, knon
	317	ygamq(i,2) = 0.0
	318	ygamt(i,2) = -2.5E-03
	319	ENDDO
[39]	320	ENDIF
[2]	321	ENDIF
	322	c
	323	c calculer la diffusion de "q" et de "h"
	324	CALL clqh(knon, dtime, yu1, yv1,
	325	e ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
	326	e ycal,ybeta,ydif,ygamt,ygamq,
	327	s y_d_t, y_d_q, y_d_ts,
	328	s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
	329	c
	330	c calculer la diffusion des vitesses "u" et "v"
	331	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
	332	s y_d_u,y_flux_u)
	333	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
	334	s y_d_v,y_flux_v)
	335	c
	336	c calculer la longueur de rugosite sur ocean
	337	IF (nsrf.EQ.is_oce) THEN
	338	DO j = 1, knon
	339	yrugm(j) = 0.018ycoefm(j,1) (yu1(j)2+yv1(j)2)/RG
	340	yrugm(j) = MAX(1.5e-05,yrugm(j))
	341	ENDDO
	342	ENDIF
	343	c
	344	cAA MAINTENANT DANS PHYTRAC
	345	cAAc calculer la diffusion des traceurs
	346	cAA IF (itr.GE.1) THEN
	347	cAA DO it = 1, itr
	348	cAA CALL cltrac(knon,dtime,ycoefh, yt, ytr(1,1,it), yflxsrf(1,it),
	349	cAA e ypaprs, ypplay, ydelp,
	350	cAA s y_d_tr(1,1,it))
	351	cAA ENDDO
	352	cAA ENDIF
	353	c
	354	DO j = 1, knon
	355	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
	356	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
	357	yu1(j) = yu1(j) * ypct(j)
	358	yv1(j) = yv1(j) * ypct(j)
	359	ENDDO
	360	c
	361	DO k = 1, klev
	362	DO j = 1, knon
	363	ycoefh(j,k) = ycoefh(j,k) * ypct(j)
	364	ycoefm(j,k) = ycoefm(j,k) * ypct(j)
	365	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
	366	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
	367	y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
	368	y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
	369	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
	370	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
	371	y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
	372	y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
	373	ENDDO
	374	ENDDO
	375	c
	376	DO j = 1, knon
	377	i = ni(j)
	378	d_ts(i,nsrf) = y_d_ts(j)
	379	rugmer(i) = yrugm(j)
	380	cdragh(i) = cdragh(i) + ycoefh(j,1)
	381	cdragm(i) = cdragm(i) + ycoefm(j,1)
	382	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
	383	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
	384	zu1(i) = zu1(i) + yu1(j)
	385	zv1(i) = zv1(i) + yv1(j)
	386	ENDDO
	387	c
	388	#ifdef CRAY
	389	DO k = 1, klev
	390	DO j = 1, knon
	391	i = ni(j)
	392	#else
	393	DO j = 1, knon
	394	i = ni(j)
	395	DO k = 1, klev
	396	#endif
	397	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
	398	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
	399	flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
	400	flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
	401	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
	402	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
	403	flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
	404	flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
	405	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
	406	ENDDO
	407	ENDDO
	408	c
	409	cAA IF (itr.GE.1) THEN
	410	cAA DO it = 1, itr
	411	cAA DO k = 1, klev
	412	cAA DO j = 1, knon
	413	cAA y_d_tr(j,k,it) = y_d_tr(j,k,it) * ypct(j)
	414	cAA ENDDO
	415	cAA ENDDO
	416	cAA ENDDO
	417	cAA DO j = 1, knon
	418	cAA i = ni(j)
	419	cAA DO it = 1, itr
	420	cAA DO k = 1, klev
	421	cAA d_tr(i,k,it) = d_tr(i,k,it) + y_d_tr(j,k,it)
	422	cAA ENDDO
	423	cAA ENDDO
	424	cAA ENDDO
	425	cAA ENDIF
	426	c
	427	99999 CONTINUE
	428	c
	429	RETURN
	430	END
	431	SUBROUTINE clqh(knon,dtime,u1lay,v1lay,coef,t,q,ts,paprs,pplay,
	432	e delp,radsol,cal,beta,dif_grnd, gamt,gamq,
	433	s d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
	434	IMPLICIT none
	435	c======================================================================
	436	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	437	c Objet: diffusion verticale de "q" et de "h"
	438	c======================================================================
	439	#include "dimensions.h"
	440	#include "dimphy.h"
	441	c Arguments:
	442	INTEGER knon
	443	REAL dtime ! intervalle du temps (s)
	444	REAL u1lay(klon) ! vitesse u de la 1ere couche (m/s)
	445	REAL v1lay(klon) ! vitesse v de la 1ere couche (m/s)
	446	REAL coef(klon,klev) ! le coefficient d'echange (m**2/s)
	447	c multiplie par le cisaillement du
	448	c vent (dV/dz); la premiere valeur
	449	c indique la valeur de Cdrag (sans unite)
	450	REAL cal(klon) ! Cp/cal indique la capacite calorifique
	451	c surfacique du sol
	452	REAL beta(klon) ! evap. reelle / evapotranspiration
	453	REAL dif_grnd(klon) ! coeff. diffusion vers le sol profond
	454	REAL t(klon,klev) ! temperature (K)
	455	REAL q(klon,klev) ! humidite specifique (kg/kg)
	456	REAL ts(klon) ! temperature du sol (K)
	457	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
	458	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
	459	REAL delp(klon,klev) ! epaisseur de couche en pression (Pa)
	460	REAL radsol(klon) ! ray. net au sol (Solaire+IR) W/m2
	461	c
	462	REAL d_t(klon,klev) ! incrementation de "t"
	463	REAL d_q(klon,klev) ! incrementation de "q"
	464	REAL d_ts(klon) ! incrementation de "ts"
	465	REAL flux_t(klon,klev) ! (diagnostic) flux de la chaleur
	466	c sensible, flux de Cp*T, positif vers
	467	c le bas: j/(m**2 s) c.a.d.: W/m2
	468	REAL flux_q(klon,klev) ! flux de la vapeur d'eau:kg/(m**2 s)
	469	REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
	470	REAL dflux_l(klon) ! derivee du flux latent dF/dTs
	471	REAL dflux_g(klon) ! derivee du flux du sol profond dF/dTs
	472	c======================================================================
	473	REAL t_grnd ! temperature de rappel pour glace de mer
	474	PARAMETER (t_grnd=271.35)
	475	REAL t_coup
	476	PARAMETER(t_coup=273.15)
	477	c======================================================================
	478	INTEGER i, k
	479	REAL zx_a
	480	REAL zx_b
	481	REAL zx_qs
	482	REAL zx_dq_s_dh
	483	REAL zx_h_grnd
	484	REAL zx_cq0(klon)
	485	REAL zx_dq0(klon)
	486	REAL zx_cq(klon,klev)
	487	REAL zx_dq(klon,klev)
	488	REAL zx_ch(klon,klev)
	489	REAL zx_dh(klon,klev)
	490	REAL zx_buf1(klon)
	491	REAL zx_buf2(klon)
	492	REAL zx_coef(klon,klev)
	493	REAL zx_q_0(klon)
	494	REAL zx_h_ts(klon)
	495	REAL zx_sl(klon)
	496	REAL local_h(klon,klev) ! enthalpie potentielle
	497	REAL local_q(klon,klev)
	498	REAL local_ts(klon)
	499	REAL psref(klon) ! pression de reference pour temperature potent.
	500	REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
	501	c======================================================================
	502	c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
	503	REAL gamq(klon,2:klev)
	504	c contre-gradient pour la chaleur sensible: Kelvin/metre
	505	REAL gamt(klon,2:klev)
	506	REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
	507	REAL zdelz
	508	c======================================================================
[77]	509	C Variables intermediaires pour le calcul des fluxs a la surface
	510	real zx_mh(klon), zx_nh(klon), zx_oh(klon)
	511	real zx_mq(klon), zx_nq(klon), zx_oq(klon)
	512	real zx_k1(klon), zx_dq_s_dt(klon)
	513	real zx_qsat(klon)
	514	c======================================================================
[2]	515	REAL zcor, zdelta, zcvm5
	516	#include "YOMCST.h"
	517	#include "YOETHF.h"
	518	#include "FCTTRE.h"
	519	c======================================================================
	520	DO i = 1, knon
	521	psref(i) = paprs(i,1) !pression de reference est celle au sol
	522	local_ts(i) = ts(i)
	523	ENDDO
	524	DO k = 1, klev
	525	DO i = 1, knon
	526	zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
	527	zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
	528	local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
	529	local_q(i,k) = q(i,k)
	530	ENDDO
	531	ENDDO
	532	c
	533	c Convertir les coefficients en variables convenables au calcul:
	534	c
	535	DO i = 1, knon
	536	zx_coef(i,1) = coef(i,1)
	537	. * (1.0+SQRT(u1lay(i)2+v1lay(i)2))
	538	. * pplay(i,1)/(RD*t(i,1))
	539	ENDDO
	540	c
	541	DO k = 2, klev
	542	DO i = 1, knon
	543	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
	544	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
	545	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
	546	ENDDO
	547	ENDDO
	548	c
	549	c Preparer les flux lies aux contre-gardients
	550	c
	551	DO k = 2, klev
	552	DO i = 1, knon
	553	zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
	554	. *(pplay(i,k-1)-pplay(i,k))
	555	z_gamaq(i,k) = gamq(i,k) * zdelz
	556	z_gamah(i,k) = gamt(i,k) * zdelz RCPD zx_pkh(i,k)
	557	ENDDO
	558	ENDDO
	559	c======================================================================
[77]	560	c
	561	c zx_qs = qsat en kg/kg
	562	c
[2]	563	DO i = 1, knon
	564	IF (thermcep) THEN
	565	zdelta=MAX(0.,SIGN(1.,rtt-local_ts(i)))
	566	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	567	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
	568	zx_qs= r2es * FOEEW(local_ts(i),zdelta)/paprs(i,1)
	569	zx_qs=MIN(0.5,zx_qs)
	570	zcor=1./(1.-retv*zx_qs)
	571	zx_qs=zx_qs*zcor
	572	zx_dq_s_dh = FOEDE(local_ts(i),zdelta,zcvm5,zx_qs,zcor)
	573	. /RLVTT / zx_pkh(i,1)
	574	ELSE
	575	IF (local_ts(i).LT.t_coup) THEN
	576	zx_qs = qsats(local_ts(i)) / paprs(i,1)
	577	zx_dq_s_dh = dqsats(local_ts(i),zx_qs)/RLVTT
	578	. / zx_pkh(i,1)
	579	ELSE
	580	zx_qs = qsatl(local_ts(i)) / paprs(i,1)
	581	zx_dq_s_dh = dqsatl(local_ts(i),zx_qs)/RLVTT
	582	. / zx_pkh(i,1)
	583	ENDIF
	584	ENDIF
[77]	585
	586	c zx_dq_s_dh = 0.
	587	zx_dq_s_dt(i) = RCPD * zx_pkh(i,1) * zx_dq_s_dh
	588	zx_qsat(i) = zx_qs
[2]	589	c
	590	ENDDO
	591	DO i = 1, knon
	592	zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
	593	zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
	594	. -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
	595	zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
	596	c
	597	zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
	598	zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
	599	. -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
	600	zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
	601	ENDDO
	602	DO k = klev-1, 2 , -1
	603	DO i = 1, knon
	604	zx_buf1(i) = delp(i,k)+zx_coef(i,k)
	605	. +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
	606	zx_cq(i,k) = (local_q(i,k)*delp(i,k)
	607	. +zx_coef(i,k+1)*zx_cq(i,k+1)
	608	. +zx_coef(i,k+1)*z_gamaq(i,k+1)
	609	. -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
	610	zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
	611	c
	612	zx_buf2(i) = delp(i,k)+zx_coef(i,k)
	613	. +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
	614	zx_ch(i,k) = (local_h(i,k)*delp(i,k)
	615	. +zx_coef(i,k+1)*zx_ch(i,k+1)
	616	. +zx_coef(i,k+1)*z_gamah(i,k+1)
	617	. -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
	618	zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
	619	ENDDO
	620	ENDDO
[77]	621	C
	622	C nouvelle formulation JL Dufresne
	623	C
	624	C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
	625	C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
	626	C
[2]	627	DO i = 1, knon
	628	zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
	629	zx_cq(i,1) = (local_q(i,1)*delp(i,1)
[77]	630	. +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
[2]	631	. /zx_buf1(i)
[77]	632	zx_dq(i,1) = -1. * RG / zx_buf1(i)
[2]	633	c
	634	zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
	635	zx_ch(i,1) = (local_h(i,1)*delp(i,1)
[77]	636	. +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
[2]	637	. /zx_buf2(i)
[77]	638	zx_dh(i,1) = -1. * RG / zx_buf2(i)
[2]	639	ENDDO
[77]	640
	641	C === Calcul de la temperature de surface ===
	642	C
	643	C zx_sl = chaleur latente d'evaporation ou de sublimation
	644	C
	645	do i = 1, knon
	646	zx_sl(i) = RLVTT
	647	if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT
	648	zx_k1(i) = zx_coef(i,1)
	649	enddo
	650	do i = 1, knon
	651	C Q
	652	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * zx_dq(i,1) * dtime)
	653	zx_mq(i) = beta(i) * zx_k1(i) *
	654	. (zx_cq(i,1) - zx_qsat(i)
	655	. + zx_dq_s_dt(i) * local_ts(i))
	656	. / zx_oq(i)
	657	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i))
	658	. / zx_oq(i)
	659	c H
	660	zx_oh(i) = 1. - (zx_k1(i) * zx_dh(i,1) * dtime)
	661	zx_mh(i) = zx_k1(i) * zx_ch(i,1) / zx_oh(i)
	662	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i,1))/ zx_oh(i)
	663	c Tsurface
	664	local_ts(i) = (ts(i) + cal(i)/(RCPD * zx_pkh(i,1)) * dtime *
	665	. (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i))
	666	. + dif_grnd(i) * t_grnd * dtime)/
	667	. ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i,1)) * (
	668	. zx_nh(i) + zx_sl(i) * zx_nq(i))
	669	. + dtime * dif_grnd(i))
	670	zx_h_ts(i) = local_ts(i) * RCPD * zx_pkh(i,1)
	671	d_ts(i) = local_ts(i) - ts(i)
	672	zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
	673	c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
	674	c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
	675	flux_q(i,1) = zx_mq(i) + zx_nq(i) * local_ts(i)
	676	flux_t(i,1) = zx_mh(i) + zx_nh(i) * local_ts(i)
	677	c Derives des flux dF/dTs (W m-2 K-1):
	678	dflux_s(i) = zx_nh(i)
	679	dflux_l(i) = (zx_sl(i) * zx_nq(i))
[2]	680	ENDDO
[77]	681
	682
[2]	683	c==== une fois on a zx_h_ts, on peut faire l'iteration ========
	684	DO i = 1, knon
[104]	685	local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)flux_t(i,1)dtime
	686	local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)flux_q(i,1)dtime
[2]	687	ENDDO
	688	DO k = 2, klev
	689	DO i = 1, knon
	690	local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
	691	local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
	692	ENDDO
	693	ENDDO
	694	c======================================================================
	695	c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
	696	c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
	697	DO k = 2, klev
	698	DO i = 1, knon
	699	flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
	700	. * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
	701	flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
	702	. * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
	703	. / zx_pkh(i,k)
	704	ENDDO
	705	ENDDO
	706	c======================================================================
[77]	707	C Calcul tendances
[2]	708	DO k = 1, klev
	709	DO i = 1, knon
	710	d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
	711	d_q(i,k) = local_q(i,k) - q(i,k)
	712	ENDDO
	713	ENDDO
	714	c
	715	RETURN
	716	END
	717	SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
	718	e paprs,pplay,delp,
	719	s d_ven,flux_v)
	720	IMPLICIT none
	721	c======================================================================
	722	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	723	c Objet: diffusion vertical de la vitesse "ven"
	724	c======================================================================
	725	c Arguments:
	726	c dtime----input-R- intervalle du temps (en second)
	727	c u1lay----input-R- vent u de la premiere couche (m/s)
	728	c v1lay----input-R- vent v de la premiere couche (m/s)
	729	c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
	730	c le cisaillement du vent (dV/dz); la premiere
	731	c valeur indique la valeur de Cdrag (sans unite)
	732	c t--------input-R- temperature (K)
	733	c ven------input-R- vitesse horizontale (m/s)
	734	c paprs----input-R- pression a inter-couche (Pa)
	735	c pplay----input-R- pression au milieu de couche (Pa)
	736	c delp-----input-R- epaisseur de couche (Pa)
	737	c
	738	c
	739	c d_ven----output-R- le changement de "ven"
	740	c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
	741	c======================================================================
	742	#include "dimensions.h"
	743	#include "dimphy.h"
	744	INTEGER knon
	745	REAL dtime
	746	REAL u1lay(klon), v1lay(klon)
	747	REAL coef(klon,klev)
	748	REAL t(klon,klev), ven(klon,klev)
	749	REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
	750	REAL d_ven(klon,klev)
	751	REAL flux_v(klon,klev)
	752	c======================================================================
	753	#include "YOMCST.h"
	754	c======================================================================
	755	INTEGER i, k
	756	REAL zx_cv(klon,2:klev)
	757	REAL zx_dv(klon,2:klev)
	758	REAL zx_buf(klon)
	759	REAL zx_coef(klon,klev)
	760	REAL local_ven(klon,klev)
	761	REAL zx_alf1(klon), zx_alf2(klon)
	762	c======================================================================
	763	DO k = 1, klev
	764	DO i = 1, knon
	765	local_ven(i,k) = ven(i,k)
	766	ENDDO
	767	ENDDO
	768	c======================================================================
	769	DO i = 1, knon
	770	ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
	771	zx_alf1(i) = 1.0
	772	zx_alf2(i) = 1.0 - zx_alf1(i)
	773	zx_coef(i,1) = coef(i,1)
	774	. * (1.0+SQRT(u1lay(i)2+v1lay(i)2))
	775	. * pplay(i,1)/(RD*t(i,1))
	776	zx_coef(i,1) = zx_coef(i,1) * dtime*RG
	777	ENDDO
	778	c======================================================================
	779	DO k = 2, klev
	780	DO i = 1, knon
	781	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
	782	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
	783	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
	784	ENDDO
	785	ENDDO
	786	c======================================================================
	787	DO i = 1, knon
	788	zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
	789	zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
	790	zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
	791	. /zx_buf(i)
	792	ENDDO
	793	DO k = 3, klev
	794	DO i = 1, knon
	795	zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
	796	. + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
	797	zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
	798	. +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
	799	zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
	800	ENDDO
	801	ENDDO
	802	DO i = 1, knon
	803	local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
	804	. +zx_coef(i,klev)*zx_cv(i,klev) )
	805	. / ( delp(i,klev) + zx_coef(i,klev)
	806	. -zx_coef(i,klev)*zx_dv(i,klev) )
	807	ENDDO
	808	DO k = klev-1, 1, -1
	809	DO i = 1, knon
	810	local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
	811	ENDDO
	812	ENDDO
	813	c======================================================================
	814	c== flux_v est le flux de moment angulaire (positif vers bas)
	815	c== dont l'unite est: (kg m/s)/(m**2 s)
	816	DO i = 1, knon
	817	flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
	818	. (local_ven(i,1)zx_alf1(i)
	819	. +local_ven(i,2)*zx_alf2(i))
	820	ENDDO
	821	DO k = 2, klev
	822	DO i = 1, knon
	823	flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
	824	. * (local_ven(i,k)-local_ven(i,k-1))
	825	ENDDO
	826	ENDDO
	827	c
	828	DO k = 1, klev
	829	DO i = 1, knon
	830	d_ven(i,k) = local_ven(i,k) - ven(i,k)
	831	ENDDO
	832	ENDDO
	833	c
	834	RETURN
	835	END
	836	SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
	837	. ts, rugos,
	838	. u,v,t,q,
	839	. pcfm, pcfh)
	840	IMPLICIT none
	841	c======================================================================
	842	c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
	843	c (une version strictement identique a l'ancien modele)
	844	c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
	845	c coefficients d'echange turbulent dans l'atmosphere.
	846	c Arguments:
	847	c nsrf-----input-I- indicateur de la nature du sol
	848	c knon-----input-I- nombre de points a traiter
	849	c paprs----input-R- pression a chaque intercouche (en Pa)
	850	c pplay----input-R- pression au milieu de chaque couche (en Pa)
	851	c ts-------input-R- temperature du sol (en Kelvin)
	852	c rugos----input-R- longeur de rugosite (en m)
	853	c xlat-----input-R- latitude en degree
	854	c u--------input-R- vitesse u
	855	c v--------input-R- vitesse v
	856	c t--------input-R- temperature (K)
	857	c q--------input-R- vapeur d'eau (kg/kg)
	858	c
	859	c itop-----output-I- numero de couche du sommet de la couche limite
	860	c pcfm-----output-R- coefficients a calculer (vitesse)
	861	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	862	c======================================================================
	863	#include "dimensions.h"
	864	#include "dimphy.h"
	865	#include "YOMCST.h"
	866	#include "indicesol.h"
	867	c
	868	c Arguments:
	869	c
	870	INTEGER knon, nsrf
	871	REAL ts(klon)
	872	REAL paprs(klon,klev+1), pplay(klon,klev)
	873	REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
	874	REAL rugos(klon)
	875	REAL xlat(klon)
	876	c
	877	REAL pcfm(klon,klev), pcfh(klon,klev)
	878	INTEGER itop(klon)
	879	c
	880	c Quelques constantes et options:
	881	c
	882	REAL cepdu2, ckap, cb, cc, cd, clam
	883	PARAMETER (cepdu2 =(0.1)**2)
	884	PARAMETER (ckap=0.35)
	885	PARAMETER (cb=5.0)
	886	PARAMETER (cc=5.0)
	887	PARAMETER (cd=5.0)
	888	PARAMETER (clam=160.0)
	889	REAL ratqs ! largeur de distribution de vapeur d'eau
[39]	890	PARAMETER (ratqs=0.05)
[2]	891	LOGICAL richum ! utilise le nombre de Richardson humide
	892	PARAMETER (richum=.TRUE.)
	893	REAL ric ! nombre de Richardson critique
	894	PARAMETER(ric=0.4)
	895	REAL prandtl
	896	PARAMETER (prandtl=0.4)
	897	REAL kstable ! diffusion minimale (situation stable)
[39]	898	PARAMETER (kstable=1.0e-10)
[2]	899	REAL mixlen ! constante controlant longueur de melange
	900	PARAMETER (mixlen=35.0)
	901	INTEGER isommet ! le sommet de la couche limite
	902	PARAMETER (isommet=klev)
	903	LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
	904	PARAMETER (tvirtu=.TRUE.)
	905	LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
	906	PARAMETER (opt_ec=.FALSE.)
	907	LOGICAL contreg ! utiliser le contre-gradient dans Ri
	908	PARAMETER (contreg=.TRUE.)
	909	c
	910	c Variables locales:
	911	c
	912	INTEGER i, k
	913	REAL zgeop(klon,klev)
	914	REAL zmgeom(klon)
	915	REAL zri(klon)
	916	REAL zl2(klon)
	917	REAL zcfm1(klon), zcfm2(klon)
	918	REAL zcfh1(klon), zcfh2(klon)
	919	REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
	920	REAL zscf, zucf, zcr
	921	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
	922	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
	923	REAL t_coup
	924	PARAMETER (t_coup=273.15)
	925	c
	926	c contre-gradient pour la chaleur sensible: Kelvin/metre
	927	REAL gamt(2:klev)
	928	c
	929	LOGICAL appel1er
	930	SAVE appel1er
	931	c
	932	c Fonctions thermodynamiques et fonctions d'instabilite
	933	REAL fsta, fins, x
	934	LOGICAL zxli ! utiliser un jeu de fonctions simples
	935	PARAMETER (zxli=.FALSE.)
	936	c
	937	#include "YOETHF.h"
	938	#include "FCTTRE.h"
	939	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
	940	fins(x) = SQRT(1.0-18.0*x)
	941	c
	942	DATA appel1er /.TRUE./
	943	c
	944	IF (appel1er) THEN
	945	PRINT*, 'coefkz, opt_ec:', opt_ec
	946	PRINT*, 'coefkz, richum:', richum
	947	IF (richum) PRINT*, 'coefkz, ratqs:', ratqs
	948	PRINT*, 'coefkz, isommet:', isommet
	949	PRINT*, 'coefkz, tvirtu:', tvirtu
	950	appel1er = .FALSE.
	951	ENDIF
	952	c
	953	c Initialiser les sorties
	954	c
	955	DO k = 1, klev
	956	DO i = 1, knon
	957	pcfm(i,k) = 0.0
	958	pcfh(i,k) = 0.0
	959	ENDDO
	960	ENDDO
	961	DO i = 1, knon
	962	itop(i) = 0
	963	ENDDO
	964	c
	965	c Prescrire la valeur de contre-gradient
	966	c
[39]	967	IF (.NOT.contreg) THEN
[2]	968	DO k = 2, klev
	969	gamt(k) = 0.0
	970	ENDDO
	971	ELSE
	972	DO k = 3, klev
	973	gamt(k) = -1.0E-03
	974	ENDDO
	975	gamt(2) = -2.5E-03
	976	ENDIF
	977	c
	978	c Calculer les geopotentiels de chaque couche
	979	c
	980	DO i = 1, knon
	981	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
	982	. * (paprs(i,1)-pplay(i,1))
	983	ENDDO
	984	DO k = 2, klev
	985	DO i = 1, knon
	986	zgeop(i,k) = zgeop(i,k-1)
	987	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
	988	. * (pplay(i,k-1)-pplay(i,k))
	989	ENDDO
	990	ENDDO
	991	c
	992	c Calculer le frottement au sol (Cdrag)
	993	c
	994	DO i = 1, knon
	995	zdu2=max(cepdu2,u(i,1)2+v(i,1)2)
	996	zdphi=zgeop(i,1)
	997	ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
	998	ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
	999	. (1.+RETVq(i,1))
	1000	zri(i)=zgeop(i,1)(ztvd-ztsolv)/(zdu2ztvd)
	1001	zcdn = (ckap/log(1.+zgeop(i,1)/(RGrugos(i))))*2
	1002	IF (zri(i) .ge. 0.) THEN ! situation stable
	1003	IF (.NOT.zxli) THEN
	1004	zscf=SQRT(1.+cd*ABS(zri(i)))
	1005	zcfm1(i) = zcdn/(1.+2.0cbzri(i)/zscf)
	1006	zcfh1(i) = zcdn/(1.+3.0cbzri(i)*zscf)
	1007	pcfm(i,1) = zcfm1(i)
	1008	pcfh(i,1) = zcfh1(i)
	1009	ELSE
	1010	pcfm(i,1) = zcdn* fsta(zri(i))
	1011	pcfh(i,1) = zcdn* fsta(zri(i))
	1012	ENDIF
	1013	ELSE ! situation instable
	1014	IF (.NOT.zxli) THEN
	1015	zucf=1./(1.+3.0cbcczcdnSQRT(ABS(zri(i))
	1016	. (1.0+zgeop(i,1)/(RGrugos(i)))))
	1017	zcfm2(i) = zcdn(1.-2.0cbzri(i)zucf)
	1018	zcfh2(i) = zcdn(1.-3.0cbzri(i)zucf)
	1019	pcfm(i,1) = zcfm2(i)
	1020	pcfh(i,1) = zcfh2(i)
	1021	ELSE
	1022	pcfm(i,1) = zcdn* fins(zri(i))
	1023	pcfh(i,1) = zcdn* fins(zri(i))
	1024	ENDIF
	1025	zcr = (0.0016/(zcdnSQRT(zdu2)))ABS(ztvd-ztsolv)**(1./3.)
	1026	IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn(1.0+zcr1.25)*(1./1.25)
	1027	ENDIF
	1028	ENDDO
	1029	c
	1030	c Calculer les coefficients turbulents dans l'atmosphere
	1031	c
	1032	DO i = 1, knon
	1033	itop(i) = isommet
	1034	ENDDO
	1035
	1036	DO k = 2, isommet
	1037	DO i = 1, knon
	1038	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
	1039	. +(v(i,k)-v(i,k-1))**2)
	1040	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
	1041	zdphi =zmgeom(i) / 2.0
	1042	zt = (t(i,k)+t(i,k-1)) * 0.5
	1043	zq = (q(i,k)+q(i,k-1)) * 0.5
	1044	c
	1045	c calculer Qs et dQs/dT:
	1046	c
	1047	IF (thermcep) THEN
	1048	zdelta = MAX(0.,SIGN(1.,RTT-zt))
	1049	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta)
	1050	. + R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
	1051	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
	1052	zqs = MIN(0.5,zqs)
	1053	zcor = 1./(1.-RETV*zqs)
	1054	zqs = zqs*zcor
	1055	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
	1056	ELSE
	1057	IF (zt .LT. t_coup) THEN
	1058	zqs = qsats(zt) / pplay(i,k)
	1059	zdqs = dqsats(zt,zqs)
	1060	ELSE
	1061	zqs = qsatl(zt) / pplay(i,k)
	1062	zdqs = dqsatl(zt,zqs)
	1063	ENDIF
	1064	ENDIF
	1065	c
	1066	c calculer la fraction nuageuse (processus humide):
	1067	c
	1068	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
	1069	zfr = MAX(0.0,MIN(1.0,zfr))
	1070	IF (.NOT.richum) zfr = 0.0
	1071	c
	1072	c calculer le nombre de Richardson:
	1073	c
	1074	IF (tvirtu) THEN
	1075	ztvd =( t(i,k)
	1076	. + zdphi/RCPD/(1.+RVTMP2*zq)
	1077	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
	1078	. )(1.+RETVq(i,k))
	1079	ztvu =( t(i,k-1)
	1080	. - zdphi/RCPD/(1.+RVTMP2*zq)
	1081	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
	1082	. )(1.+RETVq(i,k-1))
	1083	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
	1084	zri(i) = zri(i)
	1085	. + zmgeom(i)zmgeom(i)/RGgamt(k)
	1086	. (paprs(i,k)/101325.0)*RKAPPA
	1087	. /(zdu20.5(ztvd+ztvu))
	1088	c
	1089	ELSE ! calcul de Ridchardson compatible LMD5
	1090	c
	1091	zri(i) =(RCPD*(t(i,k)-t(i,k-1))
	1092	. -RD0.5(t(i,k)+t(i,k-1))/paprs(i,k)
	1093	. *(pplay(i,k)-pplay(i,k-1))
	1094	. )zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
	1095	zri(i) = zri(i) +
	1096	. zmgeom(i)zmgeom(i)gamt(k)/RG
	1097	cSB . /(paprs(i,k)/101325.0)**RKAPPA
	1098	. (paprs(i,k)/101325.0)*RKAPPA
	1099	. /(zdu20.5(t(i,k-1)+t(i,k)))
	1100	ENDIF
	1101	c
	1102	c finalement, les coefficients d'echange sont obtenus:
	1103	c
	1104	zcdn=SQRT(zdu2) / zmgeom(i) * RG
	1105	c
	1106	IF (opt_ec) THEN
	1107	z2geomf=zgeop(i,k-1)+zgeop(i,k)
	1108	zalm2=(0.5ckap/RGz2geomf
	1109	. /(1.+0.5ckap/rg/clamz2geomf))**2
	1110	zalh2=(0.5ckap/rgz2geomf
	1111	. /(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
	1112	IF (zri(i).LT.0.0) THEN ! situation instable
	1113	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3
	1114	. / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
	1115	zscf = SQRT(-zri(i)*zscf)
	1116	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
	1117	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
	1118	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
	1119	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
	1120	ELSE ! situation stable
	1121	zscf=SQRT(1.+cd*zri(i))
	1122	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
	1123	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
	1124	ENDIF
	1125	ELSE
	1126	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
	1127	. /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
	1128	pcfm(i,k)=sqrt(max(zcdnzcdn(ric-zri(i))/ric, kstable))
	1129	pcfm(i,k)= zl2(i)* pcfm(i,k)
	1130	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	1131	ENDIF
	1132	ENDDO
	1133	ENDDO
	1134	c
	1135	c Au-dela du sommet, pas de diffusion turbulente:
	1136	c
	1137	DO i = 1, knon
	1138	IF (itop(i)+1 .LE. klev) THEN
	1139	DO k = itop(i)+1, klev
	1140	pcfh(i,k) = 0.0
	1141	pcfm(i,k) = 0.0
	1142	ENDDO
	1143	ENDIF
	1144	ENDDO
	1145	c
	1146	RETURN
	1147	END
[39]	1148	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
	1149	. pcfm, pcfh)
	1150	IMPLICIT none
	1151	c======================================================================
	1152	c J'introduit un peu de diffusion sauf dans les endroits
	1153	c ou une forte inversion est presente
	1154	c On peut dire qu'il represente la convection peu profonde
	1155	c
	1156	c Arguments:
	1157	c nsrf-----input-I- indicateur de la nature du sol
	1158	c knon-----input-I- nombre de points a traiter
	1159	c paprs----input-R- pression a chaque intercouche (en Pa)
	1160	c pplay----input-R- pression au milieu de chaque couche (en Pa)
	1161	c t--------input-R- temperature (K)
	1162	c
	1163	c pcfm-----output-R- coefficients a calculer (vitesse)
	1164	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	1165	c======================================================================
	1166	#include "dimensions.h"
	1167	#include "dimphy.h"
	1168	#include "YOMCST.h"
	1169	#include "indicesol.h"
	1170	c
	1171	c Arguments:
	1172	c
	1173	INTEGER knon, nsrf
[149]	1174	REAL paprs(klon,klev+1), pplay(klon,klev)
	1175	REAL t(klon,klev)
[39]	1176	c
[149]	1177	REAL pcfm(klon,klev), pcfh(klon,klev)
[39]	1178	c
	1179	c Quelques constantes et options:
	1180	c
	1181	REAL prandtl
	1182	PARAMETER (prandtl=0.4)
	1183	REAL kstable
	1184	PARAMETER (kstable=0.002)
	1185	ccc PARAMETER (kstable=0.001)
	1186	REAL mixlen ! constante controlant longueur de melange
	1187	PARAMETER (mixlen=35.0)
	1188	REAL seuil ! au-dela l'inversion est consideree trop faible
	1189	PARAMETER (seuil=-0.02)
	1190	ccc PARAMETER (seuil=-0.04)
	1191	ccc PARAMETER (seuil=-0.06)
	1192	ccc PARAMETER (seuil=-0.09)
	1193	c
	1194	c Variables locales:
	1195	c
[145]	1196	INTEGER i, k, invb(knon)
	1197	REAL zl2(knon)
	1198	REAL zdthmin(knon), zdthdp
[39]	1199	c
	1200	c Initialiser les sorties
	1201	c
	1202	DO k = 1, klev
	1203	DO i = 1, knon
	1204	pcfm(i,k) = 0.0
	1205	pcfh(i,k) = 0.0
	1206	ENDDO
	1207	ENDDO
	1208	c
	1209	c Chercher la zone d'inversion forte
	1210	c
[145]	1211	DO i = 1, knon
[39]	1212	invb(i) = klev
	1213	zdthmin(i)=0.0
	1214	ENDDO
	1215	DO k = 2, klev/2-1
[145]	1216	DO i = 1, knon
[39]	1217	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
	1218	. - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
	1219	zdthdp = zdthdp * 100.0
	1220	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
	1221	. zdthdp.LT.zdthmin(i) ) THEN
	1222	zdthmin(i) = zdthdp
	1223	invb(i) = k
	1224	ENDIF
	1225	ENDDO
	1226	ENDDO
	1227	c
	1228	c Introduire une diffusion:
	1229	c
	1230	DO k = 2, klev
	1231	DO i = 1, knon
	1232	IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
	1233	. (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
	1234	. (zdthmin(i).GT.seuil) )THEN ! si l'inversion est trop faible
	1235	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
	1236	. /(paprs(i,2)-paprs(i,klev+1)) ))**2
	1237	pcfm(i,k)= zl2(i)* kstable
	1238	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	1239	ENDIF
	1240	ENDDO
	1241	ENDDO
	1242	c
	1243	RETURN
	1244	END
[2]	1245	SUBROUTINE calbeta(dtime,indice,snow,qsol,
	1246	. vbeta,vcal,vdif)
	1247	IMPLICIT none
	1248	c======================================================================
	1249	c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
	1250	c date: 19940414
	1251	c======================================================================
	1252	c
	1253	c Calculer quelques parametres pour appliquer la couche limite
	1254	c ------------------------------------------------------------
	1255	#include "dimensions.h"
	1256	#include "dimphy.h"
	1257	#include "YOMCST.h"
	1258	#include "indicesol.h"
	1259	REAL tau_gl ! temps de relaxation pour la glace de mer
[39]	1260	ccc PARAMETER (tau_gl=86400.0*30.0)
	1261	PARAMETER (tau_gl=86400.0*5.0)
[2]	1262	REAL mx_eau_sol
	1263	PARAMETER (mx_eau_sol=150.0)
	1264	c
	1265	REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
	1266	PARAMETER (calsol=1.0/(2.5578E+06*0.15))
	1267	PARAMETER (calsno=1.0/(2.3867E+06*0.15))
	1268	PARAMETER (calice=1.0/(5.1444E+06*0.15))
	1269	C
	1270	INTEGER i
	1271	c
	1272	REAL dtime
	1273	REAL snow(klon,nbsrf), qsol(klon,nbsrf)
	1274	INTEGER indice
	1275	C
	1276	REAL vbeta(klon)
	1277	REAL vcal(klon)
	1278	REAL vdif(klon)
	1279	C
	1280	IF (indice.EQ.is_oce) THEN
	1281	DO i = 1, klon
	1282	vcal(i) = 0.0
	1283	vbeta(i) = 1.0
	1284	vdif(i) = 0.0
	1285	ENDDO
	1286	ENDIF
	1287	c
	1288	IF (indice.EQ.is_sic) THEN
	1289	DO i = 1, klon
	1290	vcal(i) = calice
	1291	IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
	1292	vbeta(i) = 1.0
[39]	1293	vdif(i) = 1.0/tau_gl
	1294	ccc vdif(i) = calice/tau_gl ! c'etait une erreur
[2]	1295	ENDDO
	1296	ENDIF
	1297	c
	1298	IF (indice.EQ.is_ter) THEN
	1299	DO i = 1, klon
	1300	vcal(i) = calsol
	1301	IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
	1302	vbeta(i) = MIN(2.0*qsol(i,is_ter)/mx_eau_sol, 1.0)
	1303	vdif(i) = 0.0
	1304	ENDDO
	1305	ENDIF
	1306	c
	1307	IF (indice.EQ.is_lic) THEN
	1308	DO i = 1, klon
	1309	vcal(i) = calice
	1310	IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
	1311	vbeta(i) = 1.0
	1312	vdif(i) = 0.0
	1313	ENDDO
	1314	ENDIF
	1315	c
	1316	RETURN
	1317	END
	1318	C======================================================================
	1319	SUBROUTINE nonlocal(knon, paprs, pplay,
	1320	. tsol,beta,u,v,t,q,
	1321	. cd_h, cd_m, pcfh, pcfm, cgh, cgq)
	1322	IMPLICIT none
	1323	c======================================================================
	1324	c Laurent Li (LMD/CNRS), le 30 septembre 1998
	1325	c Couche limite non-locale. Adaptation du code du CCM3.
	1326	c Code non teste, donc a ne pas utiliser.
	1327	c======================================================================
	1328	c Nonlocal scheme that determines eddy diffusivities based on a
	1329	c diagnosed boundary layer height and a turbulent velocity scale.
	1330	c Also countergradient effects for heat and moisture are included.
	1331	c
	1332	c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
	1333	c Local versus nonlocal boundary-layer diffusion in a global climate
	1334	c model. J. of Climate, vol. 6, 1825-1842.
	1335	c======================================================================
	1336	#include "dimensions.h"
	1337	#include "dimphy.h"
	1338	#include "YOMCST.h"
	1339	c
	1340	c Arguments:
	1341	c
	1342	INTEGER knon ! nombre de points a calculer
	1343	REAL tsol(klon) ! temperature du sol (K)
	1344	REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
	1345	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
	1346	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
	1347	REAL u(klon,klev) ! vitesse U (m/s)
	1348	REAL v(klon,klev) ! vitesse V (m/s)
	1349	REAL t(klon,klev) ! temperature (K)
	1350	REAL q(klon,klev) ! vapeur d'eau (kg/kg)
	1351	REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
	1352	REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
	1353	c
	1354	INTEGER isommet
	1355	PARAMETER (isommet=klev)
	1356	REAL vk
	1357	PARAMETER (vk=0.35)
	1358	REAL ricr
	1359	PARAMETER (ricr=0.4)
	1360	REAL fak
	1361	PARAMETER (fak=8.5)
	1362	REAL fakn
	1363	PARAMETER (fakn=7.2)
	1364	REAL onet
	1365	PARAMETER (onet=1.0/3.0)
	1366	REAL t_coup
	1367	PARAMETER(t_coup=273.15)
	1368	REAL zkmin
	1369	PARAMETER (zkmin=0.01)
	1370	REAL betam
	1371	PARAMETER (betam=15.0)
	1372	REAL betah
	1373	PARAMETER (betah=15.0)
	1374	REAL betas
	1375	PARAMETER (betas=5.0)
	1376	REAL sffrac
	1377	PARAMETER (sffrac=0.1)
	1378	REAL binm
	1379	PARAMETER (binm=betam*sffrac)
	1380	REAL binh
	1381	PARAMETER (binh=betah*sffrac)
	1382	REAL ccon
	1383	PARAMETER (ccon=faksffracvk)
	1384	c
	1385	REAL z(klon,klev)
	1386	REAL pcfm(klon,klev), pcfh(klon,klev)
	1387	c
	1388	INTEGER i, k
	1389	REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
	1390	REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
	1391	REAL khfs(klon) ! surface kinematic heat flux [mK/s]
	1392	REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
	1393	REAL heatv(klon) ! surface virtual heat flux
	1394	REAL ustar(klon)
	1395	REAL rino(klon,klev) ! bulk Richardon no. from level to ref lev
	1396	LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
	1397	LOGICAL stblev(klon) ! stable pbl with levels within pbl
	1398	LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
	1399	LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
	1400	LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
	1401	LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
	1402	REAL pblh(klon)
	1403	REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
	1404	REAL cgq(klon,2:klev) ! counter-gradient term for constituents
	1405	REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
	1406	REAL obklen(klon)
	1407	REAL ztvd, ztvu, zdu2
	1408	REAL therm(klon) ! thermal virtual temperature excess
	1409	REAL phiminv(klon) ! inverse phi function for momentum
	1410	REAL phihinv(klon) ! inverse phi function for heat
	1411	REAL wm(klon) ! turbulent velocity scale for momentum
	1412	REAL fak1(klon) ! kustarpblh
	1413	REAL fak2(klon) ! kwmpblh
	1414	REAL fak3(klon) ! fakn*wstr/wm
	1415	REAL pblk(klon) ! level eddy diffusivity for momentum
	1416	REAL pr(klon) ! Prandtl number for eddy diffusivities
	1417	REAL zl(klon) ! zmzp / Obukhov length
	1418	REAL zh(klon) ! zmzp / pblh
	1419	REAL zzh(klon) ! (1-(zmzp/pblh))**2
	1420	REAL wstr(klon) ! w*, convective velocity scale
	1421	REAL zm(klon) ! current level height
	1422	REAL zp(klon) ! current level height + one level up
	1423	REAL zcor, zdelta, zcvm5, zxqs
	1424	REAL fac, pblmin, zmzp, term
	1425	c
	1426	#include "YOETHF.h"
	1427	#include "FCTTRE.h"
	1428	c
	1429	c Initialisation
	1430	c
	1431	DO i = 1, klon
	1432	pcfh(i,1) = cd_h(i)
	1433	pcfm(i,1) = cd_m(i)
	1434	ENDDO
	1435	DO k = 2, klev
	1436	DO i = 1, klon
	1437	pcfh(i,k) = zkmin
	1438	pcfm(i,k) = zkmin
	1439	cgs(i,k) = 0.0
	1440	cgh(i,k) = 0.0
	1441	cgq(i,k) = 0.0
	1442	ENDDO
	1443	ENDDO
	1444	c
	1445	c Calculer les hauteurs de chaque couche
	1446	c
	1447	DO i = 1, knon
	1448	z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
	1449	. * (paprs(i,1)-pplay(i,1)) / RG
	1450	ENDDO
	1451	DO k = 2, klev
	1452	DO i = 1, knon
	1453	z(i,k) = z(i,k-1)
	1454	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
	1455	. * (pplay(i,k-1)-pplay(i,k)) / RG
	1456	ENDDO
	1457	ENDDO
	1458	c
	1459	DO i = 1, knon
	1460	IF (thermcep) THEN
	1461	zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
	1462	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	1463	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
	1464	zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
	1465	zxqs=MIN(0.5,zxqs)
	1466	zcor=1./(1.-retv*zxqs)
	1467	zxqs=zxqs*zcor
	1468	ELSE
	1469	IF (tsol(i).LT.t_coup) THEN
	1470	zxqs = qsats(tsol(i)) / paprs(i,1)
	1471	ELSE
	1472	zxqs = qsatl(tsol(i)) / paprs(i,1)
	1473	ENDIF
	1474	ENDIF
	1475	zx_alf1 = 1.0
	1476	zx_alf2 = 1.0 - zx_alf1
	1477	zxt = (t(i,1)+z(i,1)RG/RCPD/(1.+RVTMP2q(i,1)))
	1478	. (1.+RETVq(i,1))*zx_alf1
	1479	. + (t(i,2)+z(i,2)RG/RCPD/(1.+RVTMP2q(i,2)))
	1480	. (1.+RETVq(i,2))*zx_alf2
	1481	zxu = u(i,1)zx_alf1+u(i,2)zx_alf2
	1482	zxv = v(i,1)zx_alf1+v(i,2)zx_alf2
	1483	zxq = q(i,1)zx_alf1+q(i,2)zx_alf2
	1484	zxmod = 1.0+SQRT(zxu2+zxv2)
	1485	khfs(i) = (tsol(i)(1.+RETVq(i,1))-zxt) zxmodcd_h(i)
	1486	kqfs(i) = (zxqs-zxq) zxmodcd_h(i) * beta(i)
	1487	heatv(i) = khfs(i) + 0.61zxtkqfs(i)
	1488	taux = zxu zxmodcd_m(i)
	1489	tauy = zxv zxmodcd_m(i)
	1490	ustar(i) = SQRT(taux2+tauy2)
	1491	ustar(i) = MAX(SQRT(ustar(i)),0.01)
	1492	ENDDO
	1493	c
	1494	DO i = 1, knon
	1495	rino(i,1) = 0.0
	1496	check(i) = .TRUE.
	1497	pblh(i) = z(i,1)
	1498	obklen(i) = -t(i,1)ustar(i)3/(RGvk*heatv(i))
	1499	ENDDO
	1500
	1501	C
	1502	C PBL height calculation:
	1503	C Search for level of pbl. Scan upward until the Richardson number between
	1504	C the first level and the current level exceeds the "critical" value.
	1505	C
	1506	fac = 100.0
	1507	DO k = 1, isommet
	1508	DO i = 1, knon
	1509	IF (check(i)) THEN
	1510	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
	1511	zdu2 = max(zdu2,1.0e-20)
	1512	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
	1513	. (1.+RETVq(i,k))
	1514	ztvu =(t(i,1)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
	1515	. (1.+RETVq(i,1))
	1516	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
	1517	. /(zdu20.5(ztvd+ztvu))
	1518	IF (rino(i,k).GE.ricr) THEN
	1519	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
	1520	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
	1521	check(i) = .FALSE.
	1522	ENDIF
	1523	ENDIF
	1524	ENDDO
	1525	ENDDO
	1526
	1527	C
	1528	C Set pbl height to maximum value where computation exceeds number of
	1529	C layers allowed
	1530	C
	1531	DO i = 1, knon
	1532	if (check(i)) pblh(i) = z(i,isommet)
	1533	ENDDO
	1534	C
	1535	C Improve estimate of pbl height for the unstable points.
	1536	C Find unstable points (sensible heat flux is upward):
	1537	C
	1538	DO i = 1, knon
	1539	IF (heatv(i) .GT. 0.) THEN
	1540	unstbl(i) = .TRUE.
	1541	check(i) = .TRUE.
	1542	ELSE
	1543	unstbl(i) = .FALSE.
	1544	check(i) = .FALSE.
	1545	ENDIF
	1546	ENDDO
	1547	C
	1548	C For the unstable case, compute velocity scale and the
	1549	C convective temperature excess:
	1550	C
	1551	DO i = 1, knon
	1552	IF (check(i)) THEN
	1553	phiminv(i) = (1.-binmpblh(i)/obklen(i))*onet
	1554	wm(i)= ustar(i)*phiminv(i)
	1555	therm(i) = heatv(i)*fak/wm(i)
	1556	rino(i,1) = 0.0
	1557	ENDIF
	1558	ENDDO
	1559	C
	1560	C Improve pblh estimate for unstable conditions using the
	1561	C convective temperature excess:
	1562	C
	1563	DO k = 1, isommet
	1564	DO i = 1, knon
	1565	IF (check(i)) THEN
	1566	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
	1567	zdu2 = max(zdu2,1.0e-20)
	1568	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
	1569	. (1.+RETVq(i,k))
	1570	ztvu =(t(i,1)+therm(i)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
	1571	. (1.+RETVq(i,1))
	1572	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
	1573	. /(zdu20.5(ztvd+ztvu))
	1574	IF (rino(i,k).GE.ricr) THEN
	1575	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
	1576	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
	1577	check(i) = .FALSE.
	1578	ENDIF
	1579	ENDIF
	1580	ENDDO
	1581	ENDDO
	1582	C
	1583	C Set pbl height to maximum value where computation exceeds number of
	1584	C layers allowed
	1585	C
	1586	DO i = 1, knon
	1587	if (check(i)) pblh(i) = z(i,isommet)
	1588	ENDDO
	1589	C
	1590	C Points for which pblh exceeds number of pbl layers allowed;
	1591	C set to maximum
	1592	C
	1593	DO i = 1, knon
	1594	IF (check(i)) pblh(i) = z(i,isommet)
	1595	ENDDO
	1596	C
	1597	C PBL height must be greater than some minimum mechanical mixing depth
	1598	C Several investigators have proposed minimum mechanical mixing depth
	1599	C relationships as a function of the local friction velocity, u*. We
	1600	C make use of a linear relationship of the form h = c u* where c=700.
	1601	C The scaling arguments that give rise to this relationship most often
	1602	C represent the coefficient c as some constant over the local coriolis
	1603	C parameter. Here we make use of the experimental results of Koracin
	1604	C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
	1605	C where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
	1606	C latitude value for f so that c = 0.07/f = 700.
	1607	C
	1608	DO i = 1, knon
	1609	pblmin = 700.0*ustar(i)
	1610	pblh(i) = MAX(pblh(i),pblmin)
	1611	ENDDO
	1612	C
	1613	C pblh is now available; do preparation for diffusivity calculation:
	1614	C
	1615	DO i = 1, knon
	1616	pblk(i) = 0.0
	1617	fak1(i) = ustar(i)pblh(i)vk
	1618	C
	1619	C Do additional preparation for unstable cases only, set temperature
	1620	C and moisture perturbations depending on stability.
	1621	C
	1622	IF (unstbl(i)) THEN
	1623	zxt=(t(i,1)-z(i,1)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
	1624	. (1.+RETVq(i,1))
	1625	phiminv(i) = (1. - binmpblh(i)/obklen(i))*onet
	1626	phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
	1627	wm(i) = ustar(i)*phiminv(i)
	1628	fak2(i) = wm(i)pblh(i)vk
	1629	wstr(i) = (heatv(i)RGpblh(i)/zxt)**onet
	1630	fak3(i) = fakn*wstr(i)/wm(i)
	1631	ENDIF
	1632	ENDDO
	1633
	1634	C Main level loop to compute the diffusivities and
	1635	C counter-gradient terms:
	1636	C
	1637	DO 1000 k = 2, isommet
	1638	C
	1639	C Find levels within boundary layer:
	1640	C
	1641	DO i = 1, knon
	1642	unslev(i) = .FALSE.
	1643	stblev(i) = .FALSE.
	1644	zm(i) = z(i,k-1)
	1645	zp(i) = z(i,k)
	1646	IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
	1647	IF (zm(i) .LT. pblh(i)) THEN
	1648	zmzp = 0.5*(zm(i) + zp(i))
	1649	zh(i) = zmzp/pblh(i)
	1650	zl(i) = zmzp/obklen(i)
	1651	zzh(i) = 0.
	1652	IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
	1653	C
	1654	C stblev for points zm < plbh and stable and neutral
	1655	C unslev for points zm < plbh and unstable
	1656	C
	1657	IF (unstbl(i)) THEN
	1658	unslev(i) = .TRUE.
	1659	ELSE
	1660	stblev(i) = .TRUE.
	1661	ENDIF
	1662	ENDIF
	1663	ENDDO
	1664	C
	1665	C Stable and neutral points; set diffusivities; counter-gradient
	1666	C terms zero for stable case:
	1667	C
	1668	DO i = 1, knon
	1669	IF (stblev(i)) THEN
	1670	IF (zl(i).LE.1.) THEN
	1671	pblk(i) = fak1(i)zh(i)zzh(i)/(1. + betas*zl(i))
	1672	ELSE
	1673	pblk(i) = fak1(i)zh(i)zzh(i)/(betas + zl(i))
	1674	ENDIF
	1675	pcfm(i,k) = pblk(i)
	1676	pcfh(i,k) = pcfm(i,k)
	1677	ENDIF
	1678	ENDDO
	1679	C
	1680	C unssrf, unstable within surface layer of pbl
	1681	C unsout, unstable within outer layer of pbl
	1682	C
	1683	DO i = 1, knon
	1684	unssrf(i) = .FALSE.
	1685	unsout(i) = .FALSE.
	1686	IF (unslev(i)) THEN
	1687	IF (zh(i).lt.sffrac) THEN
	1688	unssrf(i) = .TRUE.
	1689	ELSE
	1690	unsout(i) = .TRUE.
	1691	ENDIF
	1692	ENDIF
	1693	ENDDO
	1694	C
	1695	C Unstable for surface layer; counter-gradient terms zero
	1696	C
	1697	DO i = 1, knon
	1698	IF (unssrf(i)) THEN
	1699	term = (1. - betamzl(i))*onet
	1700	pblk(i) = fak1(i)zh(i)zzh(i)*term
	1701	pr(i) = term/sqrt(1. - betah*zl(i))
	1702	ENDIF
	1703	ENDDO
	1704	C
	1705	C Unstable for outer layer; counter-gradient terms non-zero:
	1706	C
	1707	DO i = 1, knon
	1708	IF (unsout(i)) THEN
	1709	pblk(i) = fak2(i)zh(i)zzh(i)
	1710	cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
	1711	cgh(i,k) = khfs(i)*cgs(i,k)
	1712	pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
	1713	cgq(i,k) = kqfs(i)*cgs(i,k)
	1714	ENDIF
	1715	ENDDO
	1716	C
	1717	C For all unstable layers, set diffusivities
	1718	C
	1719	DO i = 1, knon
	1720	IF (unslev(i)) THEN
	1721	pcfm(i,k) = pblk(i)
	1722	pcfh(i,k) = pblk(i)/pr(i)
	1723	ENDIF
	1724	ENDDO
	1725	1000 continue ! end of level loop
	1726
	1727	RETURN
	1728	END

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

Download in other formats:

Original Format