Context Navigation

clmain.F @ 290

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

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

Download in other formats:

Original Format