Context Navigation

clmain.F @ 710

Last change on this file since 710 was 710, checked in by Laurent Fairhead, 18 years ago
Integration des modifications d'OM et JLD pour une meilleure prise en compte de la conservation du flux d'eau LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 69.2 KB

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

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

Download in other formats:

Original Format