Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 23.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |