Context Navigation

clmain.F @ 525

Last change on this file since 525 was 524, checked in by lmdzadmin, 21 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 68.9 KB

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

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

Download in other formats:

Original Format