Context Navigation

clmain.F @ 634

Last change on this file since 634 was 634, checked in by Laurent Fairhead, 19 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 72.7 KB

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

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

Download in other formats:

Original Format