Home

Context Navigation

diffusion_w.f90 @ 4867

Last change on this file since 4867 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 19.3 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 terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2021 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: diffusion_w.f90 4828 2021-01-05 11:21:41Z moh.hefny $
26	! file re-formatted to follow the PALM coding standard
27	!
28	! 4360 2020-01-07 11:25:50Z suehring
29	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
30	! information used in wall_flags_static_0
31	!
32	! 4329 2019-12-10 15:46:36Z motisi
33	! Renamed wall_flags_0 to wall_flags_static_0
34	!
35	! 4182 2019-08-22 15:20:23Z scharf
36	! Corrected "Former revisions" section
37	!
38	! 3655 2019-01-07 16:51:22Z knoop
39	! OpenACC port for SPEC
40	!
41	! Revision 1.1 1997/09/12 06:24:11 raasch
42	! Initial revision
43	!
44	!
45	! Description:
46	! ------------
47	!> Diffusion term of the w-component
48	!--------------------------------------------------------------------------------------------------!
49	MODULE diffusion_w_mod
50
51
52	PRIVATE
53	PUBLIC diffusion_w
54
55	INTERFACE diffusion_w
56	MODULE PROCEDURE diffusion_w
57	MODULE PROCEDURE diffusion_w_ij
58	END INTERFACE diffusion_w
59
60	CONTAINS
61
62
63	!--------------------------------------------------------------------------------------------------!
64	! Description:
65	! ------------
66	!> Call for all grid points
67	!--------------------------------------------------------------------------------------------------!
68	SUBROUTINE diffusion_w
69
70	USE arrays_3d, &
71	ONLY : ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
72
73	USE grid_variables, &
74	ONLY : ddx, ddy
75
76	USE indices, &
77	ONLY : nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
78
79	USE kinds
80
81	USE surface_mod, &
82	ONLY : surf_def_v, surf_lsm_v, surf_usm_v
83
84	IMPLICIT NONE
85
86	INTEGER(iwp) :: i !< running index x direction
87	INTEGER(iwp) :: j !< running index y direction
88	INTEGER(iwp) :: k !< running index z direction
89	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
90	INTEGER(iwp) :: m !< running index surface elements
91	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
92	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
93
94	REAL(wp) :: flag !< flag to mask topography grid points
95	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
96	REAL(wp) :: kmxp !<diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
97	REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
98	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
99	REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point
100	REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point
101	REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point
102	REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point
103
104
105
106	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
107	!$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
108	!$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
109	!$ACC PRESENT(wall_flags_total_0, km) &
110	!$ACC PRESENT(u, v, w) &
111	!$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
112	!$ACC PRESENT(surf_def_v(0:3)) &
113	!$ACC PRESENT(surf_lsm_v(0:3)) &
114	!$ACC PRESENT(surf_usm_v(0:3)) &
115	!$ACC PRESENT(tend)
116	DO i = nxl, nxr
117	DO j = nys, nyn
118	DO k = nzb+1, nzt-1
119	!
120	!-- Predetermine flag to mask topography and wall-bounded grid points.
121	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
122	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
123	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
124	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
125	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
126	!
127	!-- Interpolate eddy diffusivities on staggered gridpoints
128	kmxp = 0.25_wp * ( km(k,j,i) + km(k,j,i+1) + &
129	km(k+1,j,i) + km(k+1,j,i+1) )
130	kmxm = 0.25_wp * ( km(k,j,i) + km(k,j,i-1) + &
131	km(k+1,j,i) + km(k+1,j,i-1) )
132	kmyp = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &
133	km(k,j+1,i) + km(k+1,j+1,i) )
134	kmym = 0.25_wp * ( km(k,j,i) + km(k+1,j,i) + &
135	km(k,j-1,i) + km(k+1,j-1,i) )
136
137	tend(k,j,i) = tend(k,j,i) &
138	+ ( mask_east * kmxp * ( &
139	( w(k,j,i+1) - w(k,j,i) ) * ddx &
140	+ ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
141	) &
142	- mask_west * kmxm * ( &
143	( w(k,j,i) - w(k,j,i-1) ) * ddx &
144	+ ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
145	) &
146	) * ddx * flag &
147	+ ( mask_north * kmyp * ( &
148	( w(k,j+1,i) - w(k,j,i) ) * ddy &
149	+ ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
150	) &
151	- mask_south * kmym * ( &
152	( w(k,j,i) - w(k,j-1,i) ) * ddy &
153	+ ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
154	) &
155	) * ddy * flag &
156	+ 2.0_wp * ( &
157	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
158	* rho_air(k+1) &
159	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
160	* rho_air(k) &
161	) * ddzu(k+1) * drho_air_zw(k) * flag
162	ENDDO
163
164	!
165	!-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces.
166	!-- Note, in the the flat case, loops won't be entered as start_index > end_index.
167	!-- Furtermore, note, no vertical natural surfaces so far.
168	!-- Default-type surfaces
169	DO l = 0, 1
170	surf_s = surf_def_v(l)%start_index(j,i)
171	surf_e = surf_def_v(l)%end_index(j,i)
172	DO m = surf_s, surf_e
173	k = surf_def_v(l)%k(m)
174	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
175	ENDDO
176	ENDDO
177	!
178	!-- Natural-type surfaces
179	DO l = 0, 1
180	surf_s = surf_lsm_v(l)%start_index(j,i)
181	surf_e = surf_lsm_v(l)%end_index(j,i)
182	DO m = surf_s, surf_e
183	k = surf_lsm_v(l)%k(m)
184	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
185	ENDDO
186	ENDDO
187	!
188	!-- Urban-type surfaces
189	DO l = 0, 1
190	surf_s = surf_usm_v(l)%start_index(j,i)
191	surf_e = surf_usm_v(l)%end_index(j,i)
192	DO m = surf_s, surf_e
193	k = surf_usm_v(l)%k(m)
194	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
195	ENDDO
196	ENDDO
197	!
198	!-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surface.
199	!-- Default-type surfaces
200	DO l = 2, 3
201	surf_s = surf_def_v(l)%start_index(j,i)
202	surf_e = surf_def_v(l)%end_index(j,i)
203	DO m = surf_s, surf_e
204	k = surf_def_v(l)%k(m)
205	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
206	ENDDO
207	ENDDO
208	!
209	!-- Natural-type surfaces
210	DO l = 2, 3
211	surf_s = surf_lsm_v(l)%start_index(j,i)
212	surf_e = surf_lsm_v(l)%end_index(j,i)
213	DO m = surf_s, surf_e
214	k = surf_lsm_v(l)%k(m)
215	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
216	ENDDO
217	ENDDO
218	!
219	!-- Urban-type surfaces
220	DO l = 2, 3
221	surf_s = surf_usm_v(l)%start_index(j,i)
222	surf_e = surf_usm_v(l)%end_index(j,i)
223	DO m = surf_s, surf_e
224	k = surf_usm_v(l)%k(m)
225	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
226	ENDDO
227	ENDDO
228
229	ENDDO
230	ENDDO
231
232	END SUBROUTINE diffusion_w
233
234
235	!--------------------------------------------------------------------------------------------------!
236	! Description:
237	! ------------
238	!> Call for grid point i,j
239	!--------------------------------------------------------------------------------------------------!
240	SUBROUTINE diffusion_w_ij( i, j )
241
242	USE arrays_3d, &
243	ONLY : ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
244
245	USE grid_variables, &
246	ONLY : ddx, ddy
247
248	USE indices, &
249	ONLY : nzb, nzt, wall_flags_total_0
250
251	USE kinds
252
253	USE surface_mod, &
254	ONLY : surf_def_v, surf_lsm_v, surf_usm_v
255
256	IMPLICIT NONE
257
258
259	INTEGER(iwp) :: i !< running index x direction
260	INTEGER(iwp) :: j !< running index y direction
261	INTEGER(iwp) :: k !< running index z direction
262	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
263	INTEGER(iwp) :: m !< running index surface elements
264	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
265	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
266
267	REAL(wp) :: flag !< flag to mask topography grid points
268	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
269	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
270	REAL(wp) :: kmym !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
271	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
272	REAL(wp) :: mask_east !< flag to mask vertical wall east of the grid point
273	REAL(wp) :: mask_north !< flag to mask vertical wall north of the grid point
274	REAL(wp) :: mask_south !< flag to mask vertical wall south of the grid point
275	REAL(wp) :: mask_west !< flag to mask vertical wall west of the grid point
276
277
278	DO k = nzb+1, nzt-1
279	!
280	!-- Predetermine flag to mask topography and wall-bounded grid points.
281	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
282	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
283	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
284	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
285	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
286	!
287	!-- Interpolate eddy diffusivities on staggered gridpoints
288	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
289	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
290	kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
291	kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
292
293	tend(k,j,i) = tend(k,j,i) &
294	+ ( mask_east * kmxp * ( &
295	( w(k,j,i+1) - w(k,j,i) ) * ddx &
296	+ ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
297	) &
298	- mask_west * kmxm * ( &
299	( w(k,j,i) - w(k,j,i-1) ) * ddx &
300	+ ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
301	) &
302	) * ddx * flag &
303	+ ( mask_north * kmyp * ( &
304	( w(k,j+1,i) - w(k,j,i) ) * ddy &
305	+ ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
306	) &
307	- mask_south * kmym * ( &
308	( w(k,j,i) - w(k,j-1,i) ) * ddy &
309	+ ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
310	) &
311	) * ddy * flag &
312	+ 2.0_wp * ( &
313	km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
314	* rho_air(k+1) &
315	- km(k,j,i) * ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
316	* rho_air(k) &
317	) * ddzu(k+1) * drho_air_zw(k) * flag
318	ENDDO
319	!
320	!-- Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. Note, in
321	!-- the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
322	!-- vertical natural surfaces so far.
323	!-- Default-type surfaces
324	DO l = 0, 1
325	surf_s = surf_def_v(l)%start_index(j,i)
326	surf_e = surf_def_v(l)%end_index(j,i)
327	DO m = surf_s, surf_e
328	k = surf_def_v(l)%k(m)
329	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
330	ENDDO
331	ENDDO
332	!
333	!-- Natural-type surfaces
334	DO l = 0, 1
335	surf_s = surf_lsm_v(l)%start_index(j,i)
336	surf_e = surf_lsm_v(l)%end_index(j,i)
337	DO m = surf_s, surf_e
338	k = surf_lsm_v(l)%k(m)
339	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
340	ENDDO
341	ENDDO
342	!
343	!-- Urban-type surfaces
344	DO l = 0, 1
345	surf_s = surf_usm_v(l)%start_index(j,i)
346	surf_e = surf_usm_v(l)%end_index(j,i)
347	DO m = surf_s, surf_e
348	k = surf_usm_v(l)%k(m)
349	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
350	ENDDO
351	ENDDO
352	!
353	!-- Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surfaces.
354	!-- Default-type surfaces
355	DO l = 2, 3
356	surf_s = surf_def_v(l)%start_index(j,i)
357	surf_e = surf_def_v(l)%end_index(j,i)
358	DO m = surf_s, surf_e
359	k = surf_def_v(l)%k(m)
360	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
361	ENDDO
362	ENDDO
363	!
364	!-- Natural-type surfaces
365	DO l = 2, 3
366	surf_s = surf_lsm_v(l)%start_index(j,i)
367	surf_e = surf_lsm_v(l)%end_index(j,i)
368	DO m = surf_s, surf_e
369	k = surf_lsm_v(l)%k(m)
370	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
371	ENDDO
372	ENDDO
373	!
374	!-- Urban-type surfaces
375	DO l = 2, 3
376	surf_s = surf_usm_v(l)%start_index(j,i)
377	surf_e = surf_usm_v(l)%end_index(j,i)
378	DO m = surf_s, surf_e
379	k = surf_usm_v(l)%k(m)
380	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
381	ENDDO
382	ENDDO
383
384
385	END SUBROUTINE diffusion_w_ij
386
387	END MODULE diffusion_w_mod

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

Context Navigation

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

Download in other formats: