Home

Context Navigation

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

Last change on this file since 993 was 979, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 15.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |