Context Navigation

clmain.F @ 558

Last change on this file since 558 was 546, checked in by lmdzadmin, 20 years ago
Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory) LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 69.0 KB

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

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

Download in other formats:

Original Format