Context Navigation

clmain.F @ 5420

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

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

Download in other formats:

Original Format