Context Navigation

clmain.F @ 281

Last change on this file since 281 was 280, checked in by lmdzadmin, 23 years ago
Modif dans coefkz pour remplacer l'humidite du niveau 1 par qsat(tsol) Passage des deux albedos de surface (vis et nir) LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 61.6 KB

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

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

Download in other formats:

Original Format