Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |