Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |