Context Navigation

clmain.F @ 2471

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

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

Download in other formats:

Original Format