Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |