Context Navigation

clmain.F @ 1962

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

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

Download in other formats:

Original Format