Context Navigation

clmain.F @ 232

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

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

Download in other formats:

Original Format