Home

Context Navigation

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

Last change on this file since 4148 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 21.1 KB

Line
1	!> @file diffusion_w.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! 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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_w.f90 3655 2019-01-07 16:51:22Z suehring $
27	! OpenACC port for SPEC
28	!
29	! 3547 2018-11-21 13:21:24Z suehring
30	! variables documented
31	!
32	! 3241 2018-09-12 15:02:00Z raasch
33	! unused variables removed
34	!
35	! 2718 2018-01-02 08:49:38Z maronga
36	! Corrected "Former revisions" section
37	!
38	! 2696 2017-12-14 17:12:51Z kanani
39	! Change in file header (GPL part)
40	!
41	! 2233 2017-05-30 18:08:54Z suehring
42	!
43	! 2232 2017-05-30 17:47:52Z suehring
44	! Adjustments to new topography and surface concept
45	!
46	! 2118 2017-01-17 16:38:49Z raasch
47	! OpenACC version of subroutine removed
48	!
49	! 2037 2016-10-26 11:15:40Z knoop
50	! Anelastic approximation implemented
51	!
52	! 2000 2016-08-20 18:09:15Z knoop
53	! Forced header and separation lines into 80 columns
54	!
55	! 1873 2016-04-18 14:50:06Z maronga
56	! Module renamed (removed _mod)
57	!
58	! 1850 2016-04-08 13:29:27Z maronga
59	! Module renamed
60	!
61	! 1682 2015-10-07 23:56:08Z knoop
62	! Code annotations made doxygen readable
63	!
64	! 1374 2014-04-25 12:55:07Z raasch
65	! vsws + vswst removed from acc-present-list
66	!
67	! 1353 2014-04-08 15:21:23Z heinze
68	! REAL constants provided with KIND-attribute
69	!
70	! 1340 2014-03-25 19:45:13Z kanani
71	! REAL constants defined as wp-kind
72	!
73	! 1322 2014-03-20 16:38:49Z raasch
74	! REAL constants defined as wp-kind
75	!
76	! 1320 2014-03-20 08:40:49Z raasch
77	! ONLY-attribute added to USE-statements,
78	! kind-parameters added to all INTEGER and REAL declaration statements,
79	! kinds are defined in new module kinds,
80	! revision history before 2012 removed,
81	! comment fields (!:) to be used for variable explanations added to
82	! all variable declaration statements
83	!
84	! 1257 2013-11-08 15:18:40Z raasch
85	! openacc loop and loop vector clauses removed, declare create moved after
86	! the FORTRAN declaration statement
87	!
88	! 1128 2013-04-12 06:19:32Z raasch
89	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
90	! j_north
91	!
92	! 1036 2012-10-22 13:43:42Z raasch
93	! code put under GPL (PALM 3.9)
94	!
95	! 1015 2012-09-27 09:23:24Z raasch
96	! accelerator version (*_acc) added
97	!
98	! 1001 2012-09-13 14:08:46Z raasch
99	! arrays comunicated by module instead of parameter list
100	!
101	! 978 2012-08-09 08:28:32Z fricke
102	! outflow damping layer removed
103	! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
104	! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
105	!
106	! Revision 1.1 1997/09/12 06:24:11 raasch
107	! Initial revision
108	!
109	!
110	! Description:
111	! ------------
112	!> Diffusion term of the w-component
113	!------------------------------------------------------------------------------!
114	MODULE diffusion_w_mod
115
116
117	PRIVATE
118	PUBLIC diffusion_w
119
120	INTERFACE diffusion_w
121	MODULE PROCEDURE diffusion_w
122	MODULE PROCEDURE diffusion_w_ij
123	END INTERFACE diffusion_w
124
125	CONTAINS
126
127
128	!------------------------------------------------------------------------------!
129	! Description:
130	! ------------
131	!> Call for all grid points
132	!------------------------------------------------------------------------------!
133	SUBROUTINE diffusion_w
134
135	USE arrays_3d, &
136	ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
137
138	USE grid_variables, &
139	ONLY : ddx, ddy
140
141	USE indices, &
142	ONLY : nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
143
144	USE kinds
145
146	USE surface_mod, &
147	ONLY : surf_def_v, surf_lsm_v, surf_usm_v
148
149	IMPLICIT NONE
150
151	INTEGER(iwp) :: i !< running index x direction
152	INTEGER(iwp) :: j !< running index y direction
153	INTEGER(iwp) :: k !< running index z direction
154	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
155	INTEGER(iwp) :: m !< running index surface elements
156	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
157	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
158
159	REAL(wp) :: flag !< flag to mask topography grid points
160	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
161	REAL(wp) :: kmxp !<diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
162	REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
163	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
164	REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point
165	REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point
166	REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point
167	REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point
168
169
170
171	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
172	!$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
173	!$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
174	!$ACC PRESENT(wall_flags_0, km) &
175	!$ACC PRESENT(u, v, w) &
176	!$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
177	!$ACC PRESENT(surf_def_v(0:3)) &
178	!$ACC PRESENT(surf_lsm_v(0:3)) &
179	!$ACC PRESENT(surf_usm_v(0:3)) &
180	!$ACC PRESENT(tend)
181	DO i = nxl, nxr
182	DO j = nys, nyn
183	DO k = nzb+1, nzt-1
184	!
185	!-- Predetermine flag to mask topography and wall-bounded grid points.
186	flag = MERGE( 1.0_wp, 0.0_wp, &
187	BTEST( wall_flags_0(k,j,i), 3 ) )
188	mask_east = MERGE( 1.0_wp, 0.0_wp, &
189	BTEST( wall_flags_0(k,j,i+1), 3 ) )
190	mask_west = MERGE( 1.0_wp, 0.0_wp, &
191	BTEST( wall_flags_0(k,j,i-1), 3 ) )
192	mask_south = MERGE( 1.0_wp, 0.0_wp, &
193	BTEST( wall_flags_0(k,j-1,i), 3 ) )
194	mask_north = MERGE( 1.0_wp, 0.0_wp, &
195	BTEST( wall_flags_0(k,j+1,i), 3 ) )
196	!
197	!-- Interpolate eddy diffusivities on staggered gridpoints
198	kmxp = 0.25_wp * ( km(k,j,i) + km(k,j,i+1) + &
199	km(k+1,j,i) + km(k+1,j,i+1) )
200	kmxm = 0.25_wp * ( km(k,j,i) + km(k,j,i-1) + &
201	km(k+1,j,i) + km(k+1,j,i-1) )
202	kmyp = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &
203	km(k,j+1,i) + km(k+1,j+1,i) )
204	kmym = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &
205	km(k,j-1,i) + km(k+1,j-1,i) )
206
207	tend(k,j,i) = tend(k,j,i) &
208	+ ( mask_east * kmxp * ( &
209	( w(k,j,i+1) - w(k,j,i) ) * ddx &
210	+ ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
211	) &
212	- mask_west * kmxm * ( &
213	( w(k,j,i) - w(k,j,i-1) ) * ddx &
214	+ ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
215	) &
216	) * ddx * flag &
217	+ ( mask_north * kmyp * ( &
218	( w(k,j+1,i) - w(k,j,i) ) * ddy &
219	+ ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
220	) &
221	- mask_south * kmym * ( &
222	( w(k,j,i) - w(k,j-1,i) ) * ddy &
223	+ ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
224	) &
225	) * ddy * flag &
226	+ 2.0_wp * ( &
227	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
228	* rho_air(k+1) &
229	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
230	* rho_air(k) &
231	) * ddzu(k+1) * drho_air_zw(k) * flag
232	ENDDO
233
234	!
235	!-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
236	!-- surfaces. Note, in the the flat case, loops won't be entered as
237	!-- start_index > end_index. Furtermore, note, no vertical natural surfaces
238	!-- so far.
239	!-- Default-type surfaces
240	DO l = 0, 1
241	surf_s = surf_def_v(l)%start_index(j,i)
242	surf_e = surf_def_v(l)%end_index(j,i)
243	DO m = surf_s, surf_e
244	k = surf_def_v(l)%k(m)
245	tend(k,j,i) = tend(k,j,i) + &
246	surf_def_v(l)%mom_flux_w(m) * ddy
247	ENDDO
248	ENDDO
249	!
250	!-- Natural-type surfaces
251	DO l = 0, 1
252	surf_s = surf_lsm_v(l)%start_index(j,i)
253	surf_e = surf_lsm_v(l)%end_index(j,i)
254	DO m = surf_s, surf_e
255	k = surf_lsm_v(l)%k(m)
256	tend(k,j,i) = tend(k,j,i) + &
257	surf_lsm_v(l)%mom_flux_w(m) * ddy
258	ENDDO
259	ENDDO
260	!
261	!-- Urban-type surfaces
262	DO l = 0, 1
263	surf_s = surf_usm_v(l)%start_index(j,i)
264	surf_e = surf_usm_v(l)%end_index(j,i)
265	DO m = surf_s, surf_e
266	k = surf_usm_v(l)%k(m)
267	tend(k,j,i) = tend(k,j,i) + &
268	surf_usm_v(l)%mom_flux_w(m) * ddy
269	ENDDO
270	ENDDO
271	!
272	!-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
273	!-- surface.
274	!-- Default-type surfaces
275	DO l = 2, 3
276	surf_s = surf_def_v(l)%start_index(j,i)
277	surf_e = surf_def_v(l)%end_index(j,i)
278	DO m = surf_s, surf_e
279	k = surf_def_v(l)%k(m)
280	tend(k,j,i) = tend(k,j,i) + &
281	surf_def_v(l)%mom_flux_w(m) * ddx
282	ENDDO
283	ENDDO
284	!
285	!-- Natural-type surfaces
286	DO l = 2, 3
287	surf_s = surf_lsm_v(l)%start_index(j,i)
288	surf_e = surf_lsm_v(l)%end_index(j,i)
289	DO m = surf_s, surf_e
290	k = surf_lsm_v(l)%k(m)
291	tend(k,j,i) = tend(k,j,i) + &
292	surf_lsm_v(l)%mom_flux_w(m) * ddx
293	ENDDO
294	ENDDO
295	!
296	!-- Urban-type surfaces
297	DO l = 2, 3
298	surf_s = surf_usm_v(l)%start_index(j,i)
299	surf_e = surf_usm_v(l)%end_index(j,i)
300	DO m = surf_s, surf_e
301	k = surf_usm_v(l)%k(m)
302	tend(k,j,i) = tend(k,j,i) + &
303	surf_usm_v(l)%mom_flux_w(m) * ddx
304	ENDDO
305	ENDDO
306
307	ENDDO
308	ENDDO
309
310	END SUBROUTINE diffusion_w
311
312
313	!------------------------------------------------------------------------------!
314	! Description:
315	! ------------
316	!> Call for grid point i,j
317	!------------------------------------------------------------------------------!
318	SUBROUTINE diffusion_w_ij( i, j )
319
320	USE arrays_3d, &
321	ONLY : ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
322
323	USE grid_variables, &
324	ONLY : ddx, ddy
325
326	USE indices, &
327	ONLY : nzb, nzt, wall_flags_0
328
329	USE kinds
330
331	USE surface_mod, &
332	ONLY : surf_def_v, surf_lsm_v, surf_usm_v
333
334	IMPLICIT NONE
335
336
337	INTEGER(iwp) :: i !< running index x direction
338	INTEGER(iwp) :: j !< running index y direction
339	INTEGER(iwp) :: k !< running index z direction
340	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
341	INTEGER(iwp) :: m !< running index surface elements
342	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
343	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
344
345	REAL(wp) :: flag !< flag to mask topography grid points
346	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
347	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
348	REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
349	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
350	REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point
351	REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point
352	REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point
353	REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point
354
355
356	DO k = nzb+1, nzt-1
357	!
358	!-- Predetermine flag to mask topography and wall-bounded grid points.
359	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3 ) )
360	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) )
361	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) )
362	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) )
363	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) )
364	!
365	!-- Interpolate eddy diffusivities on staggered gridpoints
366	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
367	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
368	kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
369	kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
370
371	tend(k,j,i) = tend(k,j,i) &
372	+ ( mask_east * kmxp * ( &
373	( w(k,j,i+1) - w(k,j,i) ) * ddx &
374	+ ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
375	) &
376	- mask_west * kmxm * ( &
377	( w(k,j,i) - w(k,j,i-1) ) * ddx &
378	+ ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
379	) &
380	) * ddx * flag &
381	+ ( mask_north * kmyp * ( &
382	( w(k,j+1,i) - w(k,j,i) ) * ddy &
383	+ ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
384	) &
385	- mask_south * kmym * ( &
386	( w(k,j,i) - w(k,j-1,i) ) * ddy &
387	+ ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
388	) &
389	) * ddy * flag &
390	+ 2.0_wp * ( &
391	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
392	* rho_air(k+1) &
393	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
394	* rho_air(k) &
395	) * ddzu(k+1) * drho_air_zw(k) * flag
396	ENDDO
397	!
398	!-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
399	!-- surfaces. Note, in the the flat case, loops won't be entered as
400	!-- start_index > end_index. Furtermore, note, no vertical natural surfaces
401	!-- so far.
402	!-- Default-type surfaces
403	DO l = 0, 1
404	surf_s = surf_def_v(l)%start_index(j,i)
405	surf_e = surf_def_v(l)%end_index(j,i)
406	DO m = surf_s, surf_e
407	k = surf_def_v(l)%k(m)
408	tend(k,j,i) = tend(k,j,i) + &
409	surf_def_v(l)%mom_flux_w(m) * ddy
410	ENDDO
411	ENDDO
412	!
413	!-- Natural-type surfaces
414	DO l = 0, 1
415	surf_s = surf_lsm_v(l)%start_index(j,i)
416	surf_e = surf_lsm_v(l)%end_index(j,i)
417	DO m = surf_s, surf_e
418	k = surf_lsm_v(l)%k(m)
419	tend(k,j,i) = tend(k,j,i) + &
420	surf_lsm_v(l)%mom_flux_w(m) * ddy
421	ENDDO
422	ENDDO
423	!
424	!-- Urban-type surfaces
425	DO l = 0, 1
426	surf_s = surf_usm_v(l)%start_index(j,i)
427	surf_e = surf_usm_v(l)%end_index(j,i)
428	DO m = surf_s, surf_e
429	k = surf_usm_v(l)%k(m)
430	tend(k,j,i) = tend(k,j,i) + &
431	surf_usm_v(l)%mom_flux_w(m) * ddy
432	ENDDO
433	ENDDO
434	!
435	!-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
436	!-- surfaces.
437	!-- Default-type surfaces
438	DO l = 2, 3
439	surf_s = surf_def_v(l)%start_index(j,i)
440	surf_e = surf_def_v(l)%end_index(j,i)
441	DO m = surf_s, surf_e
442	k = surf_def_v(l)%k(m)
443	tend(k,j,i) = tend(k,j,i) + &
444	surf_def_v(l)%mom_flux_w(m) * ddx
445	ENDDO
446	ENDDO
447	!
448	!-- Natural-type surfaces
449	DO l = 2, 3
450	surf_s = surf_lsm_v(l)%start_index(j,i)
451	surf_e = surf_lsm_v(l)%end_index(j,i)
452	DO m = surf_s, surf_e
453	k = surf_lsm_v(l)%k(m)
454	tend(k,j,i) = tend(k,j,i) + &
455	surf_lsm_v(l)%mom_flux_w(m) * ddx
456	ENDDO
457	ENDDO
458	!
459	!-- Urban-type surfaces
460	DO l = 2, 3
461	surf_s = surf_usm_v(l)%start_index(j,i)
462	surf_e = surf_usm_v(l)%end_index(j,i)
463	DO m = surf_s, surf_e
464	k = surf_usm_v(l)%k(m)
465	tend(k,j,i) = tend(k,j,i) + &
466	surf_usm_v(l)%mom_flux_w(m) * ddx
467	ENDDO
468	ENDDO
469
470
471	END SUBROUTINE diffusion_w_ij
472
473	END MODULE diffusion_w_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |