Context Navigation

clmain.F @ 222

Last change on this file since 222 was 205, checked in by lmdzadmin, 23 years ago
Rajout date0 dans physiq pour interface avec ORCHIDEE Fermeture des fichiers restart par restclo dans abort_gcm LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 58.6 KB

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

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

Download in other formats:

Original Format