Context Navigation

clmain.F @ 705

Last change on this file since 705 was 705, checked in by Laurent Fairhead, 18 years ago
Pour la conservation du flux d'eau OM, JLD LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 69.7 KB

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

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

Download in other formats:

Original Format