Context Navigation

clmain.F @ 3698

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

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

Download in other formats:

Original Format