Context Navigation

clmain.F @ 180

Last change on this file since 180 was 177, checked in by lmdzadmin, 23 years ago

Lots of stuff, plus particulierement:

appel a ORCHIDEE en etat de marche (pb de grille subsiste)
modifs de Pascale sur soil dans le cas ou ok_veget=false -

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

File size: 58.5 KB

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

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

Download in other formats:

Original Format