Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 18.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |