Context Navigation

clmain.F @ 361

Last change on this file since 361 was 353, checked in by lmdzadmin, 23 years ago

2 changements pour les fichiers histoire:

utilisation de l'entree "rectilineaire" de IOIPSL pour ne plus avoir

lancer ncregular a chaque fois

le calendrier des fichiers histoire est maintenant base sur la date d'initialisation de la simulation plutot que sur la date de depart du

job

en cours

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 63.6 KB

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

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

Download in other formats:

Original Format