Context Navigation

clmain.F @ 258

Last change on this file since 258 was 258, checked in by lmdzadmin, 23 years ago
Phasage avec la version de PB pour le sol, dlw (juillet 2001) LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 61.1 KB

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

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

Download in other formats:

Original Format