Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_w.f90 @ 1035

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

Line
1	MODULE diffusion_w_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: diffusion_w.f90 1017 2012-09-27 11:28:50Z raasch $
11	!
12	! 1015 2012-09-27 09:23:24Z raasch
13	! accelerator version (*_acc) added
14	!
15	! 1001 2012-09-13 14:08:46Z raasch
16	! arrays comunicated by module instead of parameter list
17	!
18	! 978 2012-08-09 08:28:32Z fricke
19	! outflow damping layer removed
20	! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
21	! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
22	!
23	! 667 2010-12-23 12:06:00Z suehring/gryschka
24	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
25	!
26	! 366 2009-08-25 08:06:27Z raasch
27	! bc_lr/bc_ns replaced by bc_lr_cyc/bc_ns_cyc
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
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.12 2006/02/23 10:38:03 raasch
39	! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
40	! +z0 in argument list
41	! WARNING: loops containing the MAX function are still not properly vectorized!
42	!
43	! Revision 1.1 1997/09/12 06:24:11 raasch
44	! Initial revision
45	!
46	!
47	! Description:
48	! ------------
49	! Diffusion term of the w-component
50	!------------------------------------------------------------------------------!
51
52	USE wall_fluxes_mod
53
54	PRIVATE
55	PUBLIC diffusion_w, diffusion_w_acc
56
57	INTERFACE diffusion_w
58	MODULE PROCEDURE diffusion_w
59	MODULE PROCEDURE diffusion_w_ij
60	END INTERFACE diffusion_w
61
62	INTERFACE diffusion_w_acc
63	MODULE PROCEDURE diffusion_w_acc
64	END INTERFACE diffusion_w_acc
65
66	CONTAINS
67
68
69	!------------------------------------------------------------------------------!
70	! Call for all grid points
71	!------------------------------------------------------------------------------!
72	SUBROUTINE diffusion_w
73
74	USE arrays_3d
75	USE control_parameters
76	USE grid_variables
77	USE indices
78
79	IMPLICIT NONE
80
81	INTEGER :: i, j, k
82	REAL :: kmxm, kmxp, kmym, kmyp
83
84	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
85
86
87	!
88	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
89	!-- walls, if neccessary
90	IF ( topography /= 'flat' ) THEN
91	CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
92	nzb_w_outer, wall_w_x )
93	CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
94	nzb_w_outer, wall_w_y )
95	ENDIF
96
97	DO i = nxl, nxr
98	DO j = nys, nyn
99	DO k = nzb_w_outer(j,i)+1, nzt-1
100	!
101	!-- Interpolate eddy diffusivities on staggered gridpoints
102	kmxp = 0.25 * &
103	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
104	kmxm = 0.25 * &
105	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
106	kmyp = 0.25 * &
107	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
108	kmym = 0.25 * &
109	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
110
111	tend(k,j,i) = tend(k,j,i) &
112	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
113	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
114	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
115	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
116	& ) * ddx &
117	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
118	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
119	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
120	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
121	& ) * ddy &
122	& + 2.0 * ( &
123	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
124	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
125	& ) * ddzu(k+1)
126	ENDDO
127
128	!
129	!-- Wall functions at all vertical walls, where necessary
130	IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN
131
132	DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
133	!
134	!-- Interpolate eddy diffusivities on staggered gridpoints
135	kmxp = 0.25 * &
136	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
137	kmxm = 0.25 * &
138	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
139	kmyp = 0.25 * &
140	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
141	kmym = 0.25 * &
142	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
143
144	tend(k,j,i) = tend(k,j,i) &
145	+ ( fwxp(j,i) * ( &
146	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
147	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
148	) &
149	- fwxm(j,i) * ( &
150	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
151	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
152	) &
153	+ wall_w_x(j,i) * wsus(k,j,i) &
154	) * ddx &
155	+ ( fwyp(j,i) * ( &
156	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
157	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
158	) &
159	- fwym(j,i) * ( &
160	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
161	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
162	) &
163	+ wall_w_y(j,i) * wsvs(k,j,i) &
164	) * ddy &
165	+ 2.0 * ( &
166	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
167	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
168	) * ddzu(k+1)
169	ENDDO
170	ENDIF
171
172	ENDDO
173	ENDDO
174
175	END SUBROUTINE diffusion_w
176
177
178	!------------------------------------------------------------------------------!
179	! Call for all grid points - accelerator version
180	!------------------------------------------------------------------------------!
181	SUBROUTINE diffusion_w_acc
182
183	USE arrays_3d
184	USE control_parameters
185	USE grid_variables
186	USE indices
187
188	IMPLICIT NONE
189
190	INTEGER :: i, j, k
191	REAL :: kmxm, kmxp, kmym, kmyp
192
193	!$acc declare create ( wsus, wsvs )
194	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs
195
196
197	!
198	!-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical
199	!-- walls, if neccessary
200	IF ( topography /= 'flat' ) THEN
201	CALL wall_fluxes_acc( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
202	nzb_w_outer, wall_w_x )
203	CALL wall_fluxes_acc( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
204	nzb_w_outer, wall_w_y )
205	ENDIF
206
207	!$acc kernels present ( u, v, w, km, tend, vsws, vswst ) &
208	!$acc present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y ) &
209	!$acc present ( nzb_w_inner, nzb_w_outer )
210	!$acc loop
211	DO i = nxl, nxr
212	DO j = nys, nyn
213	!$acc loop vector( 32 )
214	DO k = 1, nzt
215	IF ( k > nzb_w_outer(j,i) ) THEN
216	!
217	!-- Interpolate eddy diffusivities on staggered gridpoints
218	kmxp = 0.25 * &
219	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
220	kmxm = 0.25 * &
221	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
222	kmyp = 0.25 * &
223	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
224	kmym = 0.25 * &
225	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
226
227	tend(k,j,i) = tend(k,j,i) &
228	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
229	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
230	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
231	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
232	& ) * ddx &
233	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
234	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
235	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
236	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
237	& ) * ddy &
238	& + 2.0 * ( &
239	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
240	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
241	& ) * ddzu(k+1)
242	ENDIF
243	ENDDO
244
245	!
246	!-- Wall functions at all vertical walls, where necessary
247	!$acc loop vector( 32 )
248	DO k = 1,nzt
249
250	IF ( k > nzb_w_inner(j,i) .AND. k <= nzb_w_outer(j,i) .AND. &
251	wall_w_x(j,i) /= 0.0 .AND. wall_w_y(j,i) /= 0.0 ) THEN
252	!
253	!-- Interpolate eddy diffusivities on staggered gridpoints
254	kmxp = 0.25 * &
255	( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
256	kmxm = 0.25 * &
257	( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
258	kmyp = 0.25 * &
259	( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
260	kmym = 0.25 * &
261	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
262
263	tend(k,j,i) = tend(k,j,i) &
264	+ ( fwxp(j,i) * ( &
265	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
266	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
267	) &
268	- fwxm(j,i) * ( &
269	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
270	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
271	) &
272	+ wall_w_x(j,i) * wsus(k,j,i) &
273	) * ddx &
274	+ ( fwyp(j,i) * ( &
275	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
276	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
277	) &
278	- fwym(j,i) * ( &
279	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
280	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
281	) &
282	+ wall_w_y(j,i) * wsvs(k,j,i) &
283	) * ddy &
284	+ 2.0 * ( &
285	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
286	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
287	) * ddzu(k+1)
288	ENDIF
289	ENDDO
290
291	ENDDO
292	ENDDO
293	!$acc end kernels
294
295	END SUBROUTINE diffusion_w_acc
296
297
298	!------------------------------------------------------------------------------!
299	! Call for grid point i,j
300	!------------------------------------------------------------------------------!
301	SUBROUTINE diffusion_w_ij( i, j )
302
303	USE arrays_3d
304	USE control_parameters
305	USE grid_variables
306	USE indices
307
308	IMPLICIT NONE
309
310	INTEGER :: i, j, k
311	REAL :: kmxm, kmxp, kmym, kmyp
312
313	REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs
314
315
316	DO k = nzb_w_outer(j,i)+1, nzt-1
317	!
318	!-- Interpolate eddy diffusivities on staggered gridpoints
319	kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
320	kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
321	kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
322	kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
323
324	tend(k,j,i) = tend(k,j,i) &
325	& + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
326	& + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
327	& - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
328	& - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
329	& ) * ddx &
330	& + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
331	& + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
332	& - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
333	& - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
334	& ) * ddy &
335	& + 2.0 * ( &
336	& km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
337	& - km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
338	& ) * ddzu(k+1)
339	ENDDO
340
341	!
342	!-- Wall functions at all vertical walls, where necessary
343	IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN
344
345	!
346	!-- Calculate the horizontal momentum fluxes w'u' and/or w'v'
347	IF ( wall_w_x(j,i) /= 0.0 ) THEN
348	CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
349	wsus, 0.0, 0.0, 0.0, 1.0 )
350	ELSE
351	wsus = 0.0
352	ENDIF
353
354	IF ( wall_w_y(j,i) /= 0.0 ) THEN
355	CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
356	wsvs, 0.0, 0.0, 1.0, 0.0 )
357	ELSE
358	wsvs = 0.0
359	ENDIF
360
361	DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
362	!
363	!-- Interpolate eddy diffusivities on staggered gridpoints
364	kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
365	kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
366	kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
367	kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
368
369	tend(k,j,i) = tend(k,j,i) &
370	+ ( fwxp(j,i) * ( &
371	kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx &
372	+ kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
373	) &
374	- fwxm(j,i) * ( &
375	kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
376	+ kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
377	) &
378	+ wall_w_x(j,i) * wsus(k) &
379	) * ddx &
380	+ ( fwyp(j,i) * ( &
381	kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy &
382	+ kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
383	) &
384	- fwym(j,i) * ( &
385	kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy &
386	+ kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
387	) &
388	+ wall_w_y(j,i) * wsvs(k) &
389	) * ddy &
390	+ 2.0 * ( &
391	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
392	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
393	) * ddzu(k+1)
394	ENDDO
395	ENDIF
396
397	END SUBROUTINE diffusion_w_ij
398
399	END MODULE diffusion_w_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |