Home

Context Navigation

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

Last change on this file since 1003 was 1002, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 15.4 KB

Line
1	MODULE diffusion_u_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: diffusion_u.f90 1002 2012-09-13 15:12: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
61
62	INTERFACE diffusion_u
63	MODULE PROCEDURE diffusion_u
64	MODULE PROCEDURE diffusion_u_ij
65	END INTERFACE diffusion_u
66
67	CONTAINS
68
69
70	!------------------------------------------------------------------------------!
71	! Call for all grid points
72	!------------------------------------------------------------------------------!
73	SUBROUTINE diffusion_u
74
75	USE arrays_3d
76	USE control_parameters
77	USE grid_variables
78	USE indices
79
80	IMPLICIT NONE
81
82	INTEGER :: i, j, k
83	REAL :: kmym, kmyp, kmzm, kmzp
84
85	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
86
87	!
88	!-- First calculate horizontal momentum flux u'v' at vertical walls,
89	!-- if neccessary
90	IF ( topography /= 'flat' ) THEN
91	CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
92	nzb_u_outer, wall_u )
93	ENDIF
94
95	DO i = nxlu, nxr
96	DO j = nys, nyn
97	!
98	!-- Compute horizontal diffusion
99	DO k = nzb_u_outer(j,i)+1, nzt
100	!
101	!-- Interpolate eddy diffusivities on staggered gridpoints
102	kmyp = 0.25 * &
103	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
104	kmym = 0.25 * &
105	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
106
107	tend(k,j,i) = tend(k,j,i) &
108	& + 2.0 * ( &
109	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
110	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
111	& ) * ddx2 &
112	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
113	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
114	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
115	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
116	& ) * ddy
117	ENDDO
118
119	!
120	!-- Wall functions at the north and south walls, respectively
121	IF ( wall_u(j,i) /= 0.0 ) THEN
122
123	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
124	kmyp = 0.25 * &
125	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
126	kmym = 0.25 * &
127	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
128
129	tend(k,j,i) = tend(k,j,i) &
130	+ 2.0 * ( &
131	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
132	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
133	) * ddx2 &
134	+ ( fyp(j,i) * ( &
135	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
136	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
137	) &
138	- fym(j,i) * ( &
139	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
140	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
141	) &
142	+ wall_u(j,i) * usvs(k,j,i) &
143	) * ddy
144	ENDDO
145	ENDIF
146
147	!
148	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
149	!-- index k starts at nzb_u_inner+2.
150	DO k = nzb_diff_u(j,i), nzt_diff
151	!
152	!-- Interpolate eddy diffusivities on staggered gridpoints
153	kmzp = 0.25 * &
154	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
155	kmzm = 0.25 * &
156	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
157
158	tend(k,j,i) = tend(k,j,i) &
159	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
160	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
161	& ) &
162	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
163	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
164	& ) &
165	& ) * ddzw(k)
166	ENDDO
167
168	!
169	!-- Vertical diffusion at the first grid point above the surface,
170	!-- if the momentum flux at the bottom is given by the Prandtl law or
171	!-- if it is prescribed by the user.
172	!-- Difference quotient of the momentum flux is not formed over half
173	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
174	!-- with other (LES) modell showed that the values of the momentum
175	!-- flux becomes too large in this case.
176	!-- The term containing w(k-1,..) (see above equation) is removed here
177	!-- because the vertical velocity is assumed to be zero at the surface.
178	IF ( use_surface_fluxes ) THEN
179	k = nzb_u_inner(j,i)+1
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 * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
189	& ) * ddzw(k) &
190	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
191	& + usws(j,i) &
192	& ) * ddzw(k)
193	ENDIF
194
195	!
196	!-- Vertical diffusion at the first gridpoint below the top boundary,
197	!-- if the momentum flux at the top is prescribed by the user
198	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
199	k = nzt
200	!
201	!-- Interpolate eddy diffusivities on staggered gridpoints
202	kmzp = 0.25 * &
203	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
204	kmzm = 0.25 * &
205	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
206
207	tend(k,j,i) = tend(k,j,i) &
208	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
209	& ) * ddzw(k) &
210	& + ( -uswst(j,i) &
211	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
212	& ) * ddzw(k)
213	ENDIF
214
215	ENDDO
216	ENDDO
217
218	END SUBROUTINE diffusion_u
219
220
221	!------------------------------------------------------------------------------!
222	! Call for grid point i,j
223	!------------------------------------------------------------------------------!
224	SUBROUTINE diffusion_u_ij( i, j )
225
226	USE arrays_3d
227	USE control_parameters
228	USE grid_variables
229	USE indices
230
231	IMPLICIT NONE
232
233	INTEGER :: i, j, k
234	REAL :: kmym, kmyp, kmzm, kmzp
235
236	REAL, DIMENSION(nzb:nzt+1) :: usvs
237
238	!
239	!-- Compute horizontal diffusion
240	DO k = nzb_u_outer(j,i)+1, nzt
241	!
242	!-- Interpolate eddy diffusivities on staggered gridpoints
243	kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
244	kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
245
246	tend(k,j,i) = tend(k,j,i) &
247	& + 2.0 * ( &
248	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
249	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
250	& ) * ddx2 &
251	& + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
252	& + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
253	& - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
254	& - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
255	& ) * ddy
256	ENDDO
257
258	!
259	!-- Wall functions at the north and south walls, respectively
260	IF ( wall_u(j,i) .NE. 0.0 ) THEN
261
262	!
263	!-- Calculate the horizontal momentum flux u'v'
264	CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i), &
265	usvs, 1.0, 0.0, 0.0, 0.0 )
266
267	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
268	kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
269	kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
270
271	tend(k,j,i) = tend(k,j,i) &
272	+ 2.0 * ( &
273	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
274	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
275	) * ddx2 &
276	+ ( fyp(j,i) * ( &
277	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
278	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
279	) &
280	- fym(j,i) * ( &
281	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
282	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
283	) &
284	+ wall_u(j,i) * usvs(k) &
285	) * ddy
286	ENDDO
287	ENDIF
288
289	!
290	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
291	!-- index k starts at nzb_u_inner+2.
292	DO k = nzb_diff_u(j,i), nzt_diff
293	!
294	!-- Interpolate eddy diffusivities on staggered gridpoints
295	kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
296	kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
297
298	tend(k,j,i) = tend(k,j,i) &
299	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
300	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
301	& ) &
302	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
303	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
304	& ) &
305	& ) * ddzw(k)
306	ENDDO
307
308	!
309	!-- Vertical diffusion at the first grid point above the surface, if the
310	!-- momentum flux at the bottom is given by the Prandtl law or if it is
311	!-- prescribed by the user.
312	!-- Difference quotient of the momentum flux is not formed over half of
313	!-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with
314	!-- other (LES) modell showed that the values of the momentum flux becomes
315	!-- too large in this case.
316	!-- The term containing w(k-1,..) (see above equation) is removed here
317	!-- because the vertical velocity is assumed to be zero at the surface.
318	IF ( use_surface_fluxes ) THEN
319	k = nzb_u_inner(j,i)+1
320	!
321	!-- Interpolate eddy diffusivities on staggered gridpoints
322	kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
323	kmzm = 0.25 * ( 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 * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
327	& ) * ddzw(k) &
328	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
329	& + usws(j,i) &
330	& ) * ddzw(k)
331	ENDIF
332
333	!
334	!-- Vertical diffusion at the first gridpoint below the top boundary,
335	!-- if the momentum flux at the top is prescribed by the user
336	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
337	k = nzt
338	!
339	!-- Interpolate eddy diffusivities on staggered gridpoints
340	kmzp = 0.25 * &
341	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
342	kmzm = 0.25 * &
343	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
344
345	tend(k,j,i) = tend(k,j,i) &
346	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
347	& ) * ddzw(k) &
348	& + ( -uswst(j,i) &
349	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
350	& ) * ddzw(k)
351	ENDIF
352
353	END SUBROUTINE diffusion_u_ij
354
355	END MODULE diffusion_u_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |