Context Navigation

clmain.F @ 545

Last change on this file since 545 was 522, checked in by lmdzadmin, 21 years ago
Modifications pour la fermeture en eau (fonte des glaciers) JLD, OM LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 68.9 KB

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

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

Download in other formats:

Original Format