Home

Context Navigation

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

Last change on this file since 978 was 978, checked in by fricke, 12 years ago
merge fricke branch back into trunk
Property svn:keywords set to `Id`
File size: 15.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |