Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_u.f90 @ 1227

Last change on this file since 1227 was 1132, checked in by raasch, 12 years ago
r1028 documented
Property svn:keywords set to `Id`
File size: 24.4 KB

Line
1	MODULE diffusion_u_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_u.f90 1132 2013-04-12 14:35:30Z fuhrmann $
27	!
28	! 1128 2013-04-12 06:19:32Z raasch
29	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
30	! j_north
31	!
32	! 1036 2012-10-22 13:43:42Z raasch
33	! code put under GPL (PALM 3.9)
34	!
35	! 1015 2012-09-27 09:23:24Z raasch
36	! accelerator version (*_acc) added
37	!
38	! 1001 2012-09-13 14:08:46Z raasch
39	! arrays comunicated by module instead of parameter list
40	!
41	! 978 2012-08-09 08:28:32Z fricke
42	! outflow damping layer removed
43	! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
44	!
45	! 667 2010-12-23 12:06:00Z suehring/gryschka
46	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
47	!
48	! 366 2009-08-25 08:06:27Z raasch
49	! bc_ns replaced by bc_ns_cyc
50	!
51	! 106 2007-08-16 14:30:26Z raasch
52	! Momentumflux at top (uswst) included as boundary condition,
53	! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
54	!
55	! 75 2007-03-22 09:54:05Z raasch
56	! Wall functions now include diabatic conditions, call of routine wall_fluxes,
57	! z0 removed from argument list, uxrp eliminated
58	!
59	! 20 2007-02-26 00:12:32Z raasch
60	! Bugfix: ddzw dimensioned 1:nzt"+1"
61	!
62	! RCS Log replace by Id keyword, revision history cleaned up
63	!
64	! Revision 1.15 2006/02/23 10:35:35 raasch
65	! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
66	! or nzb_diff_u, respectively, in vertical diffusion,
67	! wall functions added for north and south walls, +z0 in argument list,
68	! terms containing w(k-1,..) are removed from the Prandtl-layer equation
69	! because they cause errors at the edges of topography
70	! WARNING: loops containing the MAX function are still not properly vectorized!
71	!
72	! Revision 1.1 1997/09/12 06:23:51 raasch
73	! Initial revision
74	!
75	!
76	! Description:
77	! ------------
78	! Diffusion term of the u-component
79	! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
80	! and slows down the speed on NEC about 5-10%
81	!------------------------------------------------------------------------------!
82
83	USE wall_fluxes_mod
84
85	PRIVATE
86	PUBLIC diffusion_u, diffusion_u_acc
87
88	INTERFACE diffusion_u
89	MODULE PROCEDURE diffusion_u
90	MODULE PROCEDURE diffusion_u_ij
91	END INTERFACE diffusion_u
92
93	INTERFACE diffusion_u_acc
94	MODULE PROCEDURE diffusion_u_acc
95	END INTERFACE diffusion_u_acc
96
97	CONTAINS
98
99
100	!------------------------------------------------------------------------------!
101	! Call for all grid points
102	!------------------------------------------------------------------------------!
103	SUBROUTINE diffusion_u
104
105	USE arrays_3d
106	USE control_parameters
107	USE grid_variables
108	USE indices
109
110	IMPLICIT NONE
111
112	INTEGER :: i, j, k
113	REAL :: kmym, kmyp, kmzm, kmzp
114
115	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
116
117	!
118	!-- First calculate horizontal momentum flux u'v' at vertical walls,
119	!-- if neccessary
120	IF ( topography /= 'flat' ) THEN
121	CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
122	nzb_u_outer, wall_u )
123	ENDIF
124
125	DO i = nxlu, nxr
126	DO j = nys, nyn
127	!
128	!-- Compute horizontal diffusion
129	DO k = nzb_u_outer(j,i)+1, nzt
130	!
131	!-- Interpolate eddy diffusivities on staggered gridpoints
132	kmyp = 0.25 * &
133	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
134	kmym = 0.25 * &
135	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
136
137	tend(k,j,i) = tend(k,j,i) &
138	& + 2.0 * ( &
139	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
140	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
141	& ) * ddx2 &
142	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
143	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
144	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
145	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
146	& ) * ddy
147	ENDDO
148
149	!
150	!-- Wall functions at the north and south walls, respectively
151	IF ( wall_u(j,i) /= 0.0 ) THEN
152
153	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
154	kmyp = 0.25 * &
155	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
156	kmym = 0.25 * &
157	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
158
159	tend(k,j,i) = tend(k,j,i) &
160	+ 2.0 * ( &
161	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
162	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
163	) * ddx2 &
164	+ ( fyp(j,i) * ( &
165	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
166	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
167	) &
168	- fym(j,i) * ( &
169	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
170	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
171	) &
172	+ wall_u(j,i) * usvs(k,j,i) &
173	) * ddy
174	ENDDO
175	ENDIF
176
177	!
178	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
179	!-- index k starts at nzb_u_inner+2.
180	DO k = nzb_diff_u(j,i), nzt_diff
181	!
182	!-- Interpolate eddy diffusivities on staggered gridpoints
183	kmzp = 0.25 * &
184	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
185	kmzm = 0.25 * &
186	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
187
188	tend(k,j,i) = tend(k,j,i) &
189	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
190	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
191	& ) &
192	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
193	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
194	& ) &
195	& ) * ddzw(k)
196	ENDDO
197
198	!
199	!-- Vertical diffusion at the first grid point above the surface,
200	!-- if the momentum flux at the bottom is given by the Prandtl law or
201	!-- if it is prescribed by the user.
202	!-- Difference quotient of the momentum flux is not formed over half
203	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
204	!-- with other (LES) modell showed that the values of the momentum
205	!-- flux becomes too large in this case.
206	!-- The term containing w(k-1,..) (see above equation) is removed here
207	!-- because the vertical velocity is assumed to be zero at the surface.
208	IF ( use_surface_fluxes ) THEN
209	k = nzb_u_inner(j,i)+1
210	!
211	!-- Interpolate eddy diffusivities on staggered gridpoints
212	kmzp = 0.25 * &
213	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
214	kmzm = 0.25 * &
215	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
216
217	tend(k,j,i) = tend(k,j,i) &
218	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
219	& ) * ddzw(k) &
220	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
221	& + usws(j,i) &
222	& ) * ddzw(k)
223	ENDIF
224
225	!
226	!-- Vertical diffusion at the first gridpoint below the top boundary,
227	!-- if the momentum flux at the top is prescribed by the user
228	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
229	k = nzt
230	!
231	!-- Interpolate eddy diffusivities on staggered gridpoints
232	kmzp = 0.25 * &
233	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
234	kmzm = 0.25 * &
235	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
236
237	tend(k,j,i) = tend(k,j,i) &
238	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
239	& ) * ddzw(k) &
240	& + ( -uswst(j,i) &
241	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
242	& ) * ddzw(k)
243	ENDIF
244
245	ENDDO
246	ENDDO
247
248	END SUBROUTINE diffusion_u
249
250
251	!------------------------------------------------------------------------------!
252	! Call for all grid points - accelerator version
253	!------------------------------------------------------------------------------!
254	SUBROUTINE diffusion_u_acc
255
256	USE arrays_3d
257	USE control_parameters
258	USE grid_variables
259	USE indices
260
261	IMPLICIT NONE
262
263	INTEGER :: i, j, k
264	REAL :: kmym, kmyp, kmzm, kmzp
265
266	!$acc declare create ( usvs )
267	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
268
269	!
270	!-- First calculate horizontal momentum flux u'v' at vertical walls,
271	!-- if neccessary
272	IF ( topography /= 'flat' ) THEN
273	CALL wall_fluxes_acc( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
274	nzb_u_outer, wall_u )
275	ENDIF
276
277	!$acc kernels present ( u, v, w, km, tend, usws, uswst ) &
278	!$acc present ( ddzu, ddzw, fym, fyp, wall_u ) &
279	!$acc present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
280	!$acc loop
281	DO i = i_left, i_right
282	DO j = j_south, j_north
283	!
284	!-- Compute horizontal diffusion
285	!$acc loop vector(32)
286	DO k = 1, nzt
287	IF ( k > nzb_u_outer(j,i) ) THEN
288	!
289	!-- Interpolate eddy diffusivities on staggered gridpoints
290	kmyp = 0.25 * &
291	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
292	kmym = 0.25 * &
293	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
294
295	tend(k,j,i) = tend(k,j,i) &
296	& + 2.0 * ( &
297	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
298	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
299	& ) * ddx2 &
300	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
301	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
302	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
303	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
304	& ) * ddy
305	ENDIF
306	ENDDO
307
308	!
309	!-- Wall functions at the north and south walls, respectively
310	!$acc loop vector(32)
311	DO k = 1, nzt
312	IF( k > nzb_u_inner(j,i) .AND. k <= nzb_u_outer(j,i) .AND. &
313	wall_u(j,i) /= 0.0 ) THEN
314
315	kmyp = 0.25 * &
316	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
317	kmym = 0.25 * &
318	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
319
320	tend(k,j,i) = tend(k,j,i) &
321	+ 2.0 * ( &
322	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
323	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
324	) * ddx2 &
325	+ ( fyp(j,i) * ( &
326	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
327	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
328	) &
329	- fym(j,i) * ( &
330	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
331	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
332	) &
333	+ wall_u(j,i) * usvs(k,j,i) &
334	) * ddy
335	ENDIF
336	ENDDO
337
338	!
339	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
340	!-- index k starts at nzb_u_inner+2.
341	!$acc loop vector(32)
342	DO k = 1, nzt_diff
343	IF ( k >= nzb_diff_u(j,i) ) THEN
344	!
345	!-- Interpolate eddy diffusivities on staggered gridpoints
346	kmzp = 0.25 * &
347	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
348	kmzm = 0.25 * &
349	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
350
351	tend(k,j,i) = tend(k,j,i) &
352	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&
353	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
354	& ) &
355	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)&
356	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
357	& ) &
358	& ) * ddzw(k)
359	ENDIF
360	ENDDO
361
362	ENDDO
363	ENDDO
364
365	!
366	!-- Vertical diffusion at the first grid point above the surface,
367	!-- if the momentum flux at the bottom is given by the Prandtl law or
368	!-- if it is prescribed by the user.
369	!-- Difference quotient of the momentum flux is not formed over half
370	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
371	!-- with other (LES) modell showed that the values of the momentum
372	!-- flux becomes too large in this case.
373	!-- The term containing w(k-1,..) (see above equation) is removed here
374	!-- because the vertical velocity is assumed to be zero at the surface.
375	IF ( use_surface_fluxes ) THEN
376
377	!$acc loop
378	DO i = i_left, i_right
379	!$acc loop vector(32)
380	DO j = j_south, j_north
381
382	k = nzb_u_inner(j,i)+1
383	!
384	!-- Interpolate eddy diffusivities on staggered gridpoints
385	kmzp = 0.25 * &
386	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
387	kmzm = 0.25 * &
388	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
389
390	tend(k,j,i) = tend(k,j,i) &
391	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
392	& ) * ddzw(k) &
393	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
394	& + usws(j,i) &
395	& ) * ddzw(k)
396	ENDDO
397	ENDDO
398
399	ENDIF
400
401	!
402	!-- Vertical diffusion at the first gridpoint below the top boundary,
403	!-- if the momentum flux at the top is prescribed by the user
404	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
405
406	k = nzt
407
408	!$acc loop
409	DO i = i_left, i_right
410	!$acc loop vector(32)
411	DO j = j_south, j_north
412
413	!
414	!-- Interpolate eddy diffusivities on staggered gridpoints
415	kmzp = 0.25 * &
416	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
417	kmzm = 0.25 * &
418	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
419
420	tend(k,j,i) = tend(k,j,i) &
421	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
422	& ) * ddzw(k) &
423	& + ( -uswst(j,i) &
424	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
425	& ) * ddzw(k)
426	ENDDO
427	ENDDO
428
429	ENDIF
430	!$acc end kernels
431
432	END SUBROUTINE diffusion_u_acc
433
434
435	!------------------------------------------------------------------------------!
436	! Call for grid point i,j
437	!------------------------------------------------------------------------------!
438	SUBROUTINE diffusion_u_ij( i, j )
439
440	USE arrays_3d
441	USE control_parameters
442	USE grid_variables
443	USE indices
444
445	IMPLICIT NONE
446
447	INTEGER :: i, j, k
448	REAL :: kmym, kmyp, kmzm, kmzp
449
450	REAL, DIMENSION(nzb:nzt+1) :: usvs
451
452	!
453	!-- Compute horizontal diffusion
454	DO k = nzb_u_outer(j,i)+1, nzt
455	!
456	!-- Interpolate eddy diffusivities on staggered gridpoints
457	kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
458	kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
459
460	tend(k,j,i) = tend(k,j,i) &
461	& + 2.0 * ( &
462	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
463	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
464	& ) * ddx2 &
465	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
466	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
467	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
468	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
469	& ) * ddy
470	ENDDO
471
472	!
473	!-- Wall functions at the north and south walls, respectively
474	IF ( wall_u(j,i) .NE. 0.0 ) THEN
475
476	!
477	!-- Calculate the horizontal momentum flux u'v'
478	CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i), &
479	usvs, 1.0, 0.0, 0.0, 0.0 )
480
481	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
482	kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
483	kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
484
485	tend(k,j,i) = tend(k,j,i) &
486	+ 2.0 * ( &
487	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
488	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
489	) * ddx2 &
490	+ ( fyp(j,i) * ( &
491	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
492	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
493	) &
494	- fym(j,i) * ( &
495	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
496	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
497	) &
498	+ wall_u(j,i) * usvs(k) &
499	) * ddy
500	ENDDO
501	ENDIF
502
503	!
504	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
505	!-- index k starts at nzb_u_inner+2.
506	DO k = nzb_diff_u(j,i), nzt_diff
507	!
508	!-- Interpolate eddy diffusivities on staggered gridpoints
509	kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
510	kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
511
512	tend(k,j,i) = tend(k,j,i) &
513	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
514	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
515	& ) &
516	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
517	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
518	& ) &
519	& ) * ddzw(k)
520	ENDDO
521
522	!
523	!-- Vertical diffusion at the first grid point above the surface, if the
524	!-- momentum flux at the bottom is given by the Prandtl law or if it is
525	!-- prescribed by the user.
526	!-- Difference quotient of the momentum flux is not formed over half of
527	!-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with
528	!-- other (LES) modell showed that the values of the momentum flux becomes
529	!-- too large in this case.
530	!-- The term containing w(k-1,..) (see above equation) is removed here
531	!-- because the vertical velocity is assumed to be zero at the surface.
532	IF ( use_surface_fluxes ) THEN
533	k = nzb_u_inner(j,i)+1
534	!
535	!-- Interpolate eddy diffusivities on staggered gridpoints
536	kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
537	kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
538
539	tend(k,j,i) = tend(k,j,i) &
540	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
541	& ) * ddzw(k) &
542	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
543	& + usws(j,i) &
544	& ) * ddzw(k)
545	ENDIF
546
547	!
548	!-- Vertical diffusion at the first gridpoint below the top boundary,
549	!-- if the momentum flux at the top is prescribed by the user
550	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
551	k = nzt
552	!
553	!-- Interpolate eddy diffusivities on staggered gridpoints
554	kmzp = 0.25 * &
555	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
556	kmzm = 0.25 * &
557	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
558
559	tend(k,j,i) = tend(k,j,i) &
560	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
561	& ) * ddzw(k) &
562	& + ( -uswst(j,i) &
563	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
564	& ) * ddzw(k)
565	ENDIF
566
567	END SUBROUTINE diffusion_u_ij
568
569	END MODULE diffusion_u_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |