Context Navigation

clmain.F @ 5373

Last change on this file since 5373 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 72.7 KB

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

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

Download in other formats:

Original Format