Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |