Context Navigation

clmain.F @ 407

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

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

Download in other formats:

Original Format