Home

Context Navigation

diffusion_w.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 10 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:

openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:

bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: