Context Navigation

clmain.F @ 828

Last change on this file since 828 was 546, checked in by lmdzadmin, 20 years ago
Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory) LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 69.0 KB

Rev	Line
[230]	1	c
	2	c $Header$
	3	c
[499]	4	SUBROUTINE clmain(dtime,itap,date0,pctsrf,pctsrf_new,
[112]	5	. t,q,u,v,
[475]	6	. jour, rmu0, co2_ppm,
[112]	7	. ok_veget, ocean, npas, nexca, ts,
[486]	8	. soil_model,cdmmax, cdhmax,
	9	. ksta, ksta_ter, ok_kzmin, ftsoil,qsol,
[438]	10	. paprs,pplay,radsol,snow,qsurf,evap,albe,alblw,
[280]	11	. fluxlat,
[177]	12	. rain_f, snow_f, solsw, sollw, sollwdown, fder,
[171]	13	. rlon, rlat, cufi, cvfi, rugos,
	14	. debut, lafin, agesno,rugoro,
[2]	15	. d_t,d_q,d_u,d_v,d_ts,
	16	. flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
[171]	17	. dflux_t,dflux_q,
[454]	18	cIM cf JLD . zcoefh,zu1,zv1, t2m, q2m, u10m, v10m)
	19	. zcoefh,zu1,zv1, t2m, q2m, u10m, v10m,
[522]	20	. fqcalving,ffonte, run_off_lic_0)
[2]	21	cAA . itr, tr, flux_surf, d_tr)
	22	cAA REM:
	23	cAA-----
	24	cAA Tout ce qui a trait au traceurs est dans phytrac maintenant
	25	cAA pour l'instant le calcul de la couche limite pour les traceurs
	26	cAA se fait avec cltrac et ne tient pas compte de la differentiation
	27	cAA des sous-fraction de sol.
	28	cAA REM bis :
	29	cAA----------
	30	cAA Pour pouvoir extraire les coefficient d'echanges et le vent
	31	cAA dans la premiere couche, 3 champs supplementaires ont ete crees
	32	cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
	33	cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
	34	cAA si les informations des subsurfaces doivent etre prises en compte
	35	cAA il faudra sortir ces memes champs en leur ajoutant une dimension,
	36	cAA c'est a dire nbsrf (nbre de subsurface).
[223]	37	USE ioipsl
	38	USE interface_surf
[2]	39	IMPLICIT none
	40	c======================================================================
	41	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	42	c Objet: interface de "couche limite" (diffusion verticale)
	43	c Arguments:
	44	c dtime----input-R- interval du temps (secondes)
[112]	45	c itap-----input-I- numero du pas de temps
[205]	46	c date0----input-R- jour initial
[2]	47	c t--------input-R- temperature (K)
	48	c q--------input-R- vapeur d'eau (kg/kg)
	49	c u--------input-R- vitesse u
	50	c v--------input-R- vitesse v
	51	c ts-------input-R- temperature du sol (en Kelvin)
	52	c paprs----input-R- pression a intercouche (Pa)
	53	c pplay----input-R- pression au milieu de couche (Pa)
	54	c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
[98]	55	c rlat-----input-R- latitude en degree
[2]	56	c rugos----input-R- longeur de rugosite (en m)
[177]	57	c cufi-----input-R- resolution des mailles en x (m)
	58	c cvfi-----input-R- resolution des mailles en y (m)
[2]	59	c
	60	c d_t------output-R- le changement pour "t"
	61	c d_q------output-R- le changement pour "q"
	62	c d_u------output-R- le changement pour "u"
	63	c d_v------output-R- le changement pour "v"
	64	c d_ts-----output-R- le changement pour "ts"
	65	c flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2)
	66	c (orientation positive vers le bas)
	67	c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
	68	c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
	69	c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
	70	c dflux_t derive du flux sensible
	71	c dflux_q derive du flux latent
[454]	72	c ffonte----Flux thermique utilise pour fondre la neige
	73	c fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
	74	c hauteur de neige, en kg/m2/s
[2]	75	cAA on rajoute en output yu1 et yv1 qui sont les vents dans
	76	cAA la premiere couche
	77	cAA ces 4 variables sont maintenant traites dans phytrac
	78	c itr--------input-I- nombre de traceurs
	79	c tr---------input-R- q. de traceurs
	80	c flux_surf--input-R- flux de traceurs a la surface
	81	c d_tr-------output-R tendance de traceurs
	82	c======================================================================
	83	#include "dimensions.h"
	84	#include "dimphy.h"
	85	#include "indicesol.h"
[177]	86	c$$$ PB ajout pour soil
	87	#include "dimsoil.h"
[2]	88	c
	89	REAL dtime
[205]	90	real date0
[112]	91	integer itap
[2]	92	REAL t(klon,klev), q(klon,klev)
	93	REAL u(klon,klev), v(klon,klev)
	94	REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
[171]	95	REAL rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
[2]	96	REAL d_t(klon, klev), d_q(klon, klev)
	97	REAL d_u(klon, klev), d_v(klon, klev)
[98]	98	REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
[2]	99	REAL dflux_t(klon), dflux_q(klon)
[454]	100	cIM cf JLD
	101	REAL y_fqcalving(klon), y_ffonte(klon)
	102	REAL fqcalving(klon,nbsrf), ffonte(klon,nbsrf)
[522]	103	REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
[454]	104
[98]	105	REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
[258]	106	REAL rugmer(klon), agesno(klon,nbsrf),rugoro(klon)
[2]	107	REAL cdragh(klon), cdragm(klon)
[109]	108	integer jour ! jour de l'annee en cours
	109	real rmu0(klon) ! cosinus de l'angle solaire zenithal
[475]	110	REAL co2_ppm ! taux CO2 atmosphere
[98]	111	LOGICAL debut, lafin, ok_veget
[112]	112	character*6 ocean
	113	integer npas, nexca
[2]	114	cAA INTEGER itr
	115	cAA REAL tr(klon,klev,nbtr)
	116	cAA REAL d_tr(klon,klev,nbtr)
	117	cAA REAL flux_surf(klon,nbtr)
	118	c
	119	REAL pctsrf(klon,nbsrf)
	120	REAL ts(klon,nbsrf)
	121	REAL d_ts(klon,nbsrf)
	122	REAL snow(klon,nbsrf)
[438]	123	REAL qsurf(klon,nbsrf)
[98]	124	REAL evap(klon,nbsrf)
	125	REAL albe(klon,nbsrf)
[280]	126	REAL alblw(klon,nbsrf)
[177]	127	c$$$ PB
	128	REAL fluxlat(klon,nbsrf)
	129	C
[98]	130	real rain_f(klon), snow_f(klon)
[101]	131	REAL fder(klon)
[433]	132	cIM cf. JLD REAL sollw(klon), solsw(klon), sollwdown(klon)
	133	REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)
[2]	134	REAL rugos(klon,nbsrf)
[105]	135	C la nouvelle repartition des surfaces sortie de l'interface
	136	REAL pctsrf_new(klon,nbsrf)
[2]	137	cAA
	138	REAL zcoefh(klon,klev)
	139	REAL zu1(klon)
	140	REAL zv1(klon)
	141	cAA
[177]	142	c$$$ PB ajout pour soil
	143	LOGICAL soil_model
[467]	144	cIM ajout seuils cdrm, cdrh
	145	REAL cdmmax, cdhmax
[486]	146	cIM: 261103
	147	REAL ksta, ksta_ter
	148	LOGICAL ok_kzmin
	149	cIM: 261103
[177]	150	REAL ftsoil(klon,nsoilmx,nbsrf)
	151	REAL ytsoil(klon,nsoilmx)
[438]	152	REAL qsol(klon)
[2]	153	c======================================================================
	154	EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
	155	c======================================================================
[109]	156	REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
[276]	157	REAL yalb(klon)
[280]	158	REAL yalblw(klon)
[2]	159	REAL yu1(klon), yv1(klon)
[438]	160	real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
[98]	161	real yrain_f(klon), ysnow_f(klon)
[177]	162	real ysollw(klon), ysolsw(klon), ysollwdown(klon)
[101]	163	real yfder(klon), ytaux(klon), ytauy(klon)
[171]	164	REAL yrugm(klon), yrads(klon),yrugoro(klon)
[177]	165	c$$$ PB
	166	REAL yfluxlat(klon)
	167	C
[2]	168	REAL y_d_ts(klon)
	169	REAL y_d_t(klon, klev), y_d_q(klon, klev)
	170	REAL y_d_u(klon, klev), y_d_v(klon, klev)
	171	REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
	172	REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
	173	REAL y_dflux_t(klon), y_dflux_q(klon)
	174	REAL ycoefh(klon,klev), ycoefm(klon,klev)
	175	REAL yu(klon,klev), yv(klon,klev)
	176	REAL yt(klon,klev), yq(klon,klev)
	177	REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
	178	cAA REAL ytr(klon,klev,nbtr)
	179	cAA REAL y_d_tr(klon,klev,nbtr)
	180	cAA REAL yflxsrf(klon,nbtr)
	181	c
[39]	182	LOGICAL contreg
	183	PARAMETER (contreg=.TRUE.)
	184	c
[2]	185	LOGICAL ok_nonloc
	186	PARAMETER (ok_nonloc=.FALSE.)
	187	REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
	188	c
	189	#include "YOMCST.h"
[407]	190	#include "YOETHF.h"
	191	#include "FCTTRE.h"
[2]	192	REAL u1lay(klon), v1lay(klon)
	193	REAL delp(klon,klev)
	194	INTEGER i, k, nsrf
	195	cAA INTEGER it
	196	INTEGER ni(klon), knon, j
[105]	197	c Introduction d'une variable "pourcentage potentiel" pour tenir compte
	198	c des eventuelles apparitions et/ou disparitions de la glace de mer
	199	REAL pctsrf_pot(klon,nbsrf)
	200
[2]	201	c======================================================================
	202	REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
	203	c======================================================================
[223]	204	c
	205	c maf pour sorties IOISPL en cas de debugagage
	206	c
	207	CHARACTER*80 cldebug
	208	SAVE cldebug
	209	CHARACTER*8 cl_surf(nbsrf)
	210	SAVE cl_surf
	211	INTEGER nhoridbg, nidbg
	212	SAVE nhoridbg, nidbg
	213	INTEGER ndexbg(iim*(jjm+1))
	214	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
	215	REAL tabindx(klon)
	216	REAL debugtab(iim,jjm+1)
	217	LOGICAL first_appel
	218	SAVE first_appel
[295]	219	DATA first_appel/.false./
[223]	220	LOGICAL debugindex
	221	SAVE debugindex
[295]	222	DATA debugindex/.false./
[353]	223	integer idayref
[223]	224	#include "temps.h"
[407]	225	REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
	226	REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
	227	c
	228	REAL yt2m(klon), yq2m(klon), yu10m(klon)
	229	c
	230	REAL uzon(klon), vmer(klon)
	231	REAL tair1(klon), qair1(klon), tairsol(klon)
	232	REAL psfce(klon), patm(klon)
	233	c
	234	REAL qairsol(klon), zgeo1(klon)
	235	REAL rugo1(klon)
	236	c
	237	LOGICAL zxli ! utiliser un jeu de fonctions simples
	238	PARAMETER (zxli=.FALSE.)
	239	c
	240	REAL zt, zqs, zdelta, zcor
	241	REAL t_coup
	242	PARAMETER(t_coup=273.15)
	243	C
[438]	244	character (len = 20) :: modname = 'clmain'
	245	LOGICAL check
[458]	246	PARAMETER (check=.false.)
[438]	247	C
	248	if (check) THEN
	249	write(,) modname,' klon=',klon
	250	call flush(6)
	251	endif
[223]	252	IF (first_appel) THEN
	253	first_appel=.false.
	254	!
	255	! initialisation sorties netcdf
	256	!
[353]	257	idayref = day_ini
	258	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
[223]	259	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
	260	DO i = 1, iim
	261	zx_lon(i,1) = rlon(i+1)
	262	zx_lon(i,jjm+1) = rlon(i+1)
	263	ENDDO
	264	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
	265	cldebug='sous_index'
[353]	266	CALL histbeg(cldebug, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),
	267	$ 1,iim,1,jjm
	268	$ +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)
[223]	269	! no vertical axis
	270	cl_surf(1)='ter'
	271	cl_surf(2)='lic'
	272	cl_surf(3)='oce'
	273	cl_surf(4)='sic'
	274	DO nsrf=1,nbsrf
	275	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,
	276	$ jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	277	END DO
	278	CALL histend(nidbg)
	279	CALL histsync(nidbg)
	280	ENDIF
	281
[2]	282	DO k = 1, klev ! epaisseur de couche
	283	DO i = 1, klon
	284	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	285	ENDDO
	286	ENDDO
	287	DO i = 1, klon ! vent de la premiere couche
	288	ccc zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
	289	zx_alf1 = 1.0
	290	zx_alf2 = 1.0 - zx_alf1
	291	u1lay(i) = u(i,1)zx_alf1 + u(i,2)zx_alf2
	292	v1lay(i) = v(i,1)zx_alf1 + v(i,2)zx_alf2
	293	ENDDO
	294	c
	295	c initialisation:
	296	c
	297	DO i = 1, klon
	298	rugmer(i) = 0.0
	299	cdragh(i) = 0.0
	300	cdragm(i) = 0.0
	301	dflux_t(i) = 0.0
	302	dflux_q(i) = 0.0
	303	zu1(i) = 0.0
	304	zv1(i) = 0.0
	305	ENDDO
[156]	306	ypct = 0.0
	307	yts = 0.0
	308	ysnow = 0.0
[438]	309	yqsurf = 0.0
[156]	310	yalb = 0.0
[280]	311	yalblw = 0.0
[156]	312	yrain_f = 0.0
	313	ysnow_f = 0.0
	314	yfder = 0.0
	315	ytaux = 0.0
	316	ytauy = 0.0
	317	ysolsw = 0.0
	318	ysollw = 0.0
[177]	319	ysollwdown = 0.0
[156]	320	yrugos = 0.0
	321	yu1 = 0.0
	322	yv1 = 0.0
	323	yrads = 0.0
	324	ypaprs = 0.0
	325	ypplay = 0.0
	326	ydelp = 0.0
	327	yu = 0.0
	328	yv = 0.0
	329	yt = 0.0
	330	yq = 0.0
	331	pctsrf_new = 0.0
	332	y_flux_u = 0.0
	333	y_flux_v = 0.0
[290]	334	C$$ PB
	335	y_dflux_t = 0.0
	336	y_dflux_q = 0.0
[235]	337	ytsoil = 999999.
[290]	338	yrugoro = 0.
[156]	339
[2]	340	DO nsrf = 1, nbsrf
	341	DO i = 1, klon
	342	d_ts(i,nsrf) = 0.0
	343	ENDDO
[98]	344	END DO
	345	C§§§ PB
[177]	346	yfluxlat=0.
[98]	347	flux_t = 0.
	348	flux_q = 0.
	349	flux_u = 0.
	350	flux_v = 0.
[2]	351	DO k = 1, klev
	352	DO i = 1, klon
	353	d_t(i,k) = 0.0
	354	d_q(i,k) = 0.0
[98]	355	c$$$ flux_t(i,k) = 0.0
	356	c$$$ flux_q(i,k) = 0.0
[2]	357	d_u(i,k) = 0.0
	358	d_v(i,k) = 0.0
[98]	359	c$$$ flux_u(i,k) = 0.0
	360	c$$$ flux_v(i,k) = 0.0
[2]	361	zcoefh(i,k) = 0.0
	362	ENDDO
	363	ENDDO
	364	cAA IF (itr.GE.1) THEN
	365	cAA DO it = 1, itr
	366	cAA DO k = 1, klev
	367	cAA DO i = 1, klon
	368	cAA d_tr(i,k,it) = 0.0
	369	cAA ENDDO
	370	cAA ENDDO
	371	cAA ENDDO
	372	cAA ENDIF
[150]	373
[2]	374	c
	375	c Boucler sur toutes les sous-fractions du sol:
	376	c
[105]	377	C Initialisation des "pourcentages potentiels". On considere ici qu'on
	378	C peut avoir potentiellementdela glace sur tout le domaine oceanique
	379	C (a affiner)
	380
	381	pctsrf_pot = pctsrf
[177]	382	pctsrf_pot(:,is_oce) = 1. - zmasq(:)
	383	pctsrf_pot(:,is_sic) = 1. - zmasq(:)
[105]	384
[2]	385	DO 99999 nsrf = 1, nbsrf
[98]	386
[2]	387	c chercher les indices:
	388	DO j = 1, klon
	389	ni(j) = 0
	390	ENDDO
	391	knon = 0
	392	DO i = 1, klon
[105]	393
	394	C pour determiner le domaine a traiter on utilise les surfaces "potentielles"
	395	C
	396	IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
[2]	397	knon = knon + 1
	398	ni(knon) = i
	399	ENDIF
	400	ENDDO
	401	c
[438]	402	if (check) THEN
	403	write(,)'CLMAIN, nsrf, knon =',nsrf, knon
	404	call flush(6)
	405	endif
[223]	406	c
	407	c variables pour avoir une sortie IOIPSL des INDEX
	408	c
	409	IF (debugindex) THEN
	410	tabindx(:)=0.
	411	c tabindx(1:knon)=(/FLOAT(i),i=1:knon/)
	412	DO i=1,knon
	413	tabindx(1:knon)=FLOAT(i)
	414	END DO
	415	debugtab(:,:)=0.
	416	ndexbg(:)=0
	417	CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)
	418	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)
	419	$ ,ndexbg)
	420	ENDIF
[2]	421	IF (knon.EQ.0) GOTO 99999
	422	DO j = 1, knon
	423	i = ni(j)
	424	ypct(j) = pctsrf(i,nsrf)
	425	yts(j) = ts(i,nsrf)
[98]	426	ysnow(j) = snow(i,nsrf)
[438]	427	yqsurf(j) = qsurf(i,nsrf)
[98]	428	yalb(j) = albe(i,nsrf)
[280]	429	yalblw(j) = alblw(i,nsrf)
[98]	430	yrain_f(j) = rain_f(i)
	431	ysnow_f(j) = snow_f(i)
[258]	432	yagesno(j) = agesno(i,nsrf)
[101]	433	yfder(j) = fder(i)
[102]	434	ytaux(j) = flux_u(i,1,nsrf)
	435	ytauy(j) = flux_v(i,1,nsrf)
[433]	436	ysolsw(j) = solsw(i,nsrf)
	437	ysollw(j) = sollw(i,nsrf)
[177]	438	ysollwdown(j) = sollwdown(i)
[2]	439	yrugos(j) = rugos(i,nsrf)
[171]	440	yrugoro(j) = rugoro(i)
[2]	441	yu1(j) = u1lay(i)
	442	yv1(j) = v1lay(i)
[433]	443	yrads(j) = ysolsw(j)+ ysollw(j)
[2]	444	ypaprs(j,klev+1) = paprs(i,klev+1)
[522]	445	y_run_off_lic_0(j) = run_off_lic_0(i)
[177]	446	END DO
[454]	447	C
	448	C IF bucket model for continent, copy soil water content
	449	IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN
[438]	450	DO j = 1, knon
	451	i = ni(j)
	452	yqsol(j) = qsol(i)
	453	END DO
	454	ELSE
	455	yqsol(:)=0.
	456	ENDIF
[177]	457	c$$$ PB ajour pour soil
	458	DO k = 1, nsoilmx
	459	DO j = 1, knon
	460	i = ni(j)
	461	ytsoil(j,k) = ftsoil(i,k,nsrf)
	462	END DO
	463	END DO
[2]	464	DO k = 1, klev
	465	DO j = 1, knon
	466	i = ni(j)
	467	ypaprs(j,k) = paprs(i,k)
	468	ypplay(j,k) = pplay(i,k)
	469	ydelp(j,k) = delp(i,k)
	470	yu(j,k) = u(i,k)
	471	yv(j,k) = v(i,k)
	472	yt(j,k) = t(i,k)
	473	yq(j,k) = q(i,k)
	474	ENDDO
	475	ENDDO
	476	c
	477	c
	478	c calculer Cdrag et les coefficients d'echange
	479	CALL coefkz(nsrf, knon, ypaprs, ypplay,
[486]	480	cIM 261103
	481	. ksta, ksta_ter,
	482	cIM 261103
[2]	483	. yts, yrugos, yu, yv, yt, yq,
[438]	484	. yqsurf,
[2]	485	. ycoefm, ycoefh)
[143]	486	CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
[39]	487	. ycoefm0, ycoefh0)
	488	DO k = 1, klev
	489	DO i = 1, knon
	490	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
	491	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
	492	ENDDO
	493	ENDDO
[2]	494	c
[467]	495	cIM cf JLD : on seuille ycoefm et ycoefh
	496	if (nsrf.eq.is_oce) then
	497	do j=1,knon
	498	c ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
	499	ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
	500	c ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
	501	ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
	502	enddo
	503	endif
	504
[486]	505	c
	506	cIM: 261103
	507	if (ok_kzmin) THEN
	508	cIM cf FH: 201103 BEG
	509	c Calcul d'une diffusion minimale pour les conditions tres stables.
	510	call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
	511	. ,ycoefm0,ycoefh0)
	512	c call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ ')
	513	c call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN ')
	514
	515	if ( 1.eq.1 ) then
	516	DO k = 1, klev
	517	DO i = 1, knon
	518	ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
	519	ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
	520	ENDDO
	521	ENDDO
	522	endif
	523	cIM cf FH: 201103 END
	524	endif !ok_kzmin
	525	cIM: 261103
	526
	527	c
[105]	528	c calculer la diffusion des vitesses "u" et "v"
	529	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
	530	s y_d_u,y_flux_u)
	531	CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
	532	s y_d_v,y_flux_v)
	533
	534	c pour le couplage
	535	ytaux = y_flux_u(:,1)
	536	ytauy = y_flux_v(:,1)
[150]	537
[363]	538	c FH modif sur le cdrag temperature
[467]	539	c$$$PB : déplace dans clcdrag
	540	c$$$ do i=1,knon
	541	c$$$ ycoefh(i,1)=ycoefm(i,1)*0.8
	542	c$$$ enddo
[363]	543
[2]	544	c calculer la diffusion de "q" et de "h"
[205]	545	CALL clqh(dtime, itap, date0,jour, debut,lafin,
[171]	546	e rlon, rlat, cufi, cvfi,
[112]	547	e knon, nsrf, ni, pctsrf,
[438]	548	e soil_model, ytsoil,yqsol,
[112]	549	e ok_veget, ocean, npas, nexca,
[475]	550	e rmu0, co2_ppm, yrugos, yrugoro,
[109]	551	e yu1, yv1, ycoefh,
	552	e yt,yq,yts,ypaprs,ypplay,
[438]	553	e ydelp,yrads,yalb, yalblw, ysnow, yqsurf,
[109]	554	e yrain_f, ysnow_f, yfder, ytaux, ytauy,
[177]	555	c$$$ e ysollw, ysolsw,
	556	e ysollw, ysollwdown, ysolsw,yfluxlat,
[258]	557	s pctsrf_new, yagesno,
[109]	558	s y_d_t, y_d_q, y_d_ts, yz0_new,
[454]	559	cIM cf JLD s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
	560	s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
[522]	561	s y_fqcalving,y_ffonte,y_run_off_lic_0)
[2]	562	c
	563	c calculer la longueur de rugosite sur ocean
[290]	564	yrugm=0.
[2]	565	IF (nsrf.EQ.is_oce) THEN
	566	DO j = 1, knon
[467]	567	yrugm(j) = 0.018ycoefm(j,1) (yu1(j)2+yv1(j)2)/RG
	568	$ + 0.1114e-6 / sqrt(ycoefm(j,1) (yu1(j)2+yv1(j)2))
[2]	569	yrugm(j) = MAX(1.5e-05,yrugm(j))
	570	ENDDO
	571	ENDIF
	572	DO j = 1, knon
	573	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
	574	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
	575	yu1(j) = yu1(j) * ypct(j)
	576	yv1(j) = yv1(j) * ypct(j)
	577	ENDDO
	578	c
	579	DO k = 1, klev
[98]	580	DO j = 1, knon
	581	i = ni(j)
	582	ycoefh(j,k) = ycoefh(j,k) * ypct(j)
	583	ycoefm(j,k) = ycoefm(j,k) * ypct(j)
	584	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
	585	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
	586	C§§§ PB
	587	flux_t(i,k,nsrf) = y_flux_t(j,k)
	588	flux_q(i,k,nsrf) = y_flux_q(j,k)
	589	flux_u(i,k,nsrf) = y_flux_u(j,k)
	590	flux_v(i,k,nsrf) = y_flux_v(j,k)
	591	c$$$ PB y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
	592	c$$$ PB y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
	593	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
	594	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
	595	c$$$ PB y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
	596	c$$$ PB y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
	597	ENDDO
[2]	598	ENDDO
[98]	599
[290]	600
[98]	601	evap(:,nsrf) = - flux_q(:,1,nsrf)
[2]	602	c
[235]	603	albe(:, nsrf) = 0.
[280]	604	alblw(:, nsrf) = 0.
[235]	605	snow(:, nsrf) = 0.
[438]	606	qsurf(:, nsrf) = 0.
[235]	607	rugos(:, nsrf) = 0.
	608	fluxlat(:,nsrf) = 0.
[2]	609	DO j = 1, knon
[98]	610	i = ni(j)
[2]	611	d_ts(i,nsrf) = y_d_ts(j)
[98]	612	albe(i,nsrf) = yalb(j)
[280]	613	alblw(i,nsrf) = yalblw(j)
[98]	614	snow(i,nsrf) = ysnow(j)
[438]	615	qsurf(i,nsrf) = yqsurf(j)
[109]	616	rugos(i,nsrf) = yz0_new(j)
[177]	617	fluxlat(i,nsrf) = yfluxlat(j)
	618	c$$$ pb rugmer(i) = yrugm(j)
[235]	619	IF (nsrf .EQ. is_oce) then
[400]	620	rugmer(i) = yrugm(j)
	621	rugos(i,nsrf) = yrugm(j)
	622	endif
[454]	623	cIM cf JLD ??
	624	fqcalving(i,nsrf) = y_fqcalving(j)
	625	ffonte(i,nsrf) = y_ffonte(j)
[2]	626	cdragh(i) = cdragh(i) + ycoefh(j,1)
	627	cdragm(i) = cdragm(i) + ycoefm(j,1)
	628	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
	629	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
	630	zu1(i) = zu1(i) + yu1(j)
	631	zv1(i) = zv1(i) + yv1(j)
[177]	632	END DO
[438]	633	IF ( nsrf .eq. is_ter ) THEN
	634	DO j = 1, knon
	635	i = ni(j)
	636	qsol(i) = yqsol(j)
	637	END DO
	638	END IF
[522]	639	IF ( nsrf .eq. is_lic ) THEN
	640	DO j = 1, knon
	641	i = ni(j)
	642	run_off_lic_0(i) = y_run_off_lic_0(j)
	643	END DO
	644	END IF
[177]	645	c$$$ PB ajout pour soil
[235]	646	ftsoil(:,:,nsrf) = 0.
[177]	647	DO k = 1, nsoilmx
	648	DO j = 1, knon
	649	i = ni(j)
	650	ftsoil(i, k, nsrf) = ytsoil(j,k)
	651	END DO
	652	END DO
[2]	653	c
	654	#ifdef CRAY
	655	DO k = 1, klev
	656	DO j = 1, knon
	657	i = ni(j)
	658	#else
	659	DO j = 1, knon
	660	i = ni(j)
	661	DO k = 1, klev
	662	#endif
	663	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
	664	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
[98]	665	c$$$ PB flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
	666	c$$$ flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
[2]	667	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
	668	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
[98]	669	c$$$ PB flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
	670	c$$$ flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
[2]	671	zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
	672	ENDDO
	673	ENDDO
	674	c
[407]	675	c
	676	#undef T2m
[433]	677	#define T2m
[407]	678	#ifdef T2m
	679	ccc diagnostic t,q a 2m et u, v a 10m
	680	c
	681	DO j=1, knon
	682	i = ni(j)
	683	uzon(j) = yu(j,1) + y_d_u(j,1)
	684	vmer(j) = yv(j,1) + y_d_v(j,1)
	685	tair1(j) = yt(j,1) + y_d_t(j,1)
	686	qair1(j) = yq(j,1) + y_d_q(j,1)
	687	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
	688	& * (ypaprs(j,1)-ypplay(j,1))
	689	tairsol(j) = yts(j) + y_d_ts(j)
	690	rugo1(j) = yrugos(j)
	691	IF(nsrf.EQ.is_oce) THEN
	692	rugo1(j) = rugos(i,nsrf)
	693	ENDIF
	694	psfce(j)=ypaprs(j,1)
	695	patm(j)=ypplay(j,1)
	696	c
[438]	697	qairsol(j) = yqsurf(j)
	698	c$$$ IF (nsrf.EQ.1) THEN
	699	c$$$ qairsol(j) = yqsurf(j)
	700	c$$$ ELSE IF(nsrf.GT.1) THEN
	701	c$$$ zt = ts(i,nsrf)
	702	c$$$ IF (thermcep) THEN
	703	c$$$ zdelta = MAX(0.,SIGN(1.,RTT-zt))
	704	c$$$ zqs = R2ES * FOEEW(zt,zdelta) / ypplay(j,1)
	705	c$$$ zqs = MIN(0.5,zqs)
	706	c$$$ zcor = 1./(1.-RETV*zqs)
	707	c$$$ zqs = zqs*zcor
	708	c$$$ ELSE
	709	c$$$ IF (zt .LT. t_coup) THEN
	710	c$$$ zqs = qsats(zt) / ypplay(j,1)
	711	c$$$ ELSE
	712	c$$$ zqs = qsatl(zt) / ypplay(j,1)
	713	c$$$ ENDIF
	714	c$$$ ENDIF
	715	c$$$ qairsol(j) = zqs
	716	c$$$ ENDIF
[407]	717	ENDDO
	718	c
[438]	719	if (check) THEN
	720	WRITE(,)' avant stdlevvar. nsrf=',nsrf
	721	IF(nsrf.EQ.3) THEN
	722	j=1465
	723	WRITE(,)' INstO',klon,knon,nsrf,zxli,uzon(j),vmer(j),
	724	& tair1(j),qair1(j),zgeo1(j),tairsol(j),qairsol(j),rugo1(j),
	725	& psfce(j),patm(j)
	726	ENDIF
	727	WRITE(,)' qairsol (min, max)'
	728	$ , minval(qairsol(1:knon)), maxval(qairsol(1:knon))
	729	call flush(6)
[407]	730	ENDIF
	731	c
	732	CALL stdlevvar(klon, knon, nsrf, zxli,
	733	& uzon, vmer, tair1, qair1, zgeo1,
	734	& tairsol, qairsol, rugo1, psfce, patm,
	735	& yt2m, yq2m, yu10m)
	736
	737	c
[438]	738	if (check) THEN
[407]	739	IF(nsrf.EQ.3) THEN
	740	j=1465
	741	WRITE(,)' OUstd',klon,knon,nsrf,zxli,uzon(j),vmer(j),
	742	& tair1(j),qair1(j),zgeo1(j),tairsol(j),qairsol(j),rugo1(j),
	743	& psfce(j),patm(j)
	744	WRITE(,)' tqu',yt2m(j),yq2m(j),yu10m(j)
[438]	745	call flush(6)
[407]	746	ENDIF
[438]	747	ENDIF
[407]	748	c
	749	DO j=1, knon
	750	i = ni(j)
	751	t2m(i,nsrf)=yt2m(j)
	752
[438]	753	if (check) THEN
	754	IF(nsrf.EQ.3 .and. j.EQ.1465) THEN
[407]	755	WRITE(,) 't2m APRES stdlev',j,i,tair1(j),t2m(i,nsrf),
	756	$ tairsol(j),rlon(i),rlat(i)
[438]	757	call flush(6)
[407]	758	ENDIF
[438]	759	ENDIF
[407]	760	c
	761	q2m(i,nsrf)=yq2m(j)
	762	c
	763	c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
	764	u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)2+vmer(j)2)
	765	v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)2+vmer(j)2)
	766	c
	767	ENDDO
	768	#else
	769	DO j=1, knon
[486]	770	i = ni(j)
[407]	771	t2m(i,nsrf)=0.
	772	q2m(i,nsrf)=0.
	773	u10m(i,nsrf)=0.
	774	v10m(i,nsrf)=0.
	775	ENDDO
	776	#endif
[2]	777	99999 CONTINUE
	778	c
[112]	779	C
	780	C On utilise les nouvelles surfaces
	781	C A rajouter: conservation de l'albedo
	782	C
[171]	783	rugos(:,is_oce) = rugmer
[112]	784	pctsrf = pctsrf_new
[98]	785
[2]	786	RETURN
	787	END
[205]	788	SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin,
[171]	789	e rlon, rlat, cufi, cvfi,
[112]	790	e knon, nisurf, knindex, pctsrf,
[438]	791	$ soil_model,tsoil,qsol,
[112]	792	e ok_veget, ocean, npas, nexca,
[475]	793	e rmu0, co2_ppm, rugos, rugoro,
[98]	794	e u1lay,v1lay,coef,
[86]	795	e t,q,ts,paprs,pplay,
[438]	796	e delp,radsol,albedo,alblw,snow,qsurf,
[101]	797	e precip_rain, precip_snow, fder, taux, tauy,
[293]	798	$ sollw, sollwdown, swnet,fluxlat,
[112]	799	s pctsrf_new, agesno,
[109]	800	s d_t, d_q, d_ts, z0_new,
[454]	801	cIM cf JLD s flux_t, flux_q,dflux_s,dflux_l)
	802	s flux_t, flux_q,dflux_s,dflux_l,
[522]	803	s fqcalving,ffonte,run_off_lic_0)
[86]	804
	805	USE interface_surf
	806
[2]	807	IMPLICIT none
	808	c======================================================================
	809	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	810	c Objet: diffusion verticale de "q" et de "h"
	811	c======================================================================
	812	#include "dimensions.h"
	813	#include "dimphy.h"
[98]	814	#include "YOMCST.h"
	815	#include "YOETHF.h"
	816	#include "FCTTRE.h"
	817	#include "indicesol.h"
[177]	818	#include "dimsoil.h"
[475]	819
[2]	820	c Arguments:
	821	INTEGER knon
	822	REAL dtime ! intervalle du temps (s)
[205]	823	real date0
[2]	824	REAL u1lay(klon) ! vitesse u de la 1ere couche (m/s)
	825	REAL v1lay(klon) ! vitesse v de la 1ere couche (m/s)
	826	REAL coef(klon,klev) ! le coefficient d'echange (m**2/s)
	827	c multiplie par le cisaillement du
	828	c vent (dV/dz); la premiere valeur
	829	c indique la valeur de Cdrag (sans unite)
	830	REAL t(klon,klev) ! temperature (K)
	831	REAL q(klon,klev) ! humidite specifique (kg/kg)
	832	REAL ts(klon) ! temperature du sol (K)
[98]	833	REAL evap(klon) ! evaporation au sol
[2]	834	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
	835	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
	836	REAL delp(klon,klev) ! epaisseur de couche en pression (Pa)
	837	REAL radsol(klon) ! ray. net au sol (Solaire+IR) W/m2
[98]	838	REAL albedo(klon) ! albedo de la surface
[280]	839	REAL alblw(klon)
[98]	840	REAL snow(klon) ! hauteur de neige
[438]	841	REAL qsurf(klon) ! humidite de l'air au dessus de la surface
[98]	842	real precip_rain(klon), precip_snow(klon)
[112]	843	REAL agesno(klon)
[171]	844	REAL rugoro(klon)
[522]	845	REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
[109]	846	integer jour ! jour de l'annee en cours
	847	real rmu0(klon) ! cosinus de l'angle solaire zenithal
[156]	848	real rugos(klon) ! rugosite
[98]	849	integer knindex(klon)
	850	real pctsrf(klon,nbsrf)
[171]	851	real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
[112]	852	logical ok_veget
[475]	853	REAL co2_ppm ! taux CO2 atmosphere
[112]	854	character*6 ocean
	855	integer npas, nexca
	856
[2]	857	c
	858	REAL d_t(klon,klev) ! incrementation de "t"
	859	REAL d_q(klon,klev) ! incrementation de "q"
	860	REAL d_ts(klon) ! incrementation de "ts"
	861	REAL flux_t(klon,klev) ! (diagnostic) flux de la chaleur
	862	c sensible, flux de Cp*T, positif vers
	863	c le bas: j/(m**2 s) c.a.d.: W/m2
	864	REAL flux_q(klon,klev) ! flux de la vapeur d'eau:kg/(m**2 s)
	865	REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
	866	REAL dflux_l(klon) ! derivee du flux latent dF/dTs
[454]	867	cIM cf JLD
	868	c Flux thermique utiliser pour fondre la neige
	869	REAL ffonte(klon)
	870	c Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
	871	c hauteur de neige, en kg/m2/s
	872	REAL fqcalving(klon)
[2]	873	c======================================================================
	874	REAL t_grnd ! temperature de rappel pour glace de mer
	875	PARAMETER (t_grnd=271.35)
	876	REAL t_coup
	877	PARAMETER(t_coup=273.15)
	878	c======================================================================
	879	INTEGER i, k
	880	REAL zx_cq(klon,klev)
	881	REAL zx_dq(klon,klev)
	882	REAL zx_ch(klon,klev)
	883	REAL zx_dh(klon,klev)
	884	REAL zx_buf1(klon)
	885	REAL zx_buf2(klon)
	886	REAL zx_coef(klon,klev)
	887	REAL local_h(klon,klev) ! enthalpie potentielle
	888	REAL local_q(klon,klev)
	889	REAL local_ts(klon)
	890	REAL psref(klon) ! pression de reference pour temperature potent.
	891	REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
	892	c======================================================================
	893	c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
	894	REAL gamq(klon,2:klev)
	895	c contre-gradient pour la chaleur sensible: Kelvin/metre
	896	REAL gamt(klon,2:klev)
	897	REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
	898	REAL zdelz
	899	c======================================================================
[98]	900	logical contreg
	901	parameter (contreg=.true.)
[2]	902	c======================================================================
[86]	903	c Rajout pour l'interface
	904	integer itime
	905	integer nisurf
[112]	906	logical debut, lafin
[98]	907	real zlev1(klon)
[101]	908	real fder(klon), taux(klon), tauy(klon)
[86]	909	real temp_air(klon), spechum(klon)
[177]	910	real epot_air(klon), ccanopy(klon)
[86]	911	real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
	912	real petBcoef(klon), peqBcoef(klon)
[177]	913	real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
[86]	914	real p1lay(klon)
[177]	915	c$$$C PB ajout pour soil
	916	LOGICAL soil_model
	917	REAL tsoil(klon, nsoilmx)
[438]	918	REAL qsol(klon)
[86]	919
	920	! Parametres de sortie
[98]	921	real fluxsens(klon), fluxlat(klon)
[86]	922	real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
	923	real emis_new(klon), z0_new(klon)
	924	real pctsrf_new(klon,nbsrf)
[393]	925	c JLD
	926	real zzpk
[438]	927	C
	928	character (len = 20) :: modname = 'Debut clqh'
	929	LOGICAL check
[458]	930	PARAMETER (check=.false.)
[438]	931	C
	932	if (check) THEN
	933	write(,) modname,' nisurf=',nisurf
	934	call flush(6)
	935	endif
[86]	936	c
[438]	937	if (check) THEN
	938	WRITE(,)' qsurf (min, max)'
	939	$ , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
	940	call flush(6)
	941	ENDIF
	942	C
[98]	943	if (.not. contreg) then
	944	do k = 2, klev
	945	do i = 1, knon
	946	gamq(i,k) = 0.0
	947	gamt(i,k) = 0.0
	948	enddo
	949	enddo
	950	else
	951	do k = 3, klev
	952	do i = 1, knon
	953	gamq(i,k)= 0.0
	954	gamt(i,k)= -1.0e-03
	955	enddo
	956	enddo
	957	do i = 1, knon
	958	gamq(i,2) = 0.0
	959	gamt(i,2) = -2.5e-03
	960	enddo
	961	endif
	962
[546]	963	psref = 0.
[2]	964	DO i = 1, knon
	965	psref(i) = paprs(i,1) !pression de reference est celle au sol
	966	local_ts(i) = ts(i)
	967	ENDDO
	968	DO k = 1, klev
	969	DO i = 1, knon
	970	zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
	971	zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
	972	local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
	973	local_q(i,k) = q(i,k)
	974	ENDDO
	975	ENDDO
	976	c
	977	c Convertir les coefficients en variables convenables au calcul:
	978	c
	979	c
	980	DO k = 2, klev
	981	DO i = 1, knon
	982	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
	983	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
	984	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
	985	ENDDO
	986	ENDDO
	987	c
	988	c Preparer les flux lies aux contre-gardients
	989	c
	990	DO k = 2, klev
	991	DO i = 1, knon
	992	zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
	993	. *(pplay(i,k-1)-pplay(i,k))
	994	z_gamaq(i,k) = gamq(i,k) * zdelz
	995	z_gamah(i,k) = gamt(i,k) * zdelz RCPD zx_pkh(i,k)
	996	ENDDO
	997	ENDDO
	998	DO i = 1, knon
	999	zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
	1000	zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
	1001	. -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
	1002	zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
	1003	c
[393]	1004	zzpk=(pplay(i,klev)/psref(i))**RKAPPA
	1005	zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
	1006	zx_ch(i,klev) = (local_h(i,klev)zzpkdelp(i,klev)
[2]	1007	. -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
	1008	zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
	1009	ENDDO
	1010	DO k = klev-1, 2 , -1
	1011	DO i = 1, knon
	1012	zx_buf1(i) = delp(i,k)+zx_coef(i,k)
	1013	. +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
	1014	zx_cq(i,k) = (local_q(i,k)*delp(i,k)
	1015	. +zx_coef(i,k+1)*zx_cq(i,k+1)
	1016	. +zx_coef(i,k+1)*z_gamaq(i,k+1)
	1017	. -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
	1018	zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
	1019	c
[393]	1020	zzpk=(pplay(i,k)/psref(i))**RKAPPA
	1021	zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k)
[2]	1022	. +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
[393]	1023	zx_ch(i,k) = (local_h(i,k)zzpkdelp(i,k)
[2]	1024	. +zx_coef(i,k+1)*zx_ch(i,k+1)
	1025	. +zx_coef(i,k+1)*z_gamah(i,k+1)
	1026	. -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
	1027	zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
	1028	ENDDO
	1029	ENDDO
[84]	1030	C
	1031	C nouvelle formulation JL Dufresne
	1032	C
	1033	C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
	1034	C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
	1035	C
[2]	1036	DO i = 1, knon
	1037	zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
	1038	zx_cq(i,1) = (local_q(i,1)*delp(i,1)
[84]	1039	. +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
[2]	1040	. /zx_buf1(i)
[84]	1041	zx_dq(i,1) = -1. * RG / zx_buf1(i)
[2]	1042	c
[393]	1043	zzpk=(pplay(i,1)/psref(i))**RKAPPA
	1044	zx_buf2(i) = zzpkdelp(i,1) + zx_coef(i,2)(1.-zx_dh(i,2))
	1045	zx_ch(i,1) = (local_h(i,1)zzpkdelp(i,1)
[84]	1046	. +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
[2]	1047	. /zx_buf2(i)
[84]	1048	zx_dh(i,1) = -1. * RG / zx_buf2(i)
[2]	1049	ENDDO
[84]	1050
[86]	1051	C Appel a interfsurf (appel generique) routine d'interface avec la surface
[84]	1052
[290]	1053	c initialisation
[546]	1054	petAcoef =0.
[290]	1055	peqAcoef = 0.
	1056	petBcoef =0.
	1057	peqBcoef = 0.
	1058	p1lay =0.
[546]	1059	zlev1(:) = 0.
	1060	temp_air(:) = 0.
	1061	spechum(:) = 0.
	1062	epot_air(:) = 0.
	1063	tq_cdrag(:) = 0.
[290]	1064
[171]	1065	c do i = 1, knon
[290]	1066	petAcoef(1:knon) = zx_ch(1:knon,1)
	1067	peqAcoef(1:knon) = zx_cq(1:knon,1)
	1068	petBcoef(1:knon) = zx_dh(1:knon,1)
	1069	peqBcoef(1:knon) = zx_dq(1:knon,1)
	1070	tq_cdrag(1:knon) =coef(1:knon,1)
	1071	temp_air(1:knon) =t(1:knon,1)
	1072	epot_air(1:knon) =local_h(1:knon,1)
	1073	spechum(1:knon)=q(1:knon,1)
	1074	p1lay(1:knon) = pplay(1:knon,1)
	1075	zlev1(1:knon) = delp(1:knon,1)
[293]	1076	c swnet = swdown * (1. - albedo)
	1077	swdown(1:knon) = swnet(1:knon)
[171]	1078	c enddo
[473]	1079	ccanopy = co2_ppm
[86]	1080
[205]	1081	CALL interfsurf(itime, dtime, date0, jour, rmu0,
[177]	1082	e klon, iim, jjm, nisurf, knon, knindex, pctsrf,
	1083	e rlon, rlat, cufi, cvfi,
[438]	1084	e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol,
[177]	1085	e zlev1, u1lay, v1lay, temp_air, spechum, epot_air, ccanopy,
[109]	1086	e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
[177]	1087	e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,
[171]	1088	e fder, taux, tauy, rugos, rugoro,
[438]	1089	e albedo, snow, qsurf,
[109]	1090	e ts, p1lay, psref, radsol,
[112]	1091	e ocean, npas, nexca, zmasq,
[109]	1092	s evap, fluxsens, fluxlat, dflux_l, dflux_s,
[280]	1093	s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new,
[454]	1094	cIM cf JLD s pctsrf_new, agesno)
[522]	1095	s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0)
[86]	1096
[177]	1097
[150]	1098	do i = 1, knon
[171]	1099	flux_t(i,1) = fluxsens(i)
	1100	flux_q(i,1) = - evap(i)
	1101	d_ts(i) = tsurf_new(i) - ts(i)
	1102	albedo(i) = alb_new(i)
[150]	1103	enddo
[86]	1104
[2]	1105	c==== une fois on a zx_h_ts, on peut faire l'iteration ========
	1106	DO i = 1, knon
[86]	1107	local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)flux_t(i,1)dtime
	1108	local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)flux_q(i,1)dtime
[2]	1109	ENDDO
	1110	DO k = 2, klev
	1111	DO i = 1, knon
	1112	local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
	1113	local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
	1114	ENDDO
	1115	ENDDO
	1116	c======================================================================
	1117	c== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
	1118	c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
	1119	DO k = 2, klev
	1120	DO i = 1, knon
	1121	flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
	1122	. * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
	1123	flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
	1124	. * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
	1125	. / zx_pkh(i,k)
	1126	ENDDO
	1127	ENDDO
	1128	c======================================================================
[84]	1129	C Calcul tendances
[2]	1130	DO k = 1, klev
	1131	DO i = 1, knon
	1132	d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
	1133	d_q(i,k) = local_q(i,k) - q(i,k)
	1134	ENDDO
	1135	ENDDO
	1136	c
[148]	1137
[2]	1138	RETURN
	1139	END
	1140	SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
	1141	e paprs,pplay,delp,
	1142	s d_ven,flux_v)
	1143	IMPLICIT none
	1144	c======================================================================
	1145	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
	1146	c Objet: diffusion vertical de la vitesse "ven"
	1147	c======================================================================
	1148	c Arguments:
	1149	c dtime----input-R- intervalle du temps (en second)
	1150	c u1lay----input-R- vent u de la premiere couche (m/s)
	1151	c v1lay----input-R- vent v de la premiere couche (m/s)
	1152	c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
	1153	c le cisaillement du vent (dV/dz); la premiere
	1154	c valeur indique la valeur de Cdrag (sans unite)
	1155	c t--------input-R- temperature (K)
	1156	c ven------input-R- vitesse horizontale (m/s)
	1157	c paprs----input-R- pression a inter-couche (Pa)
	1158	c pplay----input-R- pression au milieu de couche (Pa)
	1159	c delp-----input-R- epaisseur de couche (Pa)
	1160	c
	1161	c
	1162	c d_ven----output-R- le changement de "ven"
	1163	c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
	1164	c======================================================================
	1165	#include "dimensions.h"
	1166	#include "dimphy.h"
	1167	INTEGER knon
	1168	REAL dtime
	1169	REAL u1lay(klon), v1lay(klon)
	1170	REAL coef(klon,klev)
	1171	REAL t(klon,klev), ven(klon,klev)
	1172	REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
	1173	REAL d_ven(klon,klev)
	1174	REAL flux_v(klon,klev)
	1175	c======================================================================
	1176	#include "YOMCST.h"
	1177	c======================================================================
	1178	INTEGER i, k
	1179	REAL zx_cv(klon,2:klev)
	1180	REAL zx_dv(klon,2:klev)
	1181	REAL zx_buf(klon)
	1182	REAL zx_coef(klon,klev)
	1183	REAL local_ven(klon,klev)
	1184	REAL zx_alf1(klon), zx_alf2(klon)
	1185	c======================================================================
	1186	DO k = 1, klev
	1187	DO i = 1, knon
	1188	local_ven(i,k) = ven(i,k)
	1189	ENDDO
	1190	ENDDO
	1191	c======================================================================
	1192	DO i = 1, knon
	1193	ccc zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
	1194	zx_alf1(i) = 1.0
	1195	zx_alf2(i) = 1.0 - zx_alf1(i)
	1196	zx_coef(i,1) = coef(i,1)
	1197	. * (1.0+SQRT(u1lay(i)2+v1lay(i)2))
	1198	. * pplay(i,1)/(RD*t(i,1))
	1199	zx_coef(i,1) = zx_coef(i,1) * dtime*RG
	1200	ENDDO
	1201	c======================================================================
	1202	DO k = 2, klev
	1203	DO i = 1, knon
	1204	zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
	1205	. (paprs(i,k)2/(t(i,k)+t(i,k-1))/RD)**2
	1206	zx_coef(i,k) = zx_coef(i,k) * dtime*RG
	1207	ENDDO
	1208	ENDDO
	1209	c======================================================================
	1210	DO i = 1, knon
	1211	zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
	1212	zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
	1213	zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
	1214	. /zx_buf(i)
	1215	ENDDO
	1216	DO k = 3, klev
	1217	DO i = 1, knon
	1218	zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
	1219	. + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
	1220	zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
	1221	. +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
	1222	zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
	1223	ENDDO
	1224	ENDDO
	1225	DO i = 1, knon
	1226	local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
	1227	. +zx_coef(i,klev)*zx_cv(i,klev) )
	1228	. / ( delp(i,klev) + zx_coef(i,klev)
	1229	. -zx_coef(i,klev)*zx_dv(i,klev) )
	1230	ENDDO
	1231	DO k = klev-1, 1, -1
	1232	DO i = 1, knon
	1233	local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
	1234	ENDDO
	1235	ENDDO
	1236	c======================================================================
	1237	c== flux_v est le flux de moment angulaire (positif vers bas)
	1238	c== dont l'unite est: (kg m/s)/(m**2 s)
	1239	DO i = 1, knon
	1240	flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
	1241	. (local_ven(i,1)zx_alf1(i)
	1242	. +local_ven(i,2)*zx_alf2(i))
	1243	ENDDO
	1244	DO k = 2, klev
	1245	DO i = 1, knon
	1246	flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
	1247	. * (local_ven(i,k)-local_ven(i,k-1))
	1248	ENDDO
	1249	ENDDO
	1250	c
	1251	DO k = 1, klev
	1252	DO i = 1, knon
	1253	d_ven(i,k) = local_ven(i,k) - ven(i,k)
	1254	ENDDO
	1255	ENDDO
	1256	c
	1257	RETURN
	1258	END
	1259	SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
[486]	1260	cIM 261103
	1261	. ksta, ksta_ter,
	1262	cIM 261103
[2]	1263	. ts, rugos,
	1264	. u,v,t,q,
[438]	1265	. qsurf,
[2]	1266	. pcfm, pcfh)
	1267	IMPLICIT none
	1268	c======================================================================
	1269	c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
	1270	c (une version strictement identique a l'ancien modele)
	1271	c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
	1272	c coefficients d'echange turbulent dans l'atmosphere.
	1273	c Arguments:
	1274	c nsrf-----input-I- indicateur de la nature du sol
	1275	c knon-----input-I- nombre de points a traiter
	1276	c paprs----input-R- pression a chaque intercouche (en Pa)
	1277	c pplay----input-R- pression au milieu de chaque couche (en Pa)
	1278	c ts-------input-R- temperature du sol (en Kelvin)
	1279	c rugos----input-R- longeur de rugosite (en m)
	1280	c u--------input-R- vitesse u
	1281	c v--------input-R- vitesse v
	1282	c t--------input-R- temperature (K)
	1283	c q--------input-R- vapeur d'eau (kg/kg)
	1284	c
	1285	c itop-----output-I- numero de couche du sommet de la couche limite
	1286	c pcfm-----output-R- coefficients a calculer (vitesse)
	1287	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	1288	c======================================================================
	1289	#include "dimensions.h"
	1290	#include "dimphy.h"
	1291	#include "YOMCST.h"
	1292	#include "indicesol.h"
	1293	c
	1294	c Arguments:
	1295	c
	1296	INTEGER knon, nsrf
	1297	REAL ts(klon)
	1298	REAL paprs(klon,klev+1), pplay(klon,klev)
	1299	REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
	1300	REAL rugos(klon)
	1301	c
	1302	REAL pcfm(klon,klev), pcfh(klon,klev)
	1303	INTEGER itop(klon)
	1304	c
	1305	c Quelques constantes et options:
	1306	c
	1307	REAL cepdu2, ckap, cb, cc, cd, clam
	1308	PARAMETER (cepdu2 =(0.1)**2)
[407]	1309	PARAMETER (CKAP=0.4)
[2]	1310	PARAMETER (cb=5.0)
	1311	PARAMETER (cc=5.0)
	1312	PARAMETER (cd=5.0)
	1313	PARAMETER (clam=160.0)
	1314	REAL ratqs ! largeur de distribution de vapeur d'eau
[39]	1315	PARAMETER (ratqs=0.05)
[2]	1316	LOGICAL richum ! utilise le nombre de Richardson humide
	1317	PARAMETER (richum=.TRUE.)
	1318	REAL ric ! nombre de Richardson critique
	1319	PARAMETER(ric=0.4)
	1320	REAL prandtl
	1321	PARAMETER (prandtl=0.4)
	1322	REAL kstable ! diffusion minimale (situation stable)
[486]	1323	! GKtest
	1324	! PARAMETER (kstable=1.0e-10)
	1325	REAL ksta, ksta_ter
	1326	cIM: 261103 REAL kstable_ter, kstable_sinon
	1327	cIM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6)
	1328	cIM: 261103 PARAMETER (kstable_ter = 1.0e-8)
	1329	cIM: 261103 PARAMETER (kstable_ter = 1.0e-10)
	1330	cIM: 261103 PARAMETER (kstable_sinon = 1.0e-10)
	1331	! fin GKtest
[2]	1332	REAL mixlen ! constante controlant longueur de melange
	1333	PARAMETER (mixlen=35.0)
	1334	INTEGER isommet ! le sommet de la couche limite
	1335	PARAMETER (isommet=klev)
	1336	LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
	1337	PARAMETER (tvirtu=.TRUE.)
	1338	LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
	1339	PARAMETER (opt_ec=.FALSE.)
	1340	LOGICAL contreg ! utiliser le contre-gradient dans Ri
	1341	PARAMETER (contreg=.TRUE.)
	1342	c
	1343	c Variables locales:
	1344	c
	1345	INTEGER i, k
	1346	REAL zgeop(klon,klev)
	1347	REAL zmgeom(klon)
	1348	REAL zri(klon)
	1349	REAL zl2(klon)
[407]	1350
	1351	REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
	1352	REAL pcfm1(klon), pcfh1(klon)
	1353	c
	1354	REAL zdphi, zdu2, ztvd, ztvu, zcdn
	1355	REAL zscf
[2]	1356	REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
	1357	REAL z2geomf, zalh2, zalm2, zscfh, zscfm
	1358	REAL t_coup
	1359	PARAMETER (t_coup=273.15)
[467]	1360	cIM
	1361	LOGICAL check
	1362	PARAMETER (check=.false.)
[2]	1363	c
	1364	c contre-gradient pour la chaleur sensible: Kelvin/metre
	1365	REAL gamt(2:klev)
[438]	1366	real qsurf(klon)
[2]	1367	c
	1368	LOGICAL appel1er
	1369	SAVE appel1er
	1370	c
	1371	c Fonctions thermodynamiques et fonctions d'instabilite
	1372	REAL fsta, fins, x
	1373	LOGICAL zxli ! utiliser un jeu de fonctions simples
	1374	PARAMETER (zxli=.FALSE.)
	1375	c
	1376	#include "YOETHF.h"
	1377	#include "FCTTRE.h"
	1378	fsta(x) = 1.0 / (1.0+10.0x(1+8.0*x))
	1379	fins(x) = SQRT(1.0-18.0*x)
	1380	c
	1381	DATA appel1er /.TRUE./
	1382	c
	1383	IF (appel1er) THEN
	1384	PRINT*, 'coefkz, opt_ec:', opt_ec
	1385	PRINT*, 'coefkz, richum:', richum
	1386	IF (richum) PRINT*, 'coefkz, ratqs:', ratqs
	1387	PRINT*, 'coefkz, isommet:', isommet
	1388	PRINT*, 'coefkz, tvirtu:', tvirtu
	1389	appel1er = .FALSE.
	1390	ENDIF
	1391	c
	1392	c Initialiser les sorties
	1393	c
	1394	DO k = 1, klev
	1395	DO i = 1, knon
	1396	pcfm(i,k) = 0.0
	1397	pcfh(i,k) = 0.0
	1398	ENDDO
	1399	ENDDO
	1400	DO i = 1, knon
	1401	itop(i) = 0
	1402	ENDDO
[280]	1403
[438]	1404	c$$$ IF(nsrf.NE.1) THEN
	1405	c$$$ do i = 1, knon
	1406	c$$$ qsurf(i) = qsatl(ts(i))/paprs(i,1)
	1407	c$$$ enddo
	1408	c$$$ ENDIF
[280]	1409
[2]	1410	c
	1411	c Prescrire la valeur de contre-gradient
	1412	c
[39]	1413	IF (.NOT.contreg) THEN
[2]	1414	DO k = 2, klev
	1415	gamt(k) = 0.0
	1416	ENDDO
	1417	ELSE
	1418	DO k = 3, klev
	1419	gamt(k) = -1.0E-03
	1420	ENDDO
	1421	gamt(2) = -2.5E-03
	1422	ENDIF
[486]	1423	cIM cf JLD/ GKtest
[489]	1424	IF ( nsrf .NE. is_oce ) THEN
[486]	1425	cIM 261103 kstable = kstable_ter
	1426	kstable = ksta_ter
	1427	ELSE
	1428	cIM 261103 kstable = kstable_sinon
	1429	kstable = ksta
	1430	ENDIF
	1431	cIM cf JLD/ GKtest fin
[2]	1432	c
	1433	c Calculer les geopotentiels de chaque couche
	1434	c
	1435	DO i = 1, knon
	1436	zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
	1437	. * (paprs(i,1)-pplay(i,1))
	1438	ENDDO
	1439	DO k = 2, klev
	1440	DO i = 1, knon
	1441	zgeop(i,k) = zgeop(i,k-1)
	1442	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
	1443	. * (pplay(i,k-1)-pplay(i,k))
	1444	ENDDO
	1445	ENDDO
	1446	c
	1447	c Calculer le frottement au sol (Cdrag)
	1448	c
[433]	1449	DO i = 1, knon
[407]	1450	u1(i) = u(i,1)
	1451	v1(i) = v(i,1)
	1452	t1(i) = t(i,1)
	1453	q1(i) = q(i,1)
	1454	z1(i) = zgeop(i,1)
	1455	ENDDO
[2]	1456	c
[407]	1457	CALL clcdrag(klon, knon, nsrf, zxli,
	1458	$ u1, v1, t1, q1, z1,
[448]	1459	$ ts, qsurf, rugos,
[407]	1460	$ pcfm1, pcfh1)
[433]	1461	cIM $ ts, qsurf, rugos,
[407]	1462	C
	1463	DO i = 1, knon
	1464	pcfm(i,1)=pcfm1(i)
	1465	pcfh(i,1)=pcfh1(i)
	1466	ENDDO
	1467	c
[2]	1468	c Calculer les coefficients turbulents dans l'atmosphere
	1469	c
	1470	DO i = 1, knon
	1471	itop(i) = isommet
	1472	ENDDO
	1473
[467]	1474	IF (check) THEN
	1475	PRINT*,' isommet=',isommet,' knon=',knon
	1476	ENDIF
[407]	1477
[2]	1478	DO k = 2, isommet
	1479	DO i = 1, knon
	1480	zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
	1481	. +(v(i,k)-v(i,k-1))**2)
	1482	zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
	1483	zdphi =zmgeom(i) / 2.0
	1484	zt = (t(i,k)+t(i,k-1)) * 0.5
	1485	zq = (q(i,k)+q(i,k-1)) * 0.5
	1486	c
	1487	c calculer Qs et dQs/dT:
	1488	c
	1489	IF (thermcep) THEN
	1490	zdelta = MAX(0.,SIGN(1.,RTT-zt))
	1491	zcvm5 = R5LESRLVTT/RCPD/(1.0+RVTMP2zq)*(1.-zdelta)
	1492	. + R5IESRLSTT/RCPD/(1.0+RVTMP2zq)*zdelta
	1493	zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
	1494	zqs = MIN(0.5,zqs)
	1495	zcor = 1./(1.-RETV*zqs)
	1496	zqs = zqs*zcor
	1497	zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
	1498	ELSE
	1499	IF (zt .LT. t_coup) THEN
	1500	zqs = qsats(zt) / pplay(i,k)
	1501	zdqs = dqsats(zt,zqs)
	1502	ELSE
	1503	zqs = qsatl(zt) / pplay(i,k)
	1504	zdqs = dqsatl(zt,zqs)
	1505	ENDIF
	1506	ENDIF
	1507	c
	1508	c calculer la fraction nuageuse (processus humide):
	1509	c
	1510	zfr = (zq+ratqszq-zqs) / (2.0ratqs*zq)
	1511	zfr = MAX(0.0,MIN(1.0,zfr))
	1512	IF (.NOT.richum) zfr = 0.0
	1513	c
	1514	c calculer le nombre de Richardson:
	1515	c
	1516	IF (tvirtu) THEN
	1517	ztvd =( t(i,k)
	1518	. + zdphi/RCPD/(1.+RVTMP2*zq)
	1519	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
	1520	. )(1.+RETVq(i,k))
	1521	ztvu =( t(i,k-1)
	1522	. - zdphi/RCPD/(1.+RVTMP2*zq)
	1523	. ( (1.-zfr) + zfr(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
	1524	. )(1.+RETVq(i,k-1))
	1525	zri(i) =zmgeom(i)(ztvd-ztvu)/(zdu20.5*(ztvd+ztvu))
	1526	zri(i) = zri(i)
	1527	. + zmgeom(i)zmgeom(i)/RGgamt(k)
	1528	. (paprs(i,k)/101325.0)*RKAPPA
	1529	. /(zdu20.5(ztvd+ztvu))
	1530	c
	1531	ELSE ! calcul de Ridchardson compatible LMD5
	1532	c
	1533	zri(i) =(RCPD*(t(i,k)-t(i,k-1))
	1534	. -RD0.5(t(i,k)+t(i,k-1))/paprs(i,k)
	1535	. *(pplay(i,k)-pplay(i,k-1))
	1536	. )zmgeom(i)/(zdu20.5RCPD(t(i,k-1)+t(i,k)))
	1537	zri(i) = zri(i) +
	1538	. zmgeom(i)zmgeom(i)gamt(k)/RG
	1539	cSB . /(paprs(i,k)/101325.0)**RKAPPA
	1540	. (paprs(i,k)/101325.0)*RKAPPA
	1541	. /(zdu20.5(t(i,k-1)+t(i,k)))
	1542	ENDIF
	1543	c
	1544	c finalement, les coefficients d'echange sont obtenus:
	1545	c
	1546	zcdn=SQRT(zdu2) / zmgeom(i) * RG
	1547	c
	1548	IF (opt_ec) THEN
	1549	z2geomf=zgeop(i,k-1)+zgeop(i,k)
	1550	zalm2=(0.5ckap/RGz2geomf
	1551	. /(1.+0.5ckap/rg/clamz2geomf))**2
	1552	zalh2=(0.5ckap/rgz2geomf
	1553	. /(1.+0.5ckap/RG/(clamSQRT(1.5cd))z2geomf))**2
	1554	IF (zri(i).LT.0.0) THEN ! situation instable
	1555	zscf = ((zgeop(i,k)/zgeop(i,k-1))(1./3.)-1.)3
	1556	. / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
	1557	zscf = SQRT(-zri(i)*zscf)
	1558	zscfm = 1.0 / (1.0+3.0cbcczalm2zscf)
	1559	zscfh = 1.0 / (1.0+3.0cbcczalh2zscf)
	1560	pcfm(i,k)=zcdnzalm2(1.-2.0cbzri(i)*zscfm)
	1561	pcfh(i,k)=zcdnzalh2(1.-3.0cbzri(i)*zscfh)
	1562	ELSE ! situation stable
	1563	zscf=SQRT(1.+cd*zri(i))
	1564	pcfm(i,k)=zcdnzalm2/(1.+2.0cb*zri(i)/zscf)
	1565	pcfh(i,k)=zcdnzalh2/(1.+3.0cbzri(i)zscf)
	1566	ENDIF
	1567	ELSE
	1568	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
	1569	. /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
	1570	pcfm(i,k)=sqrt(max(zcdnzcdn(ric-zri(i))/ric, kstable))
	1571	pcfm(i,k)= zl2(i)* pcfm(i,k)
	1572	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	1573	ENDIF
	1574	ENDDO
	1575	ENDDO
	1576	c
	1577	c Au-dela du sommet, pas de diffusion turbulente:
	1578	c
	1579	DO i = 1, knon
	1580	IF (itop(i)+1 .LE. klev) THEN
	1581	DO k = itop(i)+1, klev
	1582	pcfh(i,k) = 0.0
	1583	pcfm(i,k) = 0.0
	1584	ENDDO
	1585	ENDIF
	1586	ENDDO
	1587	c
	1588	RETURN
	1589	END
[336]	1590
[39]	1591	SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
	1592	. pcfm, pcfh)
	1593	IMPLICIT none
	1594	c======================================================================
	1595	c J'introduit un peu de diffusion sauf dans les endroits
	1596	c ou une forte inversion est presente
	1597	c On peut dire qu'il represente la convection peu profonde
	1598	c
	1599	c Arguments:
	1600	c nsrf-----input-I- indicateur de la nature du sol
	1601	c knon-----input-I- nombre de points a traiter
	1602	c paprs----input-R- pression a chaque intercouche (en Pa)
	1603	c pplay----input-R- pression au milieu de chaque couche (en Pa)
	1604	c t--------input-R- temperature (K)
	1605	c
	1606	c pcfm-----output-R- coefficients a calculer (vitesse)
	1607	c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
	1608	c======================================================================
	1609	#include "dimensions.h"
	1610	#include "dimphy.h"
	1611	#include "YOMCST.h"
	1612	#include "indicesol.h"
	1613	c
	1614	c Arguments:
	1615	c
	1616	INTEGER knon, nsrf
[148]	1617	REAL paprs(klon,klev+1), pplay(klon,klev)
	1618	REAL t(klon,klev)
[39]	1619	c
[148]	1620	REAL pcfm(klon,klev), pcfh(klon,klev)
[39]	1621	c
	1622	c Quelques constantes et options:
	1623	c
	1624	REAL prandtl
	1625	PARAMETER (prandtl=0.4)
	1626	REAL kstable
	1627	PARAMETER (kstable=0.002)
	1628	ccc PARAMETER (kstable=0.001)
	1629	REAL mixlen ! constante controlant longueur de melange
	1630	PARAMETER (mixlen=35.0)
	1631	REAL seuil ! au-dela l'inversion est consideree trop faible
	1632	PARAMETER (seuil=-0.02)
	1633	ccc PARAMETER (seuil=-0.04)
	1634	ccc PARAMETER (seuil=-0.06)
	1635	ccc PARAMETER (seuil=-0.09)
	1636	c
	1637	c Variables locales:
	1638	c
[146]	1639	INTEGER i, k, invb(knon)
	1640	REAL zl2(knon)
	1641	REAL zdthmin(knon), zdthdp
[39]	1642	c
	1643	c Initialiser les sorties
	1644	c
	1645	DO k = 1, klev
	1646	DO i = 1, knon
	1647	pcfm(i,k) = 0.0
	1648	pcfh(i,k) = 0.0
	1649	ENDDO
	1650	ENDDO
	1651	c
	1652	c Chercher la zone d'inversion forte
	1653	c
[146]	1654	DO i = 1, knon
[39]	1655	invb(i) = klev
	1656	zdthmin(i)=0.0
	1657	ENDDO
	1658	DO k = 2, klev/2-1
[146]	1659	DO i = 1, knon
[39]	1660	zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
	1661	. - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
	1662	zdthdp = zdthdp * 100.0
	1663	IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
	1664	. zdthdp.LT.zdthmin(i) ) THEN
	1665	zdthmin(i) = zdthdp
	1666	invb(i) = k
	1667	ENDIF
	1668	ENDDO
	1669	ENDDO
	1670	c
	1671	c Introduire une diffusion:
	1672	c
	1673	DO k = 2, klev
	1674	DO i = 1, knon
[486]	1675	cIM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean
	1676	cIM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
	1677	!IM cf JLD/ GKtest TERkz2
	1678	! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre
	1679	! fin GKtest
	1680	IF ( (nsrf.EQ.is_oce) .AND. ! si on est sur ocean et si
	1681	. ( (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
	1682	. (zdthmin(i).GT.seuil) ) )THEN ! si l'inversion est trop faible
[39]	1683	zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
	1684	. /(paprs(i,2)-paprs(i,klev+1)) ))**2
	1685	pcfm(i,k)= zl2(i)* kstable
	1686	pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
	1687	ENDIF
	1688	ENDDO
	1689	ENDDO
	1690	c
	1691	RETURN
	1692	END
[98]	1693	SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
[2]	1694	. vbeta,vcal,vdif)
	1695	IMPLICIT none
	1696	c======================================================================
	1697	c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
	1698	c date: 19940414
	1699	c======================================================================
	1700	c
	1701	c Calculer quelques parametres pour appliquer la couche limite
	1702	c ------------------------------------------------------------
[156]	1703	#include "dimensions.h"
	1704	#include "dimphy.h"
[2]	1705	#include "YOMCST.h"
	1706	#include "indicesol.h"
	1707	REAL tau_gl ! temps de relaxation pour la glace de mer
[39]	1708	ccc PARAMETER (tau_gl=86400.0*30.0)
	1709	PARAMETER (tau_gl=86400.0*5.0)
[2]	1710	REAL mx_eau_sol
	1711	PARAMETER (mx_eau_sol=150.0)
	1712	c
	1713	REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
	1714	PARAMETER (calsol=1.0/(2.5578E+06*0.15))
	1715	PARAMETER (calsno=1.0/(2.3867E+06*0.15))
	1716	PARAMETER (calice=1.0/(5.1444E+06*0.15))
	1717	C
	1718	INTEGER i
	1719	c
	1720	REAL dtime
[156]	1721	REAL snow(klon), qsol(klon)
[98]	1722	INTEGER indice, knon
[2]	1723	C
[156]	1724	REAL vbeta(klon)
	1725	REAL vcal(klon)
	1726	REAL vdif(klon)
[2]	1727	C
[156]	1728
[2]	1729	IF (indice.EQ.is_oce) THEN
[98]	1730	DO i = 1, knon
[2]	1731	vcal(i) = 0.0
	1732	vbeta(i) = 1.0
	1733	vdif(i) = 0.0
	1734	ENDDO
	1735	ENDIF
	1736	c
	1737	IF (indice.EQ.is_sic) THEN
[98]	1738	DO i = 1, knon
[2]	1739	vcal(i) = calice
[156]	1740	IF (snow(i) .GT. 0.0) vcal(i) = calsno
[2]	1741	vbeta(i) = 1.0
[39]	1742	vdif(i) = 1.0/tau_gl
	1743	ccc vdif(i) = calice/tau_gl ! c'etait une erreur
[2]	1744	ENDDO
	1745	ENDIF
	1746	c
	1747	IF (indice.EQ.is_ter) THEN
[98]	1748	DO i = 1, knon
[2]	1749	vcal(i) = calsol
[156]	1750	IF (snow(i) .GT. 0.0) vcal(i) = calsno
	1751	vbeta(i) = MIN(2.0*qsol(i)/mx_eau_sol, 1.0)
[2]	1752	vdif(i) = 0.0
	1753	ENDDO
	1754	ENDIF
	1755	c
	1756	IF (indice.EQ.is_lic) THEN
[98]	1757	DO i = 1, knon
[2]	1758	vcal(i) = calice
[156]	1759	IF (snow(i) .GT. 0.0) vcal(i) = calsno
[2]	1760	vbeta(i) = 1.0
	1761	vdif(i) = 0.0
	1762	ENDDO
	1763	ENDIF
	1764	c
	1765	RETURN
	1766	END
	1767	C======================================================================
	1768	SUBROUTINE nonlocal(knon, paprs, pplay,
	1769	. tsol,beta,u,v,t,q,
	1770	. cd_h, cd_m, pcfh, pcfm, cgh, cgq)
	1771	IMPLICIT none
	1772	c======================================================================
	1773	c Laurent Li (LMD/CNRS), le 30 septembre 1998
	1774	c Couche limite non-locale. Adaptation du code du CCM3.
	1775	c Code non teste, donc a ne pas utiliser.
	1776	c======================================================================
	1777	c Nonlocal scheme that determines eddy diffusivities based on a
	1778	c diagnosed boundary layer height and a turbulent velocity scale.
	1779	c Also countergradient effects for heat and moisture are included.
	1780	c
	1781	c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
	1782	c Local versus nonlocal boundary-layer diffusion in a global climate
	1783	c model. J. of Climate, vol. 6, 1825-1842.
	1784	c======================================================================
	1785	#include "dimensions.h"
	1786	#include "dimphy.h"
	1787	#include "YOMCST.h"
	1788	c
	1789	c Arguments:
	1790	c
	1791	INTEGER knon ! nombre de points a calculer
	1792	REAL tsol(klon) ! temperature du sol (K)
	1793	REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
	1794	REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
	1795	REAL pplay(klon,klev) ! pression au milieu de couche (Pa)
	1796	REAL u(klon,klev) ! vitesse U (m/s)
	1797	REAL v(klon,klev) ! vitesse V (m/s)
	1798	REAL t(klon,klev) ! temperature (K)
	1799	REAL q(klon,klev) ! vapeur d'eau (kg/kg)
	1800	REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
	1801	REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
	1802	c
	1803	INTEGER isommet
	1804	PARAMETER (isommet=klev)
	1805	REAL vk
[407]	1806	PARAMETER (vk=0.40)
[2]	1807	REAL ricr
	1808	PARAMETER (ricr=0.4)
	1809	REAL fak
	1810	PARAMETER (fak=8.5)
	1811	REAL fakn
	1812	PARAMETER (fakn=7.2)
	1813	REAL onet
	1814	PARAMETER (onet=1.0/3.0)
	1815	REAL t_coup
	1816	PARAMETER(t_coup=273.15)
	1817	REAL zkmin
	1818	PARAMETER (zkmin=0.01)
	1819	REAL betam
	1820	PARAMETER (betam=15.0)
	1821	REAL betah
	1822	PARAMETER (betah=15.0)
	1823	REAL betas
	1824	PARAMETER (betas=5.0)
	1825	REAL sffrac
	1826	PARAMETER (sffrac=0.1)
	1827	REAL binm
	1828	PARAMETER (binm=betam*sffrac)
	1829	REAL binh
	1830	PARAMETER (binh=betah*sffrac)
	1831	REAL ccon
	1832	PARAMETER (ccon=faksffracvk)
	1833	c
	1834	REAL z(klon,klev)
	1835	REAL pcfm(klon,klev), pcfh(klon,klev)
	1836	c
	1837	INTEGER i, k
	1838	REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
	1839	REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
	1840	REAL khfs(klon) ! surface kinematic heat flux [mK/s]
	1841	REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
	1842	REAL heatv(klon) ! surface virtual heat flux
	1843	REAL ustar(klon)
	1844	REAL rino(klon,klev) ! bulk Richardon no. from level to ref lev
	1845	LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
	1846	LOGICAL stblev(klon) ! stable pbl with levels within pbl
	1847	LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
	1848	LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
	1849	LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
	1850	LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
	1851	REAL pblh(klon)
	1852	REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
	1853	REAL cgq(klon,2:klev) ! counter-gradient term for constituents
	1854	REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
	1855	REAL obklen(klon)
	1856	REAL ztvd, ztvu, zdu2
	1857	REAL therm(klon) ! thermal virtual temperature excess
	1858	REAL phiminv(klon) ! inverse phi function for momentum
	1859	REAL phihinv(klon) ! inverse phi function for heat
	1860	REAL wm(klon) ! turbulent velocity scale for momentum
	1861	REAL fak1(klon) ! kustarpblh
	1862	REAL fak2(klon) ! kwmpblh
	1863	REAL fak3(klon) ! fakn*wstr/wm
	1864	REAL pblk(klon) ! level eddy diffusivity for momentum
	1865	REAL pr(klon) ! Prandtl number for eddy diffusivities
	1866	REAL zl(klon) ! zmzp / Obukhov length
	1867	REAL zh(klon) ! zmzp / pblh
	1868	REAL zzh(klon) ! (1-(zmzp/pblh))**2
	1869	REAL wstr(klon) ! w*, convective velocity scale
	1870	REAL zm(klon) ! current level height
	1871	REAL zp(klon) ! current level height + one level up
	1872	REAL zcor, zdelta, zcvm5, zxqs
	1873	REAL fac, pblmin, zmzp, term
	1874	c
	1875	#include "YOETHF.h"
	1876	#include "FCTTRE.h"
	1877	c
	1878	c Initialisation
	1879	c
	1880	DO i = 1, klon
	1881	pcfh(i,1) = cd_h(i)
	1882	pcfm(i,1) = cd_m(i)
	1883	ENDDO
	1884	DO k = 2, klev
	1885	DO i = 1, klon
	1886	pcfh(i,k) = zkmin
	1887	pcfm(i,k) = zkmin
	1888	cgs(i,k) = 0.0
	1889	cgh(i,k) = 0.0
	1890	cgq(i,k) = 0.0
	1891	ENDDO
	1892	ENDDO
	1893	c
	1894	c Calculer les hauteurs de chaque couche
	1895	c
	1896	DO i = 1, knon
	1897	z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
	1898	. * (paprs(i,1)-pplay(i,1)) / RG
	1899	ENDDO
	1900	DO k = 2, klev
	1901	DO i = 1, knon
	1902	z(i,k) = z(i,k-1)
	1903	. + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
	1904	. * (pplay(i,k-1)-pplay(i,k)) / RG
	1905	ENDDO
	1906	ENDDO
	1907	c
	1908	DO i = 1, knon
	1909	IF (thermcep) THEN
	1910	zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
	1911	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	1912	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
	1913	zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
	1914	zxqs=MIN(0.5,zxqs)
	1915	zcor=1./(1.-retv*zxqs)
	1916	zxqs=zxqs*zcor
	1917	ELSE
	1918	IF (tsol(i).LT.t_coup) THEN
	1919	zxqs = qsats(tsol(i)) / paprs(i,1)
	1920	ELSE
	1921	zxqs = qsatl(tsol(i)) / paprs(i,1)
	1922	ENDIF
	1923	ENDIF
	1924	zx_alf1 = 1.0
	1925	zx_alf2 = 1.0 - zx_alf1
	1926	zxt = (t(i,1)+z(i,1)RG/RCPD/(1.+RVTMP2q(i,1)))
	1927	. (1.+RETVq(i,1))*zx_alf1
	1928	. + (t(i,2)+z(i,2)RG/RCPD/(1.+RVTMP2q(i,2)))
	1929	. (1.+RETVq(i,2))*zx_alf2
	1930	zxu = u(i,1)zx_alf1+u(i,2)zx_alf2
	1931	zxv = v(i,1)zx_alf1+v(i,2)zx_alf2
	1932	zxq = q(i,1)zx_alf1+q(i,2)zx_alf2
	1933	zxmod = 1.0+SQRT(zxu2+zxv2)
	1934	khfs(i) = (tsol(i)(1.+RETVq(i,1))-zxt) zxmodcd_h(i)
	1935	kqfs(i) = (zxqs-zxq) zxmodcd_h(i) * beta(i)
	1936	heatv(i) = khfs(i) + 0.61zxtkqfs(i)
	1937	taux = zxu zxmodcd_m(i)
	1938	tauy = zxv zxmodcd_m(i)
	1939	ustar(i) = SQRT(taux2+tauy2)
	1940	ustar(i) = MAX(SQRT(ustar(i)),0.01)
	1941	ENDDO
	1942	c
	1943	DO i = 1, knon
	1944	rino(i,1) = 0.0
	1945	check(i) = .TRUE.
	1946	pblh(i) = z(i,1)
	1947	obklen(i) = -t(i,1)ustar(i)3/(RGvk*heatv(i))
	1948	ENDDO
	1949
	1950	C
	1951	C PBL height calculation:
	1952	C Search for level of pbl. Scan upward until the Richardson number between
	1953	C the first level and the current level exceeds the "critical" value.
	1954	C
	1955	fac = 100.0
	1956	DO k = 1, isommet
	1957	DO i = 1, knon
	1958	IF (check(i)) THEN
	1959	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
	1960	zdu2 = max(zdu2,1.0e-20)
	1961	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
	1962	. (1.+RETVq(i,k))
	1963	ztvu =(t(i,1)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
	1964	. (1.+RETVq(i,1))
	1965	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
	1966	. /(zdu20.5(ztvd+ztvu))
	1967	IF (rino(i,k).GE.ricr) THEN
	1968	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
	1969	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
	1970	check(i) = .FALSE.
	1971	ENDIF
	1972	ENDIF
	1973	ENDDO
	1974	ENDDO
	1975
	1976	C
	1977	C Set pbl height to maximum value where computation exceeds number of
	1978	C layers allowed
	1979	C
	1980	DO i = 1, knon
	1981	if (check(i)) pblh(i) = z(i,isommet)
	1982	ENDDO
	1983	C
	1984	C Improve estimate of pbl height for the unstable points.
	1985	C Find unstable points (sensible heat flux is upward):
	1986	C
	1987	DO i = 1, knon
	1988	IF (heatv(i) .GT. 0.) THEN
	1989	unstbl(i) = .TRUE.
	1990	check(i) = .TRUE.
	1991	ELSE
	1992	unstbl(i) = .FALSE.
	1993	check(i) = .FALSE.
	1994	ENDIF
	1995	ENDDO
	1996	C
	1997	C For the unstable case, compute velocity scale and the
	1998	C convective temperature excess:
	1999	C
	2000	DO i = 1, knon
	2001	IF (check(i)) THEN
	2002	phiminv(i) = (1.-binmpblh(i)/obklen(i))*onet
	2003	wm(i)= ustar(i)*phiminv(i)
	2004	therm(i) = heatv(i)*fak/wm(i)
	2005	rino(i,1) = 0.0
	2006	ENDIF
	2007	ENDDO
	2008	C
	2009	C Improve pblh estimate for unstable conditions using the
	2010	C convective temperature excess:
	2011	C
	2012	DO k = 1, isommet
	2013	DO i = 1, knon
	2014	IF (check(i)) THEN
	2015	zdu2 = (u(i,k)-u(i,1))2+(v(i,k)-v(i,1))2+facustar(i)*2
	2016	zdu2 = max(zdu2,1.0e-20)
	2017	ztvd =(t(i,k)+z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,k)))
	2018	. (1.+RETVq(i,k))
	2019	ztvu =(t(i,1)+therm(i)-z(i,k)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
	2020	. (1.+RETVq(i,1))
	2021	rino(i,k) = (z(i,k)-z(i,1))RG(ztvd-ztvu)
	2022	. /(zdu20.5(ztvd+ztvu))
	2023	IF (rino(i,k).GE.ricr) THEN
	2024	pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
	2025	. (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
	2026	check(i) = .FALSE.
	2027	ENDIF
	2028	ENDIF
	2029	ENDDO
	2030	ENDDO
	2031	C
	2032	C Set pbl height to maximum value where computation exceeds number of
	2033	C layers allowed
	2034	C
	2035	DO i = 1, knon
	2036	if (check(i)) pblh(i) = z(i,isommet)
	2037	ENDDO
	2038	C
	2039	C Points for which pblh exceeds number of pbl layers allowed;
	2040	C set to maximum
	2041	C
	2042	DO i = 1, knon
	2043	IF (check(i)) pblh(i) = z(i,isommet)
	2044	ENDDO
	2045	C
	2046	C PBL height must be greater than some minimum mechanical mixing depth
	2047	C Several investigators have proposed minimum mechanical mixing depth
	2048	C relationships as a function of the local friction velocity, u*. We
	2049	C make use of a linear relationship of the form h = c u* where c=700.
	2050	C The scaling arguments that give rise to this relationship most often
	2051	C represent the coefficient c as some constant over the local coriolis
	2052	C parameter. Here we make use of the experimental results of Koracin
	2053	C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
	2054	C where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
	2055	C latitude value for f so that c = 0.07/f = 700.
	2056	C
	2057	DO i = 1, knon
	2058	pblmin = 700.0*ustar(i)
	2059	pblh(i) = MAX(pblh(i),pblmin)
	2060	ENDDO
	2061	C
	2062	C pblh is now available; do preparation for diffusivity calculation:
	2063	C
	2064	DO i = 1, knon
	2065	pblk(i) = 0.0
	2066	fak1(i) = ustar(i)pblh(i)vk
	2067	C
	2068	C Do additional preparation for unstable cases only, set temperature
	2069	C and moisture perturbations depending on stability.
	2070	C
	2071	IF (unstbl(i)) THEN
	2072	zxt=(t(i,1)-z(i,1)0.5RG/RCPD/(1.+RVTMP2*q(i,1)))
	2073	. (1.+RETVq(i,1))
	2074	phiminv(i) = (1. - binmpblh(i)/obklen(i))*onet
	2075	phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
	2076	wm(i) = ustar(i)*phiminv(i)
	2077	fak2(i) = wm(i)pblh(i)vk
	2078	wstr(i) = (heatv(i)RGpblh(i)/zxt)**onet
	2079	fak3(i) = fakn*wstr(i)/wm(i)
	2080	ENDIF
	2081	ENDDO
	2082
	2083	C Main level loop to compute the diffusivities and
	2084	C counter-gradient terms:
	2085	C
	2086	DO 1000 k = 2, isommet
	2087	C
	2088	C Find levels within boundary layer:
	2089	C
	2090	DO i = 1, knon
	2091	unslev(i) = .FALSE.
	2092	stblev(i) = .FALSE.
	2093	zm(i) = z(i,k-1)
	2094	zp(i) = z(i,k)
	2095	IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
	2096	IF (zm(i) .LT. pblh(i)) THEN
	2097	zmzp = 0.5*(zm(i) + zp(i))
	2098	zh(i) = zmzp/pblh(i)
	2099	zl(i) = zmzp/obklen(i)
	2100	zzh(i) = 0.
	2101	IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
	2102	C
	2103	C stblev for points zm < plbh and stable and neutral
	2104	C unslev for points zm < plbh and unstable
	2105	C
	2106	IF (unstbl(i)) THEN
	2107	unslev(i) = .TRUE.
	2108	ELSE
	2109	stblev(i) = .TRUE.
	2110	ENDIF
	2111	ENDIF
	2112	ENDDO
	2113	C
	2114	C Stable and neutral points; set diffusivities; counter-gradient
	2115	C terms zero for stable case:
	2116	C
	2117	DO i = 1, knon
	2118	IF (stblev(i)) THEN
	2119	IF (zl(i).LE.1.) THEN
	2120	pblk(i) = fak1(i)zh(i)zzh(i)/(1. + betas*zl(i))
	2121	ELSE
	2122	pblk(i) = fak1(i)zh(i)zzh(i)/(betas + zl(i))
	2123	ENDIF
	2124	pcfm(i,k) = pblk(i)
	2125	pcfh(i,k) = pcfm(i,k)
	2126	ENDIF
	2127	ENDDO
	2128	C
	2129	C unssrf, unstable within surface layer of pbl
	2130	C unsout, unstable within outer layer of pbl
	2131	C
	2132	DO i = 1, knon
	2133	unssrf(i) = .FALSE.
	2134	unsout(i) = .FALSE.
	2135	IF (unslev(i)) THEN
	2136	IF (zh(i).lt.sffrac) THEN
	2137	unssrf(i) = .TRUE.
	2138	ELSE
	2139	unsout(i) = .TRUE.
	2140	ENDIF
	2141	ENDIF
	2142	ENDDO
	2143	C
	2144	C Unstable for surface layer; counter-gradient terms zero
	2145	C
	2146	DO i = 1, knon
	2147	IF (unssrf(i)) THEN
	2148	term = (1. - betamzl(i))*onet
	2149	pblk(i) = fak1(i)zh(i)zzh(i)*term
	2150	pr(i) = term/sqrt(1. - betah*zl(i))
	2151	ENDIF
	2152	ENDDO
	2153	C
	2154	C Unstable for outer layer; counter-gradient terms non-zero:
	2155	C
	2156	DO i = 1, knon
	2157	IF (unsout(i)) THEN
	2158	pblk(i) = fak2(i)zh(i)zzh(i)
	2159	cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
	2160	cgh(i,k) = khfs(i)*cgs(i,k)
	2161	pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
	2162	cgq(i,k) = kqfs(i)*cgs(i,k)
	2163	ENDIF
	2164	ENDDO
	2165	C
	2166	C For all unstable layers, set diffusivities
	2167	C
	2168	DO i = 1, knon
	2169	IF (unslev(i)) THEN
	2170	pcfm(i,k) = pblk(i)
	2171	pcfh(i,k) = pblk(i)/pr(i)
	2172	ENDIF
	2173	ENDDO
	2174	1000 continue ! end of level loop
	2175
	2176	RETURN
	2177	END

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

Download in other formats:

Original Format