Context Navigation

clmain.F @ 342

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

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

Download in other formats:

Original Format