Context Navigation

clmain.F @ 633

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

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

Download in other formats:

Original Format