Context Navigation

clmain.F @ 589

Last change on this file since 589 was 541, checked in by lmdzadmin, 20 years ago

Convergence avec la version d'Olivia Coindreau incluant:

le offline
les thermiques
mellor & yamada dans la couche limite

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format