Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_v.f90 @ 4180

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 23.4 KB

Line
1	!> @file diffusion_v.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_v.f90 4180 2019-08-21 14:37:54Z scharf $
27	! OpenACC port for SPEC
28	!
29	!
30	! Description:
31	! ------------
32	!> Diffusion term of the v-component
33	!------------------------------------------------------------------------------!
34	MODULE diffusion_v_mod
35
36
37	PRIVATE
38	PUBLIC diffusion_v
39
40	INTERFACE diffusion_v
41	MODULE PROCEDURE diffusion_v
42	MODULE PROCEDURE diffusion_v_ij
43	END INTERFACE diffusion_v
44
45	CONTAINS
46
47
48	!------------------------------------------------------------------------------!
49	! Description:
50	! ------------
51	!> Call for all grid points
52	!------------------------------------------------------------------------------!
53	SUBROUTINE diffusion_v
54
55	USE arrays_3d, &
56	ONLY: ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
57
58	USE control_parameters, &
59	ONLY: constant_top_momentumflux, use_surface_fluxes, &
60	use_top_fluxes
61
62	USE grid_variables, &
63	ONLY: ddx, ddy, ddy2
64
65	USE indices, &
66	ONLY: nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_0
67
68	USE kinds
69
70	USE surface_mod, &
71	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
72	surf_usm_v
73
74	IMPLICIT NONE
75
76	INTEGER(iwp) :: i !< running index x direction
77	INTEGER(iwp) :: j !< running index y direction
78	INTEGER(iwp) :: k !< running index z direction
79	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
80	INTEGER(iwp) :: m !< running index surface elements
81	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
82	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
83
84	REAL(wp) :: flag !< flag to mask topography grid points
85	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the v-gridbox - interpolated onto xu-yv grid
86	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the v-gridbox - interpolated onto xu-yv grid
87	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid
88	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto yv-zw grid
89	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
90	REAL(wp) :: mask_east !< flag to mask vertical surface south of the grid point
91	REAL(wp) :: mask_west !< flag to mask vertical surface north of the grid point
92	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
93
94	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
95	!$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmzm, kmzp) &
96	!$ACC PRIVATE(mask_bottom, mask_east, mask_west, mask_top) &
97	!$ACC PRESENT(wall_flags_0, km) &
98	!$ACC PRESENT(u, v, w) &
99	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
100	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(2:3)) &
101	!$ACC PRESENT(surf_lsm_h, surf_lsm_v(2:3)) &
102	!$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
103	!$ACC PRESENT(tend)
104	DO i = nxl, nxr
105	DO j = nysv, nyn
106	!
107	!-- Compute horizontal diffusion
108	DO k = nzb+1, nzt
109
110	!
111	!-- Predetermine flag to mask topography and wall-bounded grid points.
112	!-- It is sufficient to masked only east- and west-facing surfaces, which
113	!-- need special treatment for the v-component.
114	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) )
115	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 2 ) )
116	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 2 ) )
117	!
118	!-- Interpolate eddy diffusivities on staggered gridpoints
119	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
120	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
121
122	tend(k,j,i) = tend(k,j,i) + ( &
123	mask_east * kmxp * ( &
124	( v(k,j,i+1) - v(k,j,i) ) * ddx &
125	+ ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
126	) &
127	- mask_west * kmxm * ( &
128	( v(k,j,i) - v(k,j,i-1) ) * ddx &
129	+ ( u(k,j,i) - u(k,j-1,i) ) * ddy &
130	) &
131	) * ddx * flag &
132	+ 2.0_wp * ( &
133	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
134	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
135	) * ddy2 * flag
136
137	ENDDO
138
139	!
140	!-- Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
141	!-- surfaces. Note, in the the flat case, loops won't be entered as
142	!-- start_index > end_index. Furtermore, note, no vertical natural surfaces
143	!-- so far.
144	!-- Default-type surfaces
145	DO l = 2, 3
146	surf_s = surf_def_v(l)%start_index(j,i)
147	surf_e = surf_def_v(l)%end_index(j,i)
148	DO m = surf_s, surf_e
149	k = surf_def_v(l)%k(m)
150	tend(k,j,i) = tend(k,j,i) + &
151	surf_def_v(l)%mom_flux_uv(m) * ddx
152	ENDDO
153	ENDDO
154	!
155	!-- Natural-type surfaces
156	DO l = 2, 3
157	surf_s = surf_lsm_v(l)%start_index(j,i)
158	surf_e = surf_lsm_v(l)%end_index(j,i)
159	DO m = surf_s, surf_e
160	k = surf_lsm_v(l)%k(m)
161	tend(k,j,i) = tend(k,j,i) + &
162	surf_lsm_v(l)%mom_flux_uv(m) * ddx
163	ENDDO
164	ENDDO
165	!
166	!-- Urban-type surfaces
167	DO l = 2, 3
168	surf_s = surf_usm_v(l)%start_index(j,i)
169	surf_e = surf_usm_v(l)%end_index(j,i)
170	DO m = surf_s, surf_e
171	k = surf_usm_v(l)%k(m)
172	tend(k,j,i) = tend(k,j,i) + &
173	surf_usm_v(l)%mom_flux_uv(m) * ddx
174	ENDDO
175	ENDDO
176	!
177	!-- Compute vertical diffusion. In case of simulating a surface layer,
178	!-- respective grid diffusive fluxes are masked (flag 10) within this
179	!-- loop, and added further below, else, simple gradient approach is
180	!-- applied. Model top is also mask if top-momentum flux is given.
181	DO k = nzb+1, nzt
182	!
183	!-- Determine flags to mask topography below and above. Flag 2 is
184	!-- used to mask topography in general, while flag 8 implies also
185	!-- information about use_surface_fluxes. Flag 9 is used to control
186	!-- momentum flux at model top.
187	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
188	BTEST( wall_flags_0(k-1,j,i), 8 ) )
189	mask_top = MERGE( 1.0_wp, 0.0_wp, &
190	BTEST( wall_flags_0(k+1,j,i), 8 ) ) * &
191	MERGE( 1.0_wp, 0.0_wp, &
192	BTEST( wall_flags_0(k+1,j,i), 9 ) )
193	flag = MERGE( 1.0_wp, 0.0_wp, &
194	BTEST( wall_flags_0(k,j,i), 2 ) )
195	!
196	!-- Interpolate eddy diffusivities on staggered gridpoints
197	kmzp = 0.25_wp * &
198	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
199	kmzm = 0.25_wp * &
200	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
201
202	tend(k,j,i) = tend(k,j,i) &
203	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
204	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
205	& ) * rho_air_zw(k) * mask_top &
206	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
207	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
208	& ) * rho_air_zw(k-1) * mask_bottom &
209	& ) * ddzw(k) * drho_air(k) * flag
210	ENDDO
211
212	!
213	!-- Vertical diffusion at the first grid point above the surface,
214	!-- if the momentum flux at the bottom is given by the Prandtl law
215	!-- or if it is prescribed by the user.
216	!-- Difference quotient of the momentum flux is not formed over
217	!-- half of the grid spacing (2.0*ddzw(k)) any more, since the
218	!-- comparison with other (LES) models showed that the values of
219	!-- the momentum flux becomes too large in this case.
220	IF ( use_surface_fluxes ) THEN
221	!
222	!-- Default-type surfaces, upward-facing
223	surf_s = surf_def_h(0)%start_index(j,i)
224	surf_e = surf_def_h(0)%end_index(j,i)
225	DO m = surf_s, surf_e
226	k = surf_def_h(0)%k(m)
227
228	tend(k,j,i) = tend(k,j,i) &
229	+ ( - ( - surf_def_h(0)%vsws(m) ) &
230	) * ddzw(k) * drho_air(k)
231	ENDDO
232	!
233	!-- Default-type surfaces, dowward-facing
234	surf_s = surf_def_h(1)%start_index(j,i)
235	surf_e = surf_def_h(1)%end_index(j,i)
236	DO m = surf_s, surf_e
237	k = surf_def_h(1)%k(m)
238
239	tend(k,j,i) = tend(k,j,i) &
240	+ ( - surf_def_h(1)%vsws(m) &
241	) * ddzw(k) * drho_air(k)
242	ENDDO
243	!
244	!-- Natural-type surfaces, upward-facing
245	surf_s = surf_lsm_h%start_index(j,i)
246	surf_e = surf_lsm_h%end_index(j,i)
247	DO m = surf_s, surf_e
248	k = surf_lsm_h%k(m)
249
250	tend(k,j,i) = tend(k,j,i) &
251	+ ( - ( - surf_lsm_h%vsws(m) ) &
252	) * ddzw(k) * drho_air(k)
253
254	ENDDO
255	!
256	!-- Urban-type surfaces, upward-facing
257	surf_s = surf_usm_h%start_index(j,i)
258	surf_e = surf_usm_h%end_index(j,i)
259	DO m = surf_s, surf_e
260	k = surf_usm_h%k(m)
261
262	tend(k,j,i) = tend(k,j,i) &
263	+ ( - ( - surf_usm_h%vsws(m) ) &
264	) * ddzw(k) * drho_air(k)
265
266	ENDDO
267	ENDIF
268	!
269	!-- Add momentum flux at model top
270	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
271	surf_s = surf_def_h(2)%start_index(j,i)
272	surf_e = surf_def_h(2)%end_index(j,i)
273	DO m = surf_s, surf_e
274
275	k = surf_def_h(2)%k(m)
276
277	tend(k,j,i) = tend(k,j,i) &
278	+ ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
279	ENDDO
280	ENDIF
281
282	ENDDO
283	ENDDO
284
285	END SUBROUTINE diffusion_v
286
287
288	!------------------------------------------------------------------------------!
289	! Description:
290	! ------------
291	!> Call for grid point i,j
292	!------------------------------------------------------------------------------!
293	SUBROUTINE diffusion_v_ij( i, j )
294
295	USE arrays_3d, &
296	ONLY: ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
297
298	USE control_parameters, &
299	ONLY: constant_top_momentumflux, use_surface_fluxes, &
300	use_top_fluxes
301
302	USE grid_variables, &
303	ONLY: ddx, ddy, ddy2
304
305	USE indices, &
306	ONLY: nzb, nzt, wall_flags_0
307
308	USE kinds
309
310	USE surface_mod, &
311	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
312	surf_usm_v
313
314	IMPLICIT NONE
315
316
317	INTEGER(iwp) :: i !< running index x direction
318	INTEGER(iwp) :: j !< running index y direction
319	INTEGER(iwp) :: k !< running index z direction
320	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
321	INTEGER(iwp) :: m !< running index surface elements
322	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
323	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
324
325	REAL(wp) :: flag !< flag to mask topography grid points
326	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the v-gridbox - interpolated onto xu-yv grid
327	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the v-gridbox - interpolated onto xu-yv grid
328	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
329	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
330	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
331	REAL(wp) :: mask_east !< flag to mask vertical surface south of the grid point
332	REAL(wp) :: mask_west !< flag to mask vertical surface north of the grid point
333	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
334
335	!
336	!-- Compute horizontal diffusion
337	DO k = nzb+1, nzt
338	!
339	!-- Predetermine flag to mask topography and wall-bounded grid points.
340	!-- It is sufficient to masked only east- and west-facing surfaces, which
341	!-- need special treatment for the v-component.
342	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) )
343	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 2 ) )
344	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 2 ) )
345	!
346	!-- Interpolate eddy diffusivities on staggered gridpoints
347	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
348	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
349
350	tend(k,j,i) = tend(k,j,i) + ( &
351	mask_east * kmxp * ( &
352	( v(k,j,i+1) - v(k,j,i) ) * ddx &
353	+ ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
354	) &
355	- mask_west * kmxm * ( &
356	( v(k,j,i) - v(k,j,i-1) ) * ddx &
357	+ ( u(k,j,i) - u(k,j-1,i) ) * ddy &
358	) &
359	) * ddx * flag &
360	+ 2.0_wp * ( &
361	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
362	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
363	) * ddy2 * flag
364	ENDDO
365
366	!
367	!-- Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
368	!-- surfaces. Note, in the the flat case, loops won't be entered as
369	!-- start_index > end_index. Furtermore, note, no vertical natural surfaces
370	!-- so far.
371	!-- Default-type surfaces
372	DO l = 2, 3
373	surf_s = surf_def_v(l)%start_index(j,i)
374	surf_e = surf_def_v(l)%end_index(j,i)
375	DO m = surf_s, surf_e
376	k = surf_def_v(l)%k(m)
377	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(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) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
388	ENDDO
389	ENDDO
390	!
391	!-- Urban-type surfaces
392	DO l = 2, 3
393	surf_s = surf_usm_v(l)%start_index(j,i)
394	surf_e = surf_usm_v(l)%end_index(j,i)
395	DO m = surf_s, surf_e
396	k = surf_usm_v(l)%k(m)
397	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
398	ENDDO
399	ENDDO
400	!
401	!-- Compute vertical diffusion. In case of simulating a surface layer,
402	!-- respective grid diffusive fluxes are masked (flag 8) within this
403	!-- loop, and added further below, else, simple gradient approach is
404	!-- applied. Model top is also mask if top-momentum flux is given.
405	DO k = nzb+1, nzt
406	!
407	!-- Determine flags to mask topography below and above. Flag 2 is
408	!-- used to mask topography in general, while flag 10 implies also
409	!-- information about use_surface_fluxes. Flag 9 is used to control
410	!-- momentum flux at model top.
411	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
412	BTEST( wall_flags_0(k-1,j,i), 8 ) )
413	mask_top = MERGE( 1.0_wp, 0.0_wp, &
414	BTEST( wall_flags_0(k+1,j,i), 8 ) ) * &
415	MERGE( 1.0_wp, 0.0_wp, &
416	BTEST( wall_flags_0(k+1,j,i), 9 ) )
417	flag = MERGE( 1.0_wp, 0.0_wp, &
418	BTEST( wall_flags_0(k,j,i), 2 ) )
419	!
420	!-- Interpolate eddy diffusivities on staggered gridpoints
421	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
422	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
423
424	tend(k,j,i) = tend(k,j,i) &
425	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
426	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
427	& ) * rho_air_zw(k) * mask_top &
428	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
429	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
430	& ) * rho_air_zw(k-1) * mask_bottom &
431	& ) * ddzw(k) * drho_air(k) * flag
432	ENDDO
433
434	!
435	!-- Vertical diffusion at the first grid point above the surface, if the
436	!-- momentum flux at the bottom is given by the Prandtl law or if it is
437	!-- prescribed by the user.
438	!-- Difference quotient of the momentum flux is not formed over half of
439	!-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with
440	!-- other (LES) models showed that the values of the momentum flux becomes
441	!-- too large in this case.
442	IF ( use_surface_fluxes ) THEN
443	!
444	!-- Default-type surfaces, upward-facing
445	surf_s = surf_def_h(0)%start_index(j,i)
446	surf_e = surf_def_h(0)%end_index(j,i)
447	DO m = surf_s, surf_e
448	k = surf_def_h(0)%k(m)
449
450	tend(k,j,i) = tend(k,j,i) &
451	+ ( - ( - surf_def_h(0)%vsws(m) ) &
452	) * ddzw(k) * drho_air(k)
453	ENDDO
454	!
455	!-- Default-type surfaces, dowward-facing
456	surf_s = surf_def_h(1)%start_index(j,i)
457	surf_e = surf_def_h(1)%end_index(j,i)
458	DO m = surf_s, surf_e
459	k = surf_def_h(1)%k(m)
460
461	tend(k,j,i) = tend(k,j,i) &
462	+ ( - surf_def_h(1)%vsws(m) &
463	) * ddzw(k) * drho_air(k)
464	ENDDO
465	!
466	!-- Natural-type surfaces, upward-facing
467	surf_s = surf_lsm_h%start_index(j,i)
468	surf_e = surf_lsm_h%end_index(j,i)
469	DO m = surf_s, surf_e
470	k = surf_lsm_h%k(m)
471
472	tend(k,j,i) = tend(k,j,i) &
473	+ ( - ( - surf_lsm_h%vsws(m) ) &
474	) * ddzw(k) * drho_air(k)
475
476	ENDDO
477	!
478	!-- Urban-type surfaces, upward-facing
479	surf_s = surf_usm_h%start_index(j,i)
480	surf_e = surf_usm_h%end_index(j,i)
481	DO m = surf_s, surf_e
482	k = surf_usm_h%k(m)
483
484	tend(k,j,i) = tend(k,j,i) &
485	+ ( - ( - surf_usm_h%vsws(m) ) &
486	) * ddzw(k) * drho_air(k)
487
488	ENDDO
489	ENDIF
490	!
491	!-- Add momentum flux at model top
492	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
493	surf_s = surf_def_h(2)%start_index(j,i)
494	surf_e = surf_def_h(2)%end_index(j,i)
495	DO m = surf_s, surf_e
496
497	k = surf_def_h(2)%k(m)
498
499	tend(k,j,i) = tend(k,j,i) &
500	+ ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
501	ENDDO
502	ENDIF
503
504	END SUBROUTINE diffusion_v_ij
505
506	END MODULE diffusion_v_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |