Home

Context Navigation

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

Last change on this file since 1046 was 1037, checked in by raasch, 12 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 1037 2012-10-22 14:10:22Z maronga $
27	!
28	! 1036 2012-10-22 13:43:42Z raasch
29	! code put under GPL (PALM 3.9)
30	!
31	! 1015 2012-09-27 09:23:24Z raasch
32	! accelerator version (*_acc) added
33	!
34	! 1001 2012-09-13 14:08:46Z raasch
35	! arrays comunicated by module instead of parameter list
36	!
37	! 978 2012-08-09 08:28:32Z fricke
38	! outflow damping layer removed
39	! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
40	!
41	! 667 2010-12-23 12:06:00Z suehring/gryschka
42	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
43	!
44	! 366 2009-08-25 08:06:27Z raasch
45	! bc_ns replaced by bc_ns_cyc
46	!
47	! 106 2007-08-16 14:30:26Z raasch
48	! Momentumflux at top (uswst) included as boundary condition,
49	! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
50	!
51	! 75 2007-03-22 09:54:05Z raasch
52	! Wall functions now include diabatic conditions, call of routine wall_fluxes,
53	! z0 removed from argument list, uxrp eliminated
54	!
55	! 20 2007-02-26 00:12:32Z raasch
56	! Bugfix: ddzw dimensioned 1:nzt"+1"
57	!
58	! RCS Log replace by Id keyword, revision history cleaned up
59	!
60	! Revision 1.15 2006/02/23 10:35:35 raasch
61	! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
62	! or nzb_diff_u, respectively, in vertical diffusion,
63	! wall functions added for north and south walls, +z0 in argument list,
64	! terms containing w(k-1,..) are removed from the Prandtl-layer equation
65	! because they cause errors at the edges of topography
66	! WARNING: loops containing the MAX function are still not properly vectorized!
67	!
68	! Revision 1.1 1997/09/12 06:23:51 raasch
69	! Initial revision
70	!
71	!
72	! Description:
73	! ------------
74	! Diffusion term of the u-component
75	! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
76	! and slows down the speed on NEC about 5-10%
77	!------------------------------------------------------------------------------!
78
79	USE wall_fluxes_mod
80
81	PRIVATE
82	PUBLIC diffusion_u, diffusion_u_acc
83
84	INTERFACE diffusion_u
85	MODULE PROCEDURE diffusion_u
86	MODULE PROCEDURE diffusion_u_ij
87	END INTERFACE diffusion_u
88
89	INTERFACE diffusion_u_acc
90	MODULE PROCEDURE diffusion_u_acc
91	END INTERFACE diffusion_u_acc
92
93	CONTAINS
94
95
96	!------------------------------------------------------------------------------!
97	! Call for all grid points
98	!------------------------------------------------------------------------------!
99	SUBROUTINE diffusion_u
100
101	USE arrays_3d
102	USE control_parameters
103	USE grid_variables
104	USE indices
105
106	IMPLICIT NONE
107
108	INTEGER :: i, j, k
109	REAL :: kmym, kmyp, kmzm, kmzp
110
111	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
112
113	!
114	!-- First calculate horizontal momentum flux u'v' at vertical walls,
115	!-- if neccessary
116	IF ( topography /= 'flat' ) THEN
117	CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
118	nzb_u_outer, wall_u )
119	ENDIF
120
121	DO i = nxlu, nxr
122	DO j = nys, nyn
123	!
124	!-- Compute horizontal diffusion
125	DO k = nzb_u_outer(j,i)+1, nzt
126	!
127	!-- Interpolate eddy diffusivities on staggered gridpoints
128	kmyp = 0.25 * &
129	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
130	kmym = 0.25 * &
131	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
132
133	tend(k,j,i) = tend(k,j,i) &
134	& + 2.0 * ( &
135	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
136	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
137	& ) * ddx2 &
138	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
139	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
140	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
141	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
142	& ) * ddy
143	ENDDO
144
145	!
146	!-- Wall functions at the north and south walls, respectively
147	IF ( wall_u(j,i) /= 0.0 ) THEN
148
149	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
150	kmyp = 0.25 * &
151	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
152	kmym = 0.25 * &
153	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
154
155	tend(k,j,i) = tend(k,j,i) &
156	+ 2.0 * ( &
157	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
158	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
159	) * ddx2 &
160	+ ( fyp(j,i) * ( &
161	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
162	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
163	) &
164	- fym(j,i) * ( &
165	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
166	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
167	) &
168	+ wall_u(j,i) * usvs(k,j,i) &
169	) * ddy
170	ENDDO
171	ENDIF
172
173	!
174	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
175	!-- index k starts at nzb_u_inner+2.
176	DO k = nzb_diff_u(j,i), nzt_diff
177	!
178	!-- Interpolate eddy diffusivities on staggered gridpoints
179	kmzp = 0.25 * &
180	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
181	kmzm = 0.25 * &
182	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
183
184	tend(k,j,i) = tend(k,j,i) &
185	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
186	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
187	& ) &
188	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
189	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
190	& ) &
191	& ) * ddzw(k)
192	ENDDO
193
194	!
195	!-- Vertical diffusion at the first grid point above the surface,
196	!-- if the momentum flux at the bottom is given by the Prandtl law or
197	!-- if it is prescribed by the user.
198	!-- Difference quotient of the momentum flux is not formed over half
199	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
200	!-- with other (LES) modell showed that the values of the momentum
201	!-- flux becomes too large in this case.
202	!-- The term containing w(k-1,..) (see above equation) is removed here
203	!-- because the vertical velocity is assumed to be zero at the surface.
204	IF ( use_surface_fluxes ) THEN
205	k = nzb_u_inner(j,i)+1
206	!
207	!-- Interpolate eddy diffusivities on staggered gridpoints
208	kmzp = 0.25 * &
209	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
210	kmzm = 0.25 * &
211	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
212
213	tend(k,j,i) = tend(k,j,i) &
214	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
215	& ) * ddzw(k) &
216	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
217	& + usws(j,i) &
218	& ) * ddzw(k)
219	ENDIF
220
221	!
222	!-- Vertical diffusion at the first gridpoint below the top boundary,
223	!-- if the momentum flux at the top is prescribed by the user
224	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
225	k = nzt
226	!
227	!-- Interpolate eddy diffusivities on staggered gridpoints
228	kmzp = 0.25 * &
229	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
230	kmzm = 0.25 * &
231	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
232
233	tend(k,j,i) = tend(k,j,i) &
234	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
235	& ) * ddzw(k) &
236	& + ( -uswst(j,i) &
237	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
238	& ) * ddzw(k)
239	ENDIF
240
241	ENDDO
242	ENDDO
243
244	END SUBROUTINE diffusion_u
245
246
247	!------------------------------------------------------------------------------!
248	! Call for all grid points - accelerator version
249	!------------------------------------------------------------------------------!
250	SUBROUTINE diffusion_u_acc
251
252	USE arrays_3d
253	USE control_parameters
254	USE grid_variables
255	USE indices
256
257	IMPLICIT NONE
258
259	INTEGER :: i, j, k
260	REAL :: kmym, kmyp, kmzm, kmzp
261
262	!$acc declare create ( usvs )
263	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
264
265	!
266	!-- First calculate horizontal momentum flux u'v' at vertical walls,
267	!-- if neccessary
268	IF ( topography /= 'flat' ) THEN
269	CALL wall_fluxes_acc( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
270	nzb_u_outer, wall_u )
271	ENDIF
272
273	!$acc kernels present ( u, v, w, km, tend, usws, uswst ) &
274	!$acc present ( ddzu, ddzw, fym, fyp, wall_u ) &
275	!$acc present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
276	!$acc loop
277	DO i = nxlu, nxr
278	DO j = nys, nyn
279	!
280	!-- Compute horizontal diffusion
281	!$acc loop vector(32)
282	DO k = 1, nzt
283	IF ( k > nzb_u_outer(j,i) ) THEN
284	!
285	!-- Interpolate eddy diffusivities on staggered gridpoints
286	kmyp = 0.25 * &
287	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
288	kmym = 0.25 * &
289	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
290
291	tend(k,j,i) = tend(k,j,i) &
292	& + 2.0 * ( &
293	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
294	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
295	& ) * ddx2 &
296	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
297	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
298	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
299	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
300	& ) * ddy
301	ENDIF
302	ENDDO
303
304	!
305	!-- Wall functions at the north and south walls, respectively
306	!$acc loop vector(32)
307	DO k = 1, nzt
308	IF( k > nzb_u_inner(j,i) .AND. k <= nzb_u_outer(j,i) .AND. &
309	wall_u(j,i) /= 0.0 ) THEN
310
311	kmyp = 0.25 * &
312	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
313	kmym = 0.25 * &
314	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
315
316	tend(k,j,i) = tend(k,j,i) &
317	+ 2.0 * ( &
318	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
319	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
320	) * ddx2 &
321	+ ( fyp(j,i) * ( &
322	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
323	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
324	) &
325	- fym(j,i) * ( &
326	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
327	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
328	) &
329	+ wall_u(j,i) * usvs(k,j,i) &
330	) * ddy
331	ENDIF
332	ENDDO
333
334	!
335	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
336	!-- index k starts at nzb_u_inner+2.
337	!$acc loop vector(32)
338	DO k = 1, nzt_diff
339	IF ( k >= nzb_diff_u(j,i) ) THEN
340	!
341	!-- Interpolate eddy diffusivities on staggered gridpoints
342	kmzp = 0.25 * &
343	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
344	kmzm = 0.25 * &
345	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
346
347	tend(k,j,i) = tend(k,j,i) &
348	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&
349	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
350	& ) &
351	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)&
352	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
353	& ) &
354	& ) * ddzw(k)
355	ENDIF
356	ENDDO
357
358	ENDDO
359	ENDDO
360
361	!
362	!-- Vertical diffusion at the first grid point above the surface,
363	!-- if the momentum flux at the bottom is given by the Prandtl law or
364	!-- if it is prescribed by the user.
365	!-- Difference quotient of the momentum flux is not formed over half
366	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
367	!-- with other (LES) modell showed that the values of the momentum
368	!-- flux becomes too large in this case.
369	!-- The term containing w(k-1,..) (see above equation) is removed here
370	!-- because the vertical velocity is assumed to be zero at the surface.
371	IF ( use_surface_fluxes ) THEN
372
373	!$acc loop
374	DO i = nxlu, nxr
375	!$acc loop vector(32)
376	DO j = nys, nyn
377
378	k = nzb_u_inner(j,i)+1
379	!
380	!-- Interpolate eddy diffusivities on staggered gridpoints
381	kmzp = 0.25 * &
382	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
383	kmzm = 0.25 * &
384	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
385
386	tend(k,j,i) = tend(k,j,i) &
387	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
388	& ) * ddzw(k) &
389	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
390	& + usws(j,i) &
391	& ) * ddzw(k)
392	ENDDO
393	ENDDO
394
395	ENDIF
396
397	!
398	!-- Vertical diffusion at the first gridpoint below the top boundary,
399	!-- if the momentum flux at the top is prescribed by the user
400	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
401
402	k = nzt
403
404	!$acc loop
405	DO i = nxlu, nxr
406	!$acc loop vector(32)
407	DO j = nys, nyn
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 |