Home

Context Navigation

source: palm/trunk/SOURCE/init_grid.f90 @ 2876

Last change on this file since 2876 was 2867, checked in by suehring, 6 years ago
Revise mapping of 3D-buildings onto orography; further bugfix concering call of user_init (see commit -r 2865)
Property svn:keywords set to `Id`
File size: 99.9 KB

Line
1	!> @file init_grid.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-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: init_grid.f90 2867 2018-03-09 09:40:23Z knoop $
27	! Revise mapping of 3D buildings onto onto orography.
28	!
29	! 2823 2018-02-20 15:31:45Z Giersch
30	! Set boundary conditions for 3D topography in case of non-cyclic boundary
31	! conditions
32	!
33	! 2796 2018-02-08 12:25:39Z suehring
34	! Bugfix in 3D building initialization
35	!
36	! 2747 2018-01-15 12:44:17Z suehring
37	! Bugfix, topography height is rounded to the nearest discrete grid level
38	!
39	! 2718 2018-01-02 08:49:38Z maronga
40	! Corrected "Former revisions" section
41	!
42	! 2701 2017-12-15 15:40:50Z suehring
43	! Changes from last commit documented
44	!
45	! 2698 2017-12-14 18:46:24Z suehring
46	! Bugfix in get_topography_top_index
47	!
48	! 2696 2017-12-14 17:12:51Z kanani
49	! Change in file header (GPL part)
50	! Revised topography input
51	! Set nzb_max not for the entire nest domain, only for boundary PEs
52	! Re-organize routine, split-up into several subroutines
53	! Modularize poismg_noopt
54	! Remove setting bit 26, 27, 28 in wall_flags_0, indicating former '_outer'
55	! arrays (not required any more).
56	! Bugfix in generic tunnel setup (MS)
57	!
58	! 2550 2017-10-16 17:12:01Z boeske
59	! Set lateral boundary conditions for topography on all three ghost layers
60	!
61	! 2478 2017-09-18 13:37:24Z suehring
62	! Bugfix, correct flag for use_top
63	!
64	! 2365 2017-08-21 14:59:59Z kanani
65	! Vertical nesting implemented (SadiqHuq)
66	!
67	! 2319 2017-07-20 17:33:17Z suehring
68	! Remove print statements
69	!
70	! 2318 2017-07-20 17:27:44Z suehring
71	! Get topography top index via Function call
72	!
73	! 2317 2017-07-20 17:27:19Z suehring
74	! Bugfixes in reading 3D topography from file
75	!
76	! 2274 2017-06-09 13:27:48Z Giersch
77	! Changed error messages
78	!
79	! 2233 2017-05-30 18:08:54Z suehring
80	!
81	! 2232 2017-05-30 17:47:52Z suehring
82	! - Adjustments according to new topography representation
83	! - Bugfix: Move determination of nzb_max behind topography modification in
84	! cell-edge case
85	! - Get rid off global arrays required for topography output
86	! - Enable topography input via netcdf
87	! - Generic tunnel set-up added
88	!
89	! 2200 2017-04-11 11:37:51Z suehring
90	! monotonic_adjustment removed
91	!
92	! 2169 2017-03-06 18:16:35Z suehring
93	! Bugfix, move setting for topography grid convention to init_grid, else, if no
94	! value is set, the simulation may abort in case of restarts
95	!
96	! 2128 2017-01-23 15:00:03Z suehring
97	! Bugfix in setting topography from file in case of ocean simulations
98	!
99	! 2088 2016-12-19 16:30:25Z suehring
100	! Bugfix in generic topography in case of ocean simulations
101	!
102	! 2037 2016-10-26 11:15:40Z knoop
103	! Anelastic approximation implemented
104	!
105	! 2021 2016-10-07 14:08:57Z suehring
106	! Bugfix: setting Neumann boundary conditions for topography required for
107	! topography flags in multigrid_noopt solver
108	!
109	! 2000 2016-08-20 18:09:15Z knoop
110	! Forced header and separation lines into 80 columns
111	!
112	! 1994 2016-08-15 09:52:21Z suehring
113	! Bugfix in definition of generic topography
114	!
115	! 1982 2016-08-01 11:04:48Z suehring
116	! Bugfix concering consistency check for topography
117	!
118	! 1968 2016-07-18 12:01:49Z suehring
119	! Changed: PE-wise reading of topography file in order to avoid global definition
120	! of arrays nzb_local and nzb_tmp. Thereby, topography definition for single
121	! buildings and street canyons has changed, as well as flag setting for
122	! multigrid scheme.
123	!
124	! Bugfix in checking l_grid anisotropy.
125	! Simplify initial computation of lwall and vertical_influence, i.e. remove
126	! nzb_s_inner as it is still zero at this point.
127	!
128	! 1942 2016-06-14 12:18:18Z suehring
129	! Topography filter implemented to fill holes resolved by only one grid point.
130	! Initialization of flags for ws-scheme moved to advec_ws.
131	!
132	! 1931 2016-06-10 12:06:59Z suehring
133	! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
134	!
135	! 1910 2016-05-26 06:49:46Z raasch
136	! Bugfix: if topography is read from file, Neumann conditions are used for the
137	! nzb_local array (instead of cyclic conditions) in case that non-cyclic
138	! boundary conditions are switched on for the run
139	!
140	! 1902 2016-05-09 11:18:56Z suehring
141	! Set topography flags for multigrid solver only (not for multigrid_fast)
142	!
143	! 1886 2016-04-21 11:20:47Z suehring
144	! Bugfix: setting advection flags near walls
145	! reformulated index values for nzb_v_inner
146	! variable discriptions added in declaration block
147	!
148	! 1845 2016-04-08 08:29:13Z raasch
149	! nzb_2d removed
150	!
151	! 1804 2016-04-05 16:30:18Z maronga
152	! Removed code for parameter file check (__check)
153	!
154	! 1779 2016-03-03 08:01:28Z raasch
155	! coupling_char is trimmed at every place it occurs, because it can have
156	! different length now
157	!
158	! 1762 2016-02-25 12:31:13Z hellstea
159	! Introduction of nested domain feature
160	!
161	! 1743 2016-01-13 10:23:51Z raasch
162	! Bugfix for calculation of nzb_s_outer and nzb_u_outer at north boundary of
163	! total domain
164	!
165	! 1691 2015-10-26 16:17:44Z maronga
166	! Renamed prandtl_layer to constant_flux_layer.
167	!
168	! 1682 2015-10-07 23:56:08Z knoop
169	! Code annotations made doxygen readable
170	!
171	! 1677 2015-10-02 13:25:23Z boeske
172	! Bugfix: Ghost points are included in wall_flags_0 and wall_flags_00
173	!
174	! 1675 2015-10-02 08:28:59Z gronemeier
175	! Bugfix: Definition of topography grid levels
176	!
177	! 1660 2015-09-21 08:15:16Z gronemeier
178	! Bugfix: Definition of topography grid levels if vertical grid stretching
179	! starts below the maximum topography height.
180	!
181	! 1580 2015-04-10 13:43:49Z suehring
182	! Bugfix: setting flags for 5th order scheme near buildings
183	!
184	! 1575 2015-03-27 09:56:27Z raasch
185	! adjustments for psolver-queries
186	!
187	! 1557 2015-03-05 16:43:04Z suehring
188	! Adjustment for monotoinic limiter
189	!
190	! 1418 2014-06-06 13:05:08Z fricke
191	! Bugfix: Change if-condition for stretched grid in the ocean, with the old
192	! condition and a negative value for dz_stretch_level the condition
193	! was always true for the whole model domain
194	!
195	! 1409 2014-05-23 12:11:32Z suehring
196	! Bugfix: set wall_flags_0 at inflow and outflow boundary also for i <= nxlu
197	! j <= nysv
198	!
199	! 1353 2014-04-08 15:21:23Z heinze
200	! REAL constants provided with KIND-attribute
201	!
202	! 1322 2014-03-20 16:38:49Z raasch
203	! REAL constants defined as wp-kind
204	!
205	! 1320 2014-03-20 08:40:49Z raasch
206	! ONLY-attribute added to USE-statements,
207	! kind-parameters added to all INTEGER and REAL declaration statements,
208	! kinds are defined in new module kinds,
209	! revision history before 2012 removed,
210	! comment fields (!:) to be used for variable explanations added to
211	! all variable declaration statements
212	!
213	! 1221 2013-09-10 08:59:13Z raasch
214	! wall_flags_00 introduced to hold bits 32-63,
215	! additional 3D-flag arrays for replacing the 2D-index array nzb_s_inner in
216	! loops optimized for openACC (pres + flow_statistics)
217	!
218	! 1092 2013-02-02 11:24:22Z raasch
219	! unused variables removed
220	!
221	! 1069 2012-11-28 16:18:43Z maronga
222	! bugfix: added coupling_char to TOPOGRAPHY_DATA to allow topography in the
223	! ocean model in case of coupled runs
224	!
225	! 1036 2012-10-22 13:43:42Z raasch
226	! code put under GPL (PALM 3.9)
227	!
228	! 1015 2012-09-27 09:23:24Z raasch
229	! lower index for calculating wall_flags_0 set to nzb_w_inner instead of
230	! nzb_w_inner+1
231	!
232	! 996 2012-09-07 10:41:47Z raasch
233	! little reformatting
234	!
235	! 978 2012-08-09 08:28:32Z fricke
236	! Bugfix: nzb_max is set to nzt at non-cyclic lateral boundaries
237	! Bugfix: Set wall_flags_0 for inflow boundary
238	!
239	! 927 2012-06-06 19:15:04Z raasch
240	! Wall flags are not set for multigrid method in case of masking method
241	!
242	! 864 2012-03-27 15:10:33Z gryschka
243	! In case of ocean and Dirichlet bottom bc for u and v dzu_mg and ddzu_pres
244	! were not correctly defined for k=1.
245	!
246	! 861 2012-03-26 14:18:34Z suehring
247	! Set wall_flags_0. The array is needed for degradation in ws-scheme near walls,
248	! inflow and outflow boundaries as well as near the bottom and the top of the
249	! model domain.!
250	! Initialization of nzb_s_inner and nzb_w_inner.
251	! gls has to be at least nbgp to do not exceed the array bounds of nzb_local
252	! while setting wall_flags_0
253	!
254	! 843 2012-02-29 15:16:21Z gryschka
255	! In case of ocean and dirichlet bc for u and v at the bottom
256	! the first u-level ist defined at same height as the first w-level
257	!
258	! 818 2012-02-08 16:11:23Z maronga
259	! Bugfix: topo_height is only required if topography is used. It is thus now
260	! allocated in the topography branch
261	!
262	! 809 2012-01-30 13:32:58Z maronga
263	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
264	!
265	! 807 2012-01-25 11:53:51Z maronga
266	! New cpp directive "__check" implemented which is used by check_namelist_files
267	!
268	! Revision 1.1 1997/08/11 06:17:45 raasch
269	! Initial revision (Testversion)
270	!
271	!
272	! Description:
273	! -----------------------------------------------------------------------------!
274	!> Creating grid depending constants
275	!> @todo: Move initialization mixing length
276	!> @todo: Rearrange topo flag list
277	!> @todo: reference 3D buildings on top of orography is not tested and may need
278	!> further improvement for steep slopes
279	!> @todo: Use more advanced setting of building type at filled holes
280	!------------------------------------------------------------------------------!
281	SUBROUTINE init_grid
282
283	USE advec_ws, &
284	ONLY: ws_init_flags
285
286	USE arrays_3d, &
287	ONLY: dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, zu, zw
288
289	USE control_parameters, &
290	ONLY: bc_lr_cyc, bc_ns_cyc, building_height, building_length_x, &
291	building_length_y, building_wall_left, building_wall_south, &
292	canyon_height, canyon_wall_left, canyon_wall_south, &
293	canyon_width_x, canyon_width_y, constant_flux_layer, &
294	dp_level_ind_b, dz, dz_max, dz_stretch_factor, &
295	dz_stretch_level, dz_stretch_level_index, grid_level, &
296	force_bound_l, force_bound_r, force_bound_n, force_bound_s, &
297	ibc_uv_b, inflow_l, inflow_n, inflow_r, inflow_s, &
298	masking_method, maximum_grid_level, message_string, &
299	momentum_advec, nest_domain, nest_bound_l, nest_bound_n, &
300	nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n, &
301	outflow_r, outflow_s, psolver, scalar_advec, topography, &
302	topography_grid_convention, tunnel_height, tunnel_length, &
303	tunnel_width_x, tunnel_width_y, tunnel_wall_depth, &
304	use_surface_fluxes, use_top_fluxes, wall_adjustment_factor
305
306	USE grid_variables, &
307	ONLY: ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner
308
309	USE indices, &
310	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, &
311	nzb, nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, &
312	nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner, &
313	nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner, &
314	nzb_w_outer, nzt
315
316	USE kinds
317
318	USE pegrid
319
320	USE poismg_noopt_mod
321
322	USE surface_mod, &
323	ONLY: get_topography_top_index, get_topography_top_index_ji, init_bc
324
325	USE vertical_nesting_mod, &
326	ONLY: vnested, vnest_init_grid
327
328	IMPLICIT NONE
329
330	INTEGER(iwp) :: i !< index variable along x
331	INTEGER(iwp) :: j !< index variable along y
332	INTEGER(iwp) :: k !< index variable along z
333	INTEGER(iwp) :: k_top !< topography top index on local PE
334	INTEGER(iwp) :: l !< loop variable
335	INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height
336	INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height
337
338	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center
339	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u- and v-grid
340
341	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo !< input array for 3D topography and dummy array for setting "outer"-flags
342
343	REAL(wp) :: dz_stretched !< stretched vertical grid spacing
344
345
346	!
347	!-- Calculation of horizontal array bounds including ghost layers
348	nxlg = nxl - nbgp
349	nxrg = nxr + nbgp
350	nysg = nys - nbgp
351	nyng = nyn + nbgp
352
353	!
354	!-- Allocate grid arrays
355	ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
356	dzw(1:nzt+1), zu(nzb:nzt+1), zw(nzb:nzt+1) )
357
358	!
359	!-- Compute height of u-levels from constant grid length and dz stretch factors
360	IF ( dz == -1.0_wp ) THEN
361	message_string = 'missing dz'
362	CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
363	ELSEIF ( dz <= 0.0_wp ) THEN
364	WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
365	CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
366	ENDIF
367
368	!
369	!-- Define the vertical grid levels
370	IF ( .NOT. ocean ) THEN
371	!
372	!-- Grid for atmosphere with surface at z=0 (k=0, w-grid).
373	!-- The second u-level (k=1) corresponds to the top of the
374	!-- Prandtl-layer.
375
376	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
377	zu(0) = 0.0_wp
378	! zu(0) = - dz * 0.5_wp
379	ELSE
380	zu(0) = - dz * 0.5_wp
381	ENDIF
382	zu(1) = dz * 0.5_wp
383
384	dz_stretch_level_index = nzt+1
385	dz_stretched = dz
386	DO k = 2, nzt+1
387	IF ( dz_stretch_level <= zu(k-1) .AND. dz_stretched < dz_max ) THEN
388	dz_stretched = dz_stretched * dz_stretch_factor
389	dz_stretched = MIN( dz_stretched, dz_max )
390	IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
391	ENDIF
392	zu(k) = zu(k-1) + dz_stretched
393	ENDDO
394
395	!
396	!-- Compute the w-levels. They are always staggered half-way between the
397	!-- corresponding u-levels. In case of dirichlet bc for u and v at the
398	!-- ground the first u- and w-level (k=0) are defined at same height (z=0).
399	!-- The top w-level is extrapolated linearly.
400	zw(0) = 0.0_wp
401	DO k = 1, nzt
402	zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
403	ENDDO
404	zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
405
406	ELSE
407	!
408	!-- Grid for ocean with free water surface is at k=nzt (w-grid).
409	!-- In case of neumann bc at the ground the first first u-level (k=0) lies
410	!-- below the first w-level (k=0). In case of dirichlet bc the first u- and
411	!-- w-level are defined at same height, but staggered from the second level.
412	!-- The second u-level (k=1) corresponds to the top of the Prandtl-layer.
413	zu(nzt+1) = dz * 0.5_wp
414	zu(nzt) = - dz * 0.5_wp
415
416	dz_stretch_level_index = 0
417	dz_stretched = dz
418	DO k = nzt-1, 0, -1
419	!
420	!-- The default value of dz_stretch_level is positive, thus the first
421	!-- condition is always true. Hence, the second condition is necessary.
422	IF ( dz_stretch_level >= zu(k+1) .AND. dz_stretch_level <= 0.0 &
423	.AND. dz_stretched < dz_max ) THEN
424	dz_stretched = dz_stretched * dz_stretch_factor
425	dz_stretched = MIN( dz_stretched, dz_max )
426	IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
427	ENDIF
428	zu(k) = zu(k+1) - dz_stretched
429	ENDDO
430
431	!
432	!-- Compute the w-levels. They are always staggered half-way between the
433	!-- corresponding u-levels, except in case of dirichlet bc for u and v
434	!-- at the ground. In this case the first u- and w-level are defined at
435	!-- same height. The top w-level (nzt+1) is not used but set for
436	!-- consistency, since w and all scalar variables are defined up tp nzt+1.
437	zw(nzt+1) = dz
438	zw(nzt) = 0.0_wp
439	DO k = 0, nzt
440	zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
441	ENDDO
442
443	!
444	!-- In case of dirichlet bc for u and v the first u- and w-level are defined
445	!-- at same height.
446	IF ( ibc_uv_b == 0 ) THEN
447	zu(0) = zw(0)
448	ENDIF
449
450	ENDIF
451
452	!
453	!-- Compute grid lengths.
454	DO k = 1, nzt+1
455	dzu(k) = zu(k) - zu(k-1)
456	ddzu(k) = 1.0_wp / dzu(k)
457	dzw(k) = zw(k) - zw(k-1)
458	ddzw(k) = 1.0_wp / dzw(k)
459	ENDDO
460
461	DO k = 1, nzt
462	dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
463	ENDDO
464
465	!
466	!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
467	!-- everywhere. For the actual grid, the grid spacing at the lowest level
468	!-- is only dz/2, but should be dz. Therefore, an additional array
469	!-- containing with appropriate grid information is created for these
470	!-- solvers.
471	IF ( psolver(1:9) /= 'multigrid' ) THEN
472	ALLOCATE( ddzu_pres(1:nzt+1) )
473	ddzu_pres = ddzu
474	ddzu_pres(1) = ddzu_pres(2) ! change for lowest level
475	ENDIF
476
477	!
478	!-- Compute the reciprocal values of the horizontal grid lengths.
479	ddx = 1.0_wp / dx
480	ddy = 1.0_wp / dy
481	dx2 = dx * dx
482	dy2 = dy * dy
483	ddx2 = 1.0_wp / dx2
484	ddy2 = 1.0_wp / dy2
485
486	!
487	!-- Allocate 3D array to set topography
488	ALLOCATE( topo(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
489	topo = 0
490	!
491	!-- Initialize topography by generic topography or read from topography from file.
492	CALL init_topo( topo )
493	!
494	!-- Set flags to mask topography on the grid.
495	CALL set_topo_flags( topo )
496	!
497	!-- Calculate wall flag arrays for the multigrid method.
498	!-- Please note, wall flags are only applied in the non-optimized version.
499	IF ( psolver == 'multigrid_noopt' ) CALL poismg_noopt_init
500
501	!
502	!-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e.
503	!-- to decrease the numerical stencil appropriately.
504	IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' ) &
505	CALL ws_init_flags
506
507	!
508	!-- Calculated grid-dependent as well as near-wall mixing length
509	CALL init_mixing_length
510
511	!
512	!-- Determine the maximum level of topography. It is used for
513	!-- steering the degradation of order of the applied advection scheme,
514	!-- as well in the lpm.
515	!-- In case of non-cyclic lateral boundaries, the order of the advection
516	!-- scheme has to be reduced up to nzt (required at the lateral boundaries).
517	k_top = 0
518	DO i = nxl, nxr
519	DO j = nys, nyn
520	DO k = nzb, nzt + 1
521	k_top = MAX( k_top, MERGE( k, 0, &
522	.NOT. BTEST( topo(k,j,i), 0 ) ) )
523	ENDDO
524	ENDDO
525	ENDDO
526	#if defined( __parallel )
527	CALL MPI_ALLREDUCE( k_top + 1, nzb_max, 1, MPI_INTEGER, & !is +1 really necessary here?
528	MPI_MAX, comm2d, ierr )
529	#else
530	nzb_max = k_top + 1
531	#endif
532	IF ( inflow_l .OR. outflow_l .OR. force_bound_l .OR. nest_bound_l .OR.&
533	inflow_r .OR. outflow_r .OR. force_bound_r .OR. nest_bound_r .OR.&
534	inflow_n .OR. outflow_n .OR. force_bound_n .OR. nest_bound_n .OR.&
535	inflow_s .OR. outflow_s .OR. force_bound_s .OR. nest_bound_s ) &
536	nzb_max = nzt
537	!
538	!-- Finally, if topography extents up to the model top, limit nzb_max to nzt.
539	nzb_max = MIN( nzb_max, nzt )
540
541	!
542	!-- Initialize boundary conditions via surface type
543	CALL init_bc
544
545	!
546	!-- Allocate and set topography height arrays required for data output
547	IF ( TRIM( topography ) /= 'flat' ) THEN
548	!
549	!-- Allocate and set the arrays containing the topography height
550	IF ( nxr == nx .AND. nyn /= ny ) THEN
551	ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn), &
552	zw_w_inner(nxl:nxr+1,nys:nyn) )
553	ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
554	ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1), &
555	zw_w_inner(nxl:nxr,nys:nyn+1) )
556	ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
557	ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1), &
558	zw_w_inner(nxl:nxr+1,nys:nyn+1) )
559	ELSE
560	ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn), &
561	zw_w_inner(nxl:nxr,nys:nyn) )
562	ENDIF
563
564	zu_s_inner = 0.0_wp
565	zw_w_inner = 0.0_wp
566	!
567	!-- Determine local topography height on scalar and w-grid. Note, setting
568	!-- lateral boundary values is not necessary, realized via wall_flags_0
569	!-- array. Further, please note that loop bounds are different from
570	!-- nxl to nxr and nys to nyn on south and right model boundary, hence,
571	!-- use intrinsic lbound and ubound functions to infer array bounds.
572	DO i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
573	DO j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
574	!
575	!-- Topography height on scalar grid. Therefore, determine index of
576	!-- upward-facing surface element on scalar grid.
577	zu_s_inner(i,j) = zu( get_topography_top_index_ji( j, i, 's' ) )
578	!
579	!-- Topography height on w grid. Therefore, determine index of
580	!-- upward-facing surface element on w grid.
581	zw_w_inner(i,j) = zw( get_topography_top_index_ji( j, i, 's' ) )
582	ENDDO
583	ENDDO
584	ENDIF
585
586	!
587	!-- In the following, calculate 2D index arrays. Note, these will be removed
588	!-- soon.
589	!-- Allocate outer and inner index arrays for topography and set
590	!-- defaults.
591	ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg), &
592	nzb_s_outer(nysg:nyng,nxlg:nxrg), &
593	nzb_u_inner(nysg:nyng,nxlg:nxrg), &
594	nzb_u_outer(nysg:nyng,nxlg:nxrg), &
595	nzb_v_inner(nysg:nyng,nxlg:nxrg), &
596	nzb_v_outer(nysg:nyng,nxlg:nxrg), &
597	nzb_w_inner(nysg:nyng,nxlg:nxrg), &
598	nzb_w_outer(nysg:nyng,nxlg:nxrg), &
599	nzb_diff_s_inner(nysg:nyng,nxlg:nxrg), &
600	nzb_diff_s_outer(nysg:nyng,nxlg:nxrg), &
601	nzb_local(nysg:nyng,nxlg:nxrg), &
602	nzb_tmp(nysg:nyng,nxlg:nxrg) )
603	!
604	!-- Initialize 2D-index arrays. Note, these will be removed soon!
605	nzb_local(nys:nyn,nxl:nxr) = get_topography_top_index( 's' )
606	CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
607
608	IF ( TRIM( topography ) /= 'flat' ) THEN
609	#if defined( __parallel )
610	CALL MPI_ALLREDUCE( MAXVAL( get_topography_top_index( 's' ) ), &
611	nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
612	CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ), &
613	nzb_local_min, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
614	#else
615	nzb_local_max = MAXVAL( get_topography_top_index( 's' ) )
616	nzb_local_min = MINVAL( get_topography_top_index( 's' ) )
617	#endif
618	!
619	!-- Consistency checks
620	IF ( nzb_local_min < 0 .OR. nzb_local_max > nz + 1 ) THEN
621	WRITE( message_string, * ) 'nzb_local values are outside the', &
622	'model domain', &
623	'&MINVAL( nzb_local ) = ', nzb_local_min, &
624	'&MAXVAL( nzb_local ) = ', nzb_local_max
625	CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
626	ENDIF
627	ENDIF
628
629	nzb_s_inner = nzb; nzb_s_outer = nzb
630	nzb_u_inner = nzb; nzb_u_outer = nzb
631	nzb_v_inner = nzb; nzb_v_outer = nzb
632	nzb_w_inner = nzb; nzb_w_outer = nzb
633
634	!
635	!-- Define vertical gridpoint from (or to) which on the usual finite difference
636	!-- form (which does not use surface fluxes) is applied
637	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
638	nzb_diff = nzb + 2
639	ELSE
640	nzb_diff = nzb + 1
641	ENDIF
642
643	nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff
644	!
645	!-- Set Neumann conditions for topography. Will be removed soon.
646	IF ( .NOT. bc_ns_cyc ) THEN
647	IF ( nys == 0 ) THEN
648	nzb_local(-1,:) = nzb_local(0,:)
649	ELSEIF ( nyn == ny ) THEN
650	nzb_local(ny+1,:) = nzb_local(ny,:)
651	ENDIF
652	ENDIF
653
654	IF ( .NOT. bc_lr_cyc ) THEN
655	IF ( nxl == 0 ) THEN
656	nzb_local(:,-1) = nzb_local(:,0)
657	nzb_local(:,-2) = nzb_local(:,0)
658	nzb_local(:,-3) = nzb_local(:,0)
659	ELSEIF ( nxr == nx ) THEN
660	nzb_local(:,nx+1) = nzb_local(:,nx)
661	nzb_local(:,nx+2) = nzb_local(:,nx)
662	nzb_local(:,nx+3) = nzb_local(:,nx)
663	ENDIF
664	ENDIF
665	!
666	!-- Initialization of 2D index arrays, will be removed soon!
667	!-- Initialize nzb_s_inner and nzb_w_inner
668	nzb_s_inner = nzb_local
669	nzb_w_inner = nzb_local
670
671	!
672	!-- Initialize remaining index arrays:
673	!-- first pre-initialize them with nzb_s_inner...
674	nzb_u_inner = nzb_s_inner
675	nzb_u_outer = nzb_s_inner
676	nzb_v_inner = nzb_s_inner
677	nzb_v_outer = nzb_s_inner
678	nzb_w_outer = nzb_s_inner
679	nzb_s_outer = nzb_s_inner
680
681	!
682	!-- nzb_s_outer:
683	!-- extend nzb_local east-/westwards first, then north-/southwards
684	nzb_tmp = nzb_local
685	DO j = nys, nyn
686	DO i = nxl, nxr
687	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), &
688	nzb_local(j,i+1) )
689	ENDDO
690	ENDDO
691
692	CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
693
694	DO i = nxl, nxr
695	DO j = nys, nyn
696	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
697	nzb_tmp(j+1,i) )
698	ENDDO
699	!
700	!-- non-cyclic boundary conditions (overwritten by call of
701	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
702	IF ( nys == 0 ) THEN
703	j = -1
704	nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
705	ENDIF
706	IF ( nyn == ny ) THEN
707	j = ny + 1
708	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
709	ENDIF
710	ENDDO
711	!
712	!-- nzb_w_outer:
713	!-- identical to nzb_s_outer
714	nzb_w_outer = nzb_s_outer
715	!
716	!-- nzb_u_inner:
717	!-- extend nzb_local rightwards only
718	nzb_tmp = nzb_local
719	DO j = nys, nyn
720	DO i = nxl, nxr
721	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
722	ENDDO
723	ENDDO
724
725	CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
726
727	nzb_u_inner = nzb_tmp
728	!
729	!-- nzb_u_outer:
730	!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
731	DO i = nxl, nxr
732	DO j = nys, nyn
733	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
734	nzb_tmp(j+1,i) )
735	ENDDO
736	!
737	!-- non-cyclic boundary conditions (overwritten by call of
738	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
739	IF ( nys == 0 ) THEN
740	j = -1
741	nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
742	ENDIF
743	IF ( nyn == ny ) THEN
744	j = ny + 1
745	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
746	ENDIF
747	ENDDO
748
749	!
750	!-- nzb_v_inner:
751	!-- extend nzb_local northwards only
752	nzb_tmp = nzb_local
753	DO i = nxl, nxr
754	DO j = nys, nyn
755	nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
756	ENDDO
757	ENDDO
758
759	CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
760	nzb_v_inner = nzb_tmp
761
762	!
763	!-- nzb_v_outer:
764	!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
765	DO j = nys, nyn
766	DO i = nxl, nxr
767	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &
768	nzb_tmp(j,i+1) )
769	ENDDO
770	!
771	!-- non-cyclic boundary conditions (overwritten by call of
772	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
773	IF ( nxl == 0 ) THEN
774	i = -1
775	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
776	ENDIF
777	IF ( nxr == nx ) THEN
778	i = nx + 1
779	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
780	ENDIF
781	ENDDO
782
783	!
784	!-- Exchange of lateral boundary values (parallel computers) and cyclic
785	!-- boundary conditions, if applicable.
786	!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
787	!-- they do not require exchange and are not included here.
788	CALL exchange_horiz_2d_int( nzb_u_inner, nys, nyn, nxl, nxr, nbgp )
789	CALL exchange_horiz_2d_int( nzb_u_outer, nys, nyn, nxl, nxr, nbgp )
790	CALL exchange_horiz_2d_int( nzb_v_inner, nys, nyn, nxl, nxr, nbgp )
791	CALL exchange_horiz_2d_int( nzb_v_outer, nys, nyn, nxl, nxr, nbgp )
792	CALL exchange_horiz_2d_int( nzb_w_outer, nys, nyn, nxl, nxr, nbgp )
793	CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp )
794
795	!
796	!-- Set the individual index arrays which define the k index from which on
797	!-- the usual finite difference form (which does not use surface fluxes) is
798	!-- applied
799	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
800	nzb_diff_s_inner = nzb_s_inner + 2
801	nzb_diff_s_outer = nzb_s_outer + 2
802	ELSE
803	nzb_diff_s_inner = nzb_s_inner + 1
804	nzb_diff_s_outer = nzb_s_outer + 1
805	ENDIF
806	!
807	!-- Vertical nesting: communicate vertical grid level arrays between fine and
808	!-- coarse grid
809	IF ( vnested ) CALL vnest_init_grid
810
811	END SUBROUTINE init_grid
812
813	! Description:
814	! -----------------------------------------------------------------------------!
815	!> Set temporary topography flags and reference buildings on top of underlying
816	!> orography.
817	!------------------------------------------------------------------------------!
818	SUBROUTINE process_topography( topo_3d )
819
820	USE arrays_3d, &
821	ONLY: zu, zw
822
823	USE control_parameters, &
824	ONLY: bc_lr_cyc, bc_ns_cyc, land_surface, ocean, urban_surface
825
826	USE indices, &
827	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, &
828	nzt
829
830	USE netcdf_data_input_mod, &
831	ONLY: buildings_f, building_id_f, input_pids_static, &
832	terrain_height_f
833
834	USE kinds
835
836	USE pegrid
837
838	IMPLICIT NONE
839
840	INTEGER(iwp) :: i !< running index along x-direction
841	INTEGER(iwp) :: j !< running index along y-direction
842	INTEGER(iwp) :: k !< running index along z-direction with respect to numeric grid
843	INTEGER(iwp) :: k2 !< running index along z-direction with respect to netcdf grid
844	INTEGER(iwp) :: nr !< index variable indication maximum terrain height for respective building ID
845	INTEGER(iwp) :: num_build !< counter for number of buildings
846	INTEGER(iwp) :: topo_top_index !< orography top index, used to map 3D buildings onto terrain
847
848	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displace_dum !< displacements of start addresses, used for MPI_ALLGATHERV
849	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids !< building IDs on entire model domain
850	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final !< building IDs on entire model domain, multiple occurences are sorted out
851	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final_tmp !< temporary array used for resizing
852	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l !< building IDs on local subdomain
853	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l_tmp !< temporary array used to resize array of building IDs
854
855	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings !< number of buildings with different ID on entire model domain
856	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings_l !< number of buildings with different ID on local subdomain
857
858	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags
859
860	REAL(wp) :: ocean_offset !< offset to consider inverse vertical coordinate at topography definition
861	REAL(wp), DIMENSION(:), ALLOCATABLE :: oro_max !< maximum terrain height occupied by an building with certain id
862	REAL(wp), DIMENSION(:), ALLOCATABLE :: oro_max_l !< maximum terrain height occupied by an building with certain id, on local subdomain
863
864	!
865	!-- In the following, buildings and orography are further preprocessed
866	!-- before they are mapped on the LES grid.
867	!-- Buildings are mapped on top of the orography by maintaining the roof
868	!-- shape of the building. This can be achieved by referencing building on
869	!-- top of the maximum terrain height within the area occupied by the
870	!-- respective building. As buildings and terrain height are defined PE-wise,
871	!-- parallelization of this referencing is required (a building can be
872	!-- distributed between different PEs).
873	!-- In a first step, determine the number of buildings with different
874	!-- building id on each PE. In a next step, all building ids are gathered
875	!-- into one array which is present to all PEs. For each building ID,
876	!-- the maximum terrain height occupied by the respective building is
877	!-- computed and distributed to each PE.
878	!-- Finally, for each building id and its respective reference orography,
879	!-- builidings are mapped on top.
880	!--
881	!-- First, pre-set topography flags, bit 1 indicates orography, bit 2
882	!-- buildings
883	!-- classify the respective surfaces.
884	topo_3d = IBSET( topo_3d, 0 )
885	topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
886	!
887	!-- Reference buildings on top of orography. This is not necessary
888	!-- if topography is read from ASCII file as no distinction between buildings
889	!-- and terrain height can be made. Moreover, this is also not necessary if
890	!-- urban-surface and land-surface model are used at the same time.
891	IF ( input_pids_static ) THEN
892	num_buildings_l = 0
893	num_buildings = 0
894	!
895	!-- Allocate at least one element for building ids,
896	ALLOCATE( build_ids_l(1) )
897	DO i = nxl, nxr
898	DO j = nys, nyn
899	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
900	IF ( num_buildings_l(myid) > 0 ) THEN
901	IF ( ANY( building_id_f%var(j,i) .EQ. build_ids_l ) ) THEN
902	CYCLE
903	ELSE
904	num_buildings_l(myid) = num_buildings_l(myid) + 1
905	!
906	!-- Resize array with different local building ids
907	ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
908	build_ids_l_tmp = build_ids_l
909	DEALLOCATE( build_ids_l )
910	ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
911	build_ids_l(1:num_buildings_l(myid)-1) = &
912	build_ids_l_tmp(1:num_buildings_l(myid)-1)
913	build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
914	DEALLOCATE( build_ids_l_tmp )
915	ENDIF
916	!
917	!-- First occuring building id on PE
918	ELSE
919	num_buildings_l(myid) = num_buildings_l(myid) + 1
920	build_ids_l(1) = building_id_f%var(j,i)
921	ENDIF
922	ENDIF
923	ENDDO
924	ENDDO
925	!
926	!-- Determine number of different building ids for the entire domain
927	#if defined( __parallel )
928	CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs, &
929	MPI_INTEGER, MPI_SUM, comm2d, ierr )
930	#else
931	num_buildings = num_buildings_l
932	#endif
933	!
934	!-- Gather all buildings ids on each PEs.
935	!-- First, allocate array encompassing all building ids in model domain.
936	ALLOCATE( build_ids(1:SUM(num_buildings)) )
937	#if defined( __parallel )
938	!
939	!-- Allocate array for displacements.
940	!-- As each PE may has a different number of buildings, so that
941	!-- the block sizes send by each PE may not be equal. Hence,
942	!-- information about the respective displacement is required, indicating
943	!-- the respective adress where each MPI-task writes into the receive
944	!-- buffer array
945	ALLOCATE( displace_dum(0:numprocs-1) )
946	displace_dum(0) = 0
947	DO i = 1, numprocs-1
948	displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
949	ENDDO
950
951	CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)), &
952	num_buildings(myid), &
953	MPI_INTEGER, &
954	build_ids, &
955	num_buildings, &
956	displace_dum, &
957	MPI_INTEGER, &
958	comm2d, ierr )
959
960	DEALLOCATE( displace_dum )
961
962	#else
963	build_ids = build_ids_l
964	#endif
965
966	!
967	!-- Note, in parallel mode building ids can occure mutliple times, as
968	!-- each PE has send its own ids. Therefore, sort out building ids which
969	!-- appear more than one time.
970	num_build = 0
971	DO nr = 1, SIZE(build_ids)
972
973	IF ( ALLOCATED(build_ids_final) ) THEN
974	IF ( ANY( build_ids(nr) .EQ. build_ids_final ) ) THEN
975	CYCLE
976	ELSE
977	num_build = num_build + 1
978	!
979	!-- Resize
980	ALLOCATE( build_ids_final_tmp(1:num_build) )
981	build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
982	DEALLOCATE( build_ids_final )
983	ALLOCATE( build_ids_final(1:num_build) )
984	build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
985	build_ids_final(num_build) = build_ids(nr)
986	DEALLOCATE( build_ids_final_tmp )
987	ENDIF
988	ELSE
989	num_build = num_build + 1
990	ALLOCATE( build_ids_final(1:num_build) )
991	build_ids_final(num_build) = build_ids(nr)
992	ENDIF
993	ENDDO
994
995	!
996	!-- Finally, determine maximumum terrain height occupied by the
997	!-- respective building.
998	ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
999	ALLOCATE( oro_max(1:SIZE(build_ids_final)) )
1000	oro_max_l = 0.0_wp
1001
1002	DO nr = 1, SIZE(build_ids_final)
1003	oro_max_l(nr) = MAXVAL( &
1004	MERGE( terrain_height_f%var, 0.0_wp, &
1005	building_id_f%var(nys:nyn,nxl:nxr) .EQ. &
1006	build_ids_final(nr) ) )
1007	ENDDO
1008
1009	#if defined( __parallel )
1010	IF ( SIZE(build_ids_final) >= 1 ) THEN
1011	CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL, &
1012	MPI_MAX, comm2d, ierr )
1013	ENDIF
1014	#else
1015	oro_max = oro_max_l
1016	#endif
1017	!
1018	!-- Map orography as well as buildings onto grid.
1019	!-- In case of ocean simulations, add an offset.
1020	ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
1021	DO i = nxl, nxr
1022	DO j = nys, nyn
1023	topo_top_index = 0
1024	DO k = nzb, nzt
1025	!
1026	!-- In a first step, if grid point is below or equal the given
1027	!-- terrain height, grid point is flagged to be of type natural.
1028	!-- Please note, in case there is also a building which is lower
1029	!-- than the vertical grid spacing, initialization of surface
1030	!-- attributes will not be correct as given surface information
1031	!-- will not be in accordance to the classified grid points.
1032	!-- Hence, in this case, de-flag the grid point and give it
1033	!-- urban type instead.
1034	IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) ) THEN
1035	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1036	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1037	topo_top_index = topo_top_index + 1
1038	ENDIF
1039	!
1040	!-- Set building grid points. Here, only consider 2D buildings.
1041	!-- 3D buildings require separate treatment.
1042	IF ( buildings_f%lod == 1 ) THEN
1043	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1044	!
1045	!-- Determine index where maximum terrain height occupied by
1046	!-- the respective building height is stored.
1047	nr = MINLOC( ABS( build_ids_final - &
1048	building_id_f%var(j,i) ), DIM = 1 )
1049
1050	IF ( zu(k) - ocean_offset <= &
1051	oro_max(nr) + buildings_f%var_2d(j,i) ) THEN
1052	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1053	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1054	!
1055	!-- De-flag grid point of type natural. See comment above.
1056	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
1057	ENDIF
1058	ENDIF
1059	ENDIF
1060	ENDDO
1061	!
1062	!-- Map 3D buildings onto terrain height.
1063	!-- In case of any slopes, map building on top of maximum terrain
1064	!-- height covered by the building. In other words, extend
1065	!-- building down to the respective local terrain-surface height.
1066	IF ( buildings_f%lod == 2 ) THEN
1067	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1068	!
1069	!-- Determine index for maximum-terrain-height array.
1070	nr = MINLOC( ABS( build_ids_final - &
1071	building_id_f%var(j,i) ), DIM = 1 )
1072	!
1073	!-- Extend building down to the terrain surface.
1074	k2 = topo_top_index
1075	DO k = topo_top_index + 1, nzt + 1
1076	IF ( zu(k) - ocean_offset <= oro_max(nr) ) THEN
1077	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1078	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1079	k2 = k2 + 1
1080	ENDIF
1081	ENDDO
1082	topo_top_index = k2
1083	!
1084	!-- Now, map building on top.
1085	k2 = 0
1086	DO k = topo_top_index, nzt + 1
1087	IF ( k2 <= buildings_f%nz-1 ) THEN
1088	IF ( buildings_f%var_3d(k2,j,i) == 1 ) THEN
1089	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1090	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
1091	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1092	ENDIF
1093	ENDIF
1094	k2 = k2 + 1
1095	ENDDO
1096	ENDIF
1097	ENDIF
1098	ENDDO
1099	ENDDO
1100	!
1101	!-- Deallocate temporary arrays required for processing and reading data
1102	IF ( ALLOCATED( oro_max ) ) DEALLOCATE( oro_max )
1103	IF ( ALLOCATED( oro_max_l ) ) DEALLOCATE( oro_max_l )
1104	IF ( ALLOCATED( build_ids_final ) ) DEALLOCATE( build_ids_final )
1105	!
1106	!-- Topography input via ASCII format.
1107	ELSE
1108	ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
1109	topo_3d = IBSET( topo_3d, 0 )
1110	topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
1111	DO i = nxl, nxr
1112	DO j = nys, nyn
1113	DO k = nzb, nzt
1114	IF ( zu(k) - ocean_offset <= buildings_f%var_2d(j,i) ) THEN
1115	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1116	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates terrain
1117	ENDIF
1118	ENDDO
1119	ENDDO
1120	ENDDO
1121	ENDIF
1122
1123	CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1124
1125	IF ( .NOT. bc_ns_cyc ) THEN
1126	IF ( nys == 0 ) topo_3d(:,-1,:) = topo_3d(:,0,:)
1127	IF ( nyn == ny ) topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
1128	ENDIF
1129
1130	IF ( .NOT. bc_lr_cyc ) THEN
1131	IF ( nxl == 0 ) topo_3d(:,:,-1) = topo_3d(:,:,0)
1132	IF ( nxr == nx ) topo_3d(:,:,nx+1) = topo_3d(:,:,nx)
1133	ENDIF
1134
1135	END SUBROUTINE process_topography
1136
1137
1138	! Description:
1139	! -----------------------------------------------------------------------------!
1140	!> Filter topography, i.e. fill holes resolved by only one grid point.
1141	!> Such holes are suspected to lead to velocity blow-ups as continuity
1142	!> equation on discrete grid cannot be fulfilled in such case.
1143	!------------------------------------------------------------------------------!
1144	SUBROUTINE filter_topography( topo_3d )
1145
1146	USE control_parameters, &
1147	ONLY: bc_lr_cyc, bc_ns_cyc, message_string
1148
1149	USE indices, &
1150	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
1151
1152	USE netcdf_data_input_mod, &
1153	ONLY: building_id_f, building_type_f
1154
1155	USE pegrid
1156
1157	IMPLICIT NONE
1158
1159	INTEGER(iwp) :: i !< running index along x-direction
1160	INTEGER(iwp) :: it !< counter for number of iterations
1161	INTEGER(iwp) :: j !< running index along y-direction
1162	INTEGER(iwp) :: k !< running index along z-direction
1163	INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point
1164	INTEGER(iwp) :: num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE
1165	INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point
1166
1167	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo_tmp !< temporary 3D-topography used to fill holes
1168	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo_3d !< 3D-topography array merging buildings and orography
1169	!
1170	!-- Before checking for holes, set lateral boundary conditions for
1171	!-- topography. After hole-filling, boundary conditions must be set again.
1172	!-- Several iterations are performed, in order to fill holes which might
1173	!-- emerge by the filling-algorithm itself.
1174	ALLOCATE( topo_tmp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1175	topo_tmp = 0
1176
1177	num_hole = 99999
1178	it = 0
1179	DO WHILE ( num_hole > 0 )
1180
1181	num_hole = 0
1182	CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1183
1184	topo_tmp = topo_3d
1185	!
1186	!-- In case of non-cyclic lateral boundaries, assume lateral boundary to be
1187	!-- a solid wall. Thus, intermediate spaces of one grid point between
1188	!-- boundary and some topographic structure will be filled.
1189	IF ( .NOT. bc_ns_cyc ) THEN
1190	IF ( nys == 0 ) topo_tmp(:,-1,:) = IBCLR( topo_tmp(:,0,:), 0 )
1191	IF ( nyn == ny ) topo_tmp(:,ny+1,:) = IBCLR( topo_tmp(:,ny,:), 0 )
1192	ENDIF
1193
1194	IF ( .NOT. bc_lr_cyc ) THEN
1195	IF ( nxl == 0 ) topo_tmp(:,:,-1) = IBCLR( topo_tmp(:,:,0), 0 )
1196	IF ( nxr == nx ) topo_tmp(:,:,nx+1) = IBCLR( topo_tmp(:,:,nx), 0 )
1197	ENDIF
1198
1199	it = it + 1
1200	num_hole_l = 0
1201	DO i = nxl, nxr
1202	DO j = nys, nyn
1203	DO k = nzb+1, nzt
1204	IF ( BTEST( topo_tmp(k,j,i), 0 ) ) THEN
1205	num_wall = 0
1206	IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) ) &
1207	num_wall = num_wall + 1
1208	IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) ) &
1209	num_wall = num_wall + 1
1210	IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) ) &
1211	num_wall = num_wall + 1
1212	IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) ) &
1213	num_wall = num_wall + 1
1214	IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) ) &
1215	num_wall = num_wall + 1
1216	IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) ) &
1217	num_wall = num_wall + 1
1218
1219	IF ( num_wall >= 4 ) THEN
1220	num_hole_l = num_hole_l + 1
1221	!
1222	!-- Clear flag 0 and set special flag ( bit 3) to indicate
1223	!-- that new topography point is a result of filtering process.
1224	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1225	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 3 )
1226	!
1227	!-- If filled grid point is occupied by a building, classify
1228	!-- it as building grid point.
1229	IF ( building_type_f%from_file ) THEN
1230	IF ( building_type_f%var(j,i) /= &
1231	building_type_f%fill .OR. &
1232	building_type_f%var(j+1,i) /= &
1233	building_type_f%fill .OR. &
1234	building_type_f%var(j-1,i) /= &
1235	building_type_f%fill .OR. &
1236	building_type_f%var(j,i+1) /= &
1237	building_type_f%fill .OR. &
1238	building_type_f%var(j,i-1) /= &
1239	building_type_f%fill ) THEN
1240	!
1241	!-- Set flag indicating building surfaces
1242	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1243	!
1244	!-- Set building_type and ID at this position if not
1245	!-- already set. This is required for proper
1246	!-- initialization of urban-surface energy balance
1247	!-- solver.
1248	IF ( building_type_f%var(j,i) == &
1249	building_type_f%fill ) THEN
1250
1251	IF ( building_type_f%var(j+1,i) /= &
1252	building_type_f%fill ) THEN
1253	building_type_f%var(j,i) = &
1254	building_type_f%var(j+1,i)
1255	building_id_f%var(j,i) = &
1256	building_id_f%var(j+1,i)
1257	ELSEIF ( building_type_f%var(j-1,i) /= &
1258	building_type_f%fill ) THEN
1259	building_type_f%var(j,i) = &
1260	building_type_f%var(j-1,i)
1261	building_id_f%var(j,i) = &
1262	building_id_f%var(j-1,i)
1263	ELSEIF ( building_type_f%var(j,i+1) /= &
1264	building_type_f%fill ) THEN
1265	building_type_f%var(j,i) = &
1266	building_type_f%var(j,i+1)
1267	building_id_f%var(j,i) = &
1268	building_id_f%var(j,i+1)
1269	ELSEIF ( building_type_f%var(j,i-1) /= &
1270	building_type_f%fill ) THEN
1271	building_type_f%var(j,i) = &
1272	building_type_f%var(j,i-1)
1273	building_id_f%var(j,i) = &
1274	building_id_f%var(j,i-1)
1275	ENDIF
1276	ENDIF
1277	ENDIF
1278	ENDIF
1279	!
1280	!-- If filled grid point is already classified as building
1281	!-- everything is fine, else classify this grid point as
1282	!-- natural type grid point. This case, values for the
1283	!-- surface type are already set.
1284	IF ( .NOT. BTEST( topo_3d(k,j,i), 2 ) ) THEN
1285	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1286	ENDIF
1287	ENDIF
1288	ENDIF
1289	ENDDO
1290	ENDDO
1291	ENDDO
1292	!
1293	!-- Count the total number of holes, required for informative message.
1294	#if defined( __parallel )
1295	CALL MPI_ALLREDUCE( num_hole_l, num_hole, 1, MPI_INTEGER, MPI_SUM, &
1296	comm2d, ierr )
1297	#else
1298	num_hole = num_hole_l
1299	#endif
1300
1301	!
1302	!-- Create an informative message if any holes were filled.
1303	IF ( num_hole > 0 ) THEN
1304	WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '// &
1305	'one grid point were filled at ', it, &
1306	'th iteration '
1307	CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 )
1308	ENDIF
1309
1310
1311	ENDDO
1312
1313	DEALLOCATE( topo_tmp )
1314	!
1315	!-- Finally, exchange topo_3d array again and if necessary set Neumann boundary
1316	!-- condition in case of non-cyclic lateral boundaries.
1317	CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1318
1319	IF ( .NOT. bc_ns_cyc ) THEN
1320	IF ( nys == 0 ) topo_3d(:,-1,:) = topo_3d(:,0,:)
1321	IF ( nyn == ny ) topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
1322	ENDIF
1323
1324	IF ( .NOT. bc_lr_cyc ) THEN
1325	IF ( nxl == 0 ) topo_3d(:,:,-1) = topo_3d(:,:,0)
1326	IF ( nxr == nx ) topo_3d(:,:,nx+1) = topo_3d(:,:,nx)
1327	ENDIF
1328
1329	END SUBROUTINE filter_topography
1330
1331	! Description:
1332	! -----------------------------------------------------------------------------!
1333	!> Pre-computation of grid-dependent and near-wall mixing length.
1334	!> @todo: move subroutine to turbulence module developed by M1 (Tobi).
1335	!------------------------------------------------------------------------------!
1336	SUBROUTINE init_mixing_length
1337
1338	USE arrays_3d, &
1339	ONLY: dzw, l_grid, l_wall, zu, zw
1340
1341	USE control_parameters, &
1342	ONLY: message_string, wall_adjustment_factor
1343
1344	USE grid_variables, &
1345	ONLY: dx, dy
1346
1347	USE indices, &
1348	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
1349	wall_flags_0
1350
1351	USE kinds
1352
1353	IMPLICIT NONE
1354
1355	INTEGER(iwp) :: i !< index variable along x
1356	INTEGER(iwp) :: j !< index variable along y
1357	INTEGER(iwp) :: k !< index variable along z
1358
1359	ALLOCATE( l_grid(1:nzt) )
1360	ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1361	!
1362	!-- Compute the grid-dependent mixing length.
1363	DO k = 1, nzt
1364	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp
1365	ENDDO
1366	!
1367	!-- Initialize near-wall mixing length l_wall only in the vertical direction
1368	!-- for the moment, multiplication with wall_adjustment_factor further below
1369	l_wall(nzb,:,:) = l_grid(1)
1370	DO k = nzb+1, nzt
1371	l_wall(k,:,:) = l_grid(k)
1372	ENDDO
1373	l_wall(nzt+1,:,:) = l_grid(nzt)
1374
1375	DO k = 1, nzt
1376	IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. &
1377	l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN
1378	WRITE( message_string, * ) 'grid anisotropy exceeds ', &
1379	'threshold given by only local', &
1380	' &horizontal reduction of near_wall ', &
1381	'mixing length l_wall', &
1382	' &starting from height level k = ', k, '.'
1383	CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
1384	EXIT
1385	ENDIF
1386	ENDDO
1387	!
1388	!-- In case of topography: limit near-wall mixing length l_wall further:
1389	!-- Go through all points of the subdomain one by one and look for the closest
1390	!-- surface.
1391	!-- Is this correct in the ocean case?
1392	DO i = nxl, nxr
1393	DO j = nys, nyn
1394	DO k = nzb+1, nzt
1395	!
1396	!-- Check if current gridpoint belongs to the atmosphere
1397	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1398	!
1399	!-- Check for neighbouring grid-points.
1400	!-- Vertical distance, down
1401	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) &
1402	l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
1403	!
1404	!-- Vertical distance, up
1405	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
1406	l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
1407	!
1408	!-- y-distance
1409	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) .OR. &
1410	.NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) &
1411	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
1412	!
1413	!-- x-distance
1414	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) .OR. &
1415	.NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) &
1416	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
1417	!
1418	!-- yz-distance (vertical edges, down)
1419	IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 ) .OR. &
1420	.NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 ) ) &
1421	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1422	SQRT( 0.25_wp * dy**2 + &
1423	( zu(k) - zw(k-1) )**2 ) )
1424	!
1425	!-- yz-distance (vertical edges, up)
1426	IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 ) .OR. &
1427	.NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 ) ) &
1428	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1429	SQRT( 0.25_wp * dy**2 + &
1430	( zw(k) - zu(k) )**2 ) )
1431	!
1432	!-- xz-distance (vertical edges, down)
1433	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 ) .OR. &
1434	.NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 ) ) &
1435	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1436	SQRT( 0.25_wp * dx**2 + &
1437	( zu(k) - zw(k-1) )**2 ) )
1438	!
1439	!-- xz-distance (vertical edges, up)
1440	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 ) .OR. &
1441	.NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 ) ) &
1442	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1443	SQRT( 0.25_wp * dx**2 + &
1444	( zw(k) - zu(k) )**2 ) )
1445	!
1446	!-- xy-distance (horizontal edges)
1447	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 ) .OR. &
1448	.NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 ) .OR. &
1449	.NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 ) .OR. &
1450	.NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) ) &
1451	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1452	SQRT( 0.25_wp * ( dx2 + dy2 ) ) )
1453	!
1454	!-- xyz distance (vertical and horizontal edges, down)
1455	IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 ) .OR. &
1456	.NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 ) .OR. &
1457	.NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 ) .OR. &
1458	.NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) ) &
1459	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1460	SQRT( 0.25_wp * ( dx2 + dy2 ) &
1461	+ ( zu(k) - zw(k-1) )**2 ) )
1462	!
1463	!-- xyz distance (vertical and horizontal edges, up)
1464	IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 ) .OR. &
1465	.NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 ) .OR. &
1466	.NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 ) .OR. &
1467	.NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) ) &
1468	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1469	SQRT( 0.25_wp * ( dx2 + dy2 ) &
1470	+ ( zw(k) - zu(k) )**2 ) )
1471
1472	ENDIF
1473	ENDDO
1474	ENDDO
1475	ENDDO
1476	!
1477	!-- Set lateral boundary conditions for l_wall
1478	CALL exchange_horiz( l_wall, nbgp )
1479
1480	END SUBROUTINE init_mixing_length
1481
1482	! Description:
1483	! -----------------------------------------------------------------------------!
1484	!> Reads topography information from file or sets generic topography. Moreover,
1485	!> all topography-relevant topography arrays are initialized, and grid flags
1486	!> are set.
1487	!------------------------------------------------------------------------------!
1488	SUBROUTINE init_topo( topo )
1489
1490	USE arrays_3d, &
1491	ONLY: zw
1492
1493	USE control_parameters, &
1494	ONLY: bc_lr_cyc, bc_ns_cyc, building_height, building_length_x, &
1495	building_length_y, building_wall_left, building_wall_south, &
1496	canyon_height, canyon_wall_left, canyon_wall_south, &
1497	canyon_width_x, canyon_width_y, dp_level_ind_b, dz, &
1498	message_string, ocean, topography, topography_grid_convention, &
1499	tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y, &
1500	tunnel_wall_depth
1501
1502	USE grid_variables, &
1503	ONLY: dx, dy
1504
1505	USE indices, &
1506	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, &
1507	nzb, nzt
1508
1509	USE kinds
1510
1511	USE pegrid
1512
1513	USE surface_mod, &
1514	ONLY: get_topography_top_index, get_topography_top_index_ji
1515
1516	IMPLICIT NONE
1517
1518	INTEGER(iwp) :: bh !< temporary vertical index of building height
1519	INTEGER(iwp) :: blx !< grid point number of building size along x
1520	INTEGER(iwp) :: bly !< grid point number of building size along y
1521	INTEGER(iwp) :: bxl !< index for left building wall
1522	INTEGER(iwp) :: bxr !< index for right building wall
1523	INTEGER(iwp) :: byn !< index for north building wall
1524	INTEGER(iwp) :: bys !< index for south building wall
1525	INTEGER(iwp) :: ch !< temporary vertical index for canyon height
1526	INTEGER(iwp) :: cwx !< grid point number of canyon size along x
1527	INTEGER(iwp) :: cwy !< grid point number of canyon size along y
1528	INTEGER(iwp) :: cxl !< index for left canyon wall
1529	INTEGER(iwp) :: cxr !< index for right canyon wall
1530	INTEGER(iwp) :: cyn !< index for north canyon wall
1531	INTEGER(iwp) :: cys !< index for south canyon wall
1532	INTEGER(iwp) :: i !< index variable along x
1533	INTEGER(iwp) :: it !< Number of iterations required to fill all problematic holes
1534	INTEGER(iwp) :: j !< index variable along y
1535	INTEGER(iwp) :: k !< index variable along z
1536	INTEGER(iwp) :: hv_in !< heavyside function to model inner tunnel surface
1537	INTEGER(iwp) :: hv_out !< heavyside function to model outer tunnel surface
1538	INTEGER(iwp) :: txe_out !< end position of outer tunnel wall in x
1539	INTEGER(iwp) :: txs_out !< start position of outer tunnel wall in x
1540	INTEGER(iwp) :: tye_out !< end position of outer tunnel wall in y
1541	INTEGER(iwp) :: tys_out !< start position of outer tunnel wall in y
1542	INTEGER(iwp) :: txe_in !< end position of inner tunnel wall in x
1543	INTEGER(iwp) :: txs_in !< start position of inner tunnel wall in x
1544	INTEGER(iwp) :: tye_in !< end position of inner tunnel wall in y
1545	INTEGER(iwp) :: tys_in !< start position of inner tunnel wall in y
1546	INTEGER(iwp) :: td !< tunnel wall depth
1547	INTEGER(iwp) :: th !< height of outer tunnel wall
1548
1549	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center
1550	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo !< input array for 3D topography and dummy array for setting "outer"-flags
1551
1552
1553	!
1554	!-- Set outer and inner index arrays for non-flat topography.
1555	!-- Here consistency checks concerning domain size and periodicity are
1556	!-- necessary.
1557	!-- Within this SELECT CASE structure only nzb_local is initialized
1558	!-- individually depending on the chosen topography type, all other index
1559	!-- arrays are initialized further below.
1560	SELECT CASE ( TRIM( topography ) )
1561
1562	CASE ( 'flat' )
1563	!
1564	!-- Initialilize 3D topography array, used later for initializing flags
1565	topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
1566
1567	CASE ( 'single_building' )
1568	!
1569	!-- Single rectangular building, by default centered in the middle of the
1570	!-- total domain
1571	blx = NINT( building_length_x / dx )
1572	bly = NINT( building_length_y / dy )
1573	bh = MINLOC( ABS( zw - building_height ), 1 ) - 1
1574	IF ( ABS( zw(bh) - building_height ) == &
1575	ABS( zw(bh+1) - building_height ) ) bh = bh + 1
1576	IF ( building_wall_left == 9999999.9_wp ) THEN
1577	building_wall_left = ( nx + 1 - blx ) / 2 * dx
1578	ENDIF
1579	bxl = NINT( building_wall_left / dx )
1580	bxr = bxl + blx
1581
1582	IF ( building_wall_south == 9999999.9_wp ) THEN
1583	building_wall_south = ( ny + 1 - bly ) / 2 * dy
1584	ENDIF
1585	bys = NINT( building_wall_south / dy )
1586	byn = bys + bly
1587
1588	!
1589	!-- Building size has to meet some requirements
1590	IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR. &
1591	( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) ) THEN
1592	WRITE( message_string, * ) 'inconsistent building parameters:', &
1593	'& bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
1594	'byn=', byn, 'nx=', nx, 'ny=', ny
1595	CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
1596	ENDIF
1597
1598	ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
1599	!
1600	!-- Define the building.
1601	IF ( bxl <= nxr .AND. bxr >= nxl .AND. &
1602	bys <= nyn .AND. byn >= nys ) &
1603	nzb_local(MAX(nys,bys):MIN(nyn,byn),MAX(nxl,bxl):MIN(nxr,bxr)) = bh
1604	!
1605	!-- Set bit array on basis of nzb_local
1606	DO i = nxl, nxr
1607	DO j = nys, nyn
1608	topo(nzb_local(j,i)+1:nzt+1,j,i) = &
1609	IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 )
1610	ENDDO
1611	ENDDO
1612
1613	DEALLOCATE( nzb_local )
1614
1615	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1616	!
1617	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
1618	!-- boundary conditions for topography.
1619	IF ( .NOT. bc_ns_cyc ) THEN
1620	IF ( nys == 0 ) THEN
1621	DO i = 1, nbgp
1622	topo(:,nys-i,:) = topo(:,nys,:)
1623	ENDDO
1624	ENDIF
1625	IF ( nyn == ny ) THEN
1626	DO i = 1, nbgp
1627	topo(:,nyn+i,:) = topo(:,nyn,:)
1628	ENDDO
1629	ENDIF
1630	ENDIF
1631	IF ( .NOT. bc_lr_cyc ) THEN
1632	IF ( nxl == 0 ) THEN
1633	DO i = 1, nbgp
1634	topo(:,:,nxl-i) = topo(:,:,nxl)
1635	ENDDO
1636	ENDIF
1637	IF ( nxr == nx ) THEN
1638	DO i = 1, nbgp
1639	topo(:,:,nxr+i) = topo(:,:,nxr)
1640	ENDDO
1641	ENDIF
1642	ENDIF
1643
1644	CASE ( 'single_street_canyon' )
1645	!
1646	!-- Single quasi-2D street canyon of infinite length in x or y direction.
1647	!-- The canyon is centered in the other direction by default.
1648	IF ( canyon_width_x /= 9999999.9_wp ) THEN
1649	!
1650	!-- Street canyon in y direction
1651	cwx = NINT( canyon_width_x / dx )
1652	IF ( canyon_wall_left == 9999999.9_wp ) THEN
1653	canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
1654	ENDIF
1655	cxl = NINT( canyon_wall_left / dx )
1656	cxr = cxl + cwx
1657	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
1658	!
1659	!-- Street canyon in x direction
1660	cwy = NINT( canyon_width_y / dy )
1661	IF ( canyon_wall_south == 9999999.9_wp ) THEN
1662	canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
1663	ENDIF
1664	cys = NINT( canyon_wall_south / dy )
1665	cyn = cys + cwy
1666
1667	ELSE
1668
1669	message_string = 'no street canyon width given'
1670	CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
1671
1672	ENDIF
1673
1674	ch = MINLOC( ABS( zw - canyon_height ), 1 ) - 1
1675	IF ( ABS( zw(ch) - canyon_height ) == &
1676	ABS( zw(ch+1) - canyon_height ) ) ch = ch + 1
1677	dp_level_ind_b = ch
1678	!
1679	!-- Street canyon size has to meet some requirements
1680	IF ( canyon_width_x /= 9999999.9_wp ) THEN
1681	IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR. &
1682	( ch < 3 ) ) THEN
1683	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
1684	'&cxl=', cxl, 'cxr=', cxr, &
1685	'cwx=', cwx, &
1686	'ch=', ch, 'nx=', nx, 'ny=', ny
1687	CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
1688	ENDIF
1689	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
1690	IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR. &
1691	( ch < 3 ) ) THEN
1692	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
1693	'&cys=', cys, 'cyn=', cyn, &
1694	'cwy=', cwy, &
1695	'ch=', ch, 'nx=', nx, 'ny=', ny
1696	CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
1697	ENDIF
1698	ENDIF
1699	IF ( canyon_width_x /= 9999999.9_wp .AND. &
1700	canyon_width_y /= 9999999.9_wp ) THEN
1701	message_string = 'inconsistent canyon parameters:' // &
1702	'&street canyon can only be oriented' // &
1703	'&either in x- or in y-direction'
1704	CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
1705	ENDIF
1706
1707	ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
1708	nzb_local = ch
1709	IF ( canyon_width_x /= 9999999.9_wp ) THEN
1710	IF ( cxl <= nxr .AND. cxr >= nxl ) &
1711	nzb_local(:,MAX(nxl,cxl+1):MIN(nxr,cxr-1)) = 0
1712	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
1713	IF ( cys <= nyn .AND. cyn >= nys ) &
1714	nzb_local(MAX(nys,cys+1):MIN(nyn,cyn-1),:) = 0
1715	ENDIF
1716	!
1717	!-- Set bit array on basis of nzb_local
1718	DO i = nxl, nxr
1719	DO j = nys, nyn
1720	topo(nzb_local(j,i)+1:nzt+1,j,i) = &
1721	IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 )
1722	ENDDO
1723	ENDDO
1724	DEALLOCATE( nzb_local )
1725
1726	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1727	!
1728	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
1729	!-- boundary conditions for topography.
1730	IF ( .NOT. bc_ns_cyc ) THEN
1731	IF ( nys == 0 ) THEN
1732	DO i = 1, nbgp
1733	topo(:,nys-i,:) = topo(:,nys,:)
1734	ENDDO
1735	ENDIF
1736	IF ( nyn == ny ) THEN
1737	DO i = 1, nbgp
1738	topo(:,nyn+i,:) = topo(:,nyn,:)
1739	ENDDO
1740	ENDIF
1741	ENDIF
1742	IF ( .NOT. bc_lr_cyc ) THEN
1743	IF ( nxl == 0 ) THEN
1744	DO i = 1, nbgp
1745	topo(:,:,nxl-i) = topo(:,:,nxl)
1746	ENDDO
1747	ENDIF
1748	IF ( nxr == nx ) THEN
1749	DO i = 1, nbgp
1750	topo(:,:,nxr+i) = topo(:,:,nxr)
1751	ENDDO
1752	ENDIF
1753	ENDIF
1754
1755	CASE ( 'tunnel' )
1756
1757	!
1758	!-- Tunnel height
1759	IF ( tunnel_height == 9999999.9_wp ) THEN
1760	th = zw( INT( 0.2 * nz) )
1761	ELSE
1762	th = tunnel_height
1763	ENDIF
1764	!
1765	!-- Tunnel-wall depth
1766	IF ( tunnel_wall_depth == 9999999.9_wp ) THEN
1767	td = MAX ( dx, dy, dz )
1768	ELSE
1769	td = tunnel_wall_depth
1770	ENDIF
1771	!
1772	!-- Check for tunnel width
1773	IF ( tunnel_width_x == 9999999.9_wp .AND. &
1774	tunnel_width_y == 9999999.9_wp ) THEN
1775	message_string = 'No tunnel width is given. '
1776	CALL message( 'init_grid', 'PA0280', 1, 2, 0, 6, 0 )
1777	ENDIF
1778	IF ( tunnel_width_x /= 9999999.9_wp .AND. &
1779	tunnel_width_y /= 9999999.9_wp ) THEN
1780	message_string = 'Inconsistent tunnel parameters:' // &
1781	'tunnel can only be oriented' // &
1782	'either in x- or in y-direction.'
1783	CALL message( 'init_grid', 'PA0281', 1, 2, 0, 6, 0 )
1784	ENDIF
1785	!
1786	!-- Tunnel axis along y
1787	IF ( tunnel_width_x /= 9999999.9_wp ) THEN
1788	IF ( tunnel_width_x > ( nx + 1 ) * dx ) THEN
1789	message_string = 'Tunnel width too large'
1790	CALL message( 'init_grid', 'PA0282', 1, 2, 0, 6, 0 )
1791	ENDIF
1792
1793	txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_width_x * 0.5_wp )
1794	txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_width_x * 0.5_wp )
1795	txs_in = INT( ( nx + 1 ) * 0.5_wp * dx - &
1796	( tunnel_width_x * 0.5_wp - td ) )
1797	txe_in = INT( ( nx + 1 ) * 0.5_wp * dx + &
1798	( tunnel_width_x * 0.5_wp - td ) )
1799
1800	tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_length * 0.5_wp )
1801	tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_length * 0.5_wp )
1802	tys_in = tys_out
1803	tye_in = tye_out
1804	ENDIF
1805	IF ( tunnel_width_x /= 9999999.9_wp .AND. &
1806	tunnel_width_x - 2.0_wp * td <= 2.0_wp * dx ) &
1807	THEN
1808	message_string = 'Tunnel width too small'
1809	CALL message( 'init_grid', 'PA0175', 1, 2, 0, 6, 0 )
1810	ENDIF
1811	IF ( tunnel_width_y /= 9999999.9_wp .AND. &
1812	tunnel_width_y - 2.0_wp * td <= 2.0_wp * dy ) &
1813	THEN
1814	message_string = 'Tunnel width too small'
1815	CALL message( 'init_grid', 'PA0455', 1, 2, 0, 6, 0 )
1816	ENDIF
1817	!
1818	!-- Tunnel axis along x
1819	IF ( tunnel_width_y /= 9999999.9_wp ) THEN
1820	IF ( tunnel_width_y > ( ny + 1 ) * dy ) THEN
1821	message_string = 'Tunnel width too large'
1822	CALL message( 'init_grid', 'PA0456', 1, 2, 0, 6, 0 )
1823	ENDIF
1824
1825	txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_length * 0.5_wp )
1826	txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_length * 0.5_wp )
1827	txs_in = txs_out
1828	txe_in = txe_out
1829
1830	tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_width_y * 0.5_wp )
1831	tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_width_y * 0.5_wp )
1832	tys_in = INT( ( ny + 1 ) * 0.5_wp * dy - &
1833	( tunnel_width_y * 0.5_wp - td ) )
1834	tye_in = INT( ( ny + 1 ) * 0.5_wp * dy + &
1835	( tunnel_width_y * 0.5_wp - td ) )
1836	ENDIF
1837
1838	topo = 0
1839	DO i = nxl, nxr
1840	DO j = nys, nyn
1841	!
1842	!-- Use heaviside function to model outer tunnel surface
1843	hv_out = th * 0.5_wp * &
1844	( ( SIGN( 1.0_wp, i * dx - txs_out ) + 1.0_wp ) &
1845	- ( SIGN( 1.0_wp, i * dx - txe_out ) + 1.0_wp ) )
1846
1847	hv_out = hv_out * 0.5_wp * &
1848	( ( SIGN( 1.0_wp, j * dy - tys_out ) + 1.0_wp ) &
1849	- ( SIGN( 1.0_wp, j * dy - tye_out ) + 1.0_wp ) )
1850	!
1851	!-- Use heaviside function to model inner tunnel surface
1852	hv_in = ( th - td ) * 0.5_wp * &
1853	( ( SIGN( 1.0_wp, i * dx - txs_in ) + 1.0_wp ) &
1854	- ( SIGN( 1.0_wp, i * dx - txe_in ) + 1.0_wp ) )
1855
1856	hv_in = hv_in * 0.5_wp * &
1857	( ( SIGN( 1.0_wp, j * dy - tys_in ) + 1.0_wp ) &
1858	- ( SIGN( 1.0_wp, j * dy - tye_in ) + 1.0_wp ) )
1859	!
1860	!-- Set flags at x-y-positions without any tunnel surface
1861	IF ( hv_out - hv_in == 0.0_wp ) THEN
1862	topo(nzb+1:nzt+1,j,i) = IBSET( topo(nzb+1:nzt+1,j,i), 0 )
1863	!
1864	!-- Set flags at x-y-positions with tunnel surfaces
1865	ELSE
1866	DO k = nzb + 1, nzt + 1
1867	!
1868	!-- Inner tunnel
1869	IF ( hv_out - hv_in == th ) THEN
1870	IF ( zw(k) <= hv_out ) THEN
1871	topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
1872	ELSE
1873	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
1874	ENDIF
1875	ENDIF
1876	!
1877	!-- Lateral tunnel walls
1878	IF ( hv_out - hv_in == td ) THEN
1879	IF ( zw(k) <= hv_in ) THEN
1880	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
1881	ELSEIF ( zw(k) > hv_in .AND. zw(k) <= hv_out ) THEN
1882	topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
1883	ELSEIF ( zw(k) > hv_out ) THEN
1884	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
1885	ENDIF
1886	ENDIF
1887	ENDDO
1888	ENDIF
1889	ENDDO
1890	ENDDO
1891
1892	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1893	!
1894	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
1895	!-- boundary conditions for topography.
1896	IF ( .NOT. bc_ns_cyc ) THEN
1897	IF ( nys == 0 ) THEN
1898	DO i = 1, nbgp
1899	topo(:,nys-i,:) = topo(:,nys,:)
1900	ENDDO
1901	ENDIF
1902	IF ( nyn == ny ) THEN
1903	DO i = 1, nbgp
1904	topo(:,nyn+i,:) = topo(:,nyn,:)
1905	ENDDO
1906	ENDIF
1907	ENDIF
1908	IF ( .NOT. bc_lr_cyc ) THEN
1909	IF ( nxl == 0 ) THEN
1910	DO i = 1, nbgp
1911	topo(:,:,nxl-i) = topo(:,:,nxl)
1912	ENDDO
1913	ENDIF
1914	IF ( nxr == nx ) THEN
1915	DO i = 1, nbgp
1916	topo(:,:,nxr+i) = topo(:,:,nxr)
1917	ENDDO
1918	ENDIF
1919	ENDIF
1920
1921	CASE ( 'read_from_file' )
1922	!
1923	!-- Note, topography information have been already read.
1924	!-- If required, further process topography, i.e. reference buildings on
1925	!-- top of orography and set temporary 3D topography array, which is
1926	!-- used later to set grid flags. Calling of this rouinte is also
1927	!-- required in case of ASCII input, even though no distinction between
1928	!-- terrain- and building height is made in this case.
1929	CALL process_topography( topo )
1930	!
1931	!-- Filter holes resolved by only one grid-point
1932	CALL filter_topography( topo )
1933	!
1934	!-- Exchange ghost-points, as well as add cyclic or Neumann boundary
1935	!-- conditions.
1936	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1937	!
1938	!-- Set lateral boundary conditions for topography on all ghost layers
1939	IF ( .NOT. bc_ns_cyc ) THEN
1940	IF ( nys == 0 ) THEN
1941	DO i = 1, nbgp
1942	topo(:,nys-i,:) = topo(:,nys,:)
1943	ENDDO
1944	ENDIF
1945	IF ( nyn == ny ) THEN
1946	DO i = 1, nbgp
1947	topo(:,nyn+i,:) = topo(:,nyn,:)
1948	ENDDO
1949	ENDIF
1950	ENDIF
1951
1952	IF ( .NOT. bc_lr_cyc ) THEN
1953	IF ( nxl == 0 ) THEN
1954	DO i = 1, nbgp
1955	topo(:,:,nxl-i) = topo(:,:,nxl)
1956	ENDDO
1957	ENDIF
1958	IF ( nxr == nx ) THEN
1959	DO i = 1, nbgp
1960	topo(:,:,nxr+i) = topo(:,:,nxr)
1961	ENDDO
1962	ENDIF
1963	ENDIF
1964
1965
1966	CASE DEFAULT
1967	!
1968	!-- The DEFAULT case is reached either if the parameter topography
1969	!-- contains a wrong character string or if the user has defined a special
1970	!-- case in the user interface. There, the subroutine user_init_grid
1971	!-- checks which of these two conditions applies.
1972	CALL user_init_grid( topo )
1973	CALL filter_topography( topo )
1974
1975	END SELECT
1976	!
1977	!-- Consistency checks and index array initialization are only required for
1978	!-- non-flat topography.
1979	IF ( TRIM( topography ) /= 'flat' ) THEN
1980	!
1981	!-- In case of non-flat topography, check whether the convention how to
1982	!-- define the topography grid has been set correctly, or whether the default
1983	!-- is applicable. If this is not possible, abort.
1984	IF ( TRIM( topography_grid_convention ) == ' ' ) THEN
1985	IF ( TRIM( topography ) /= 'single_building' .AND. &
1986	TRIM( topography ) /= 'single_street_canyon' .AND. &
1987	TRIM( topography ) /= 'tunnel' .AND. &
1988	TRIM( topography ) /= 'read_from_file') THEN
1989	!-- The default value is not applicable here, because it is only valid
1990	!-- for the two standard cases 'single_building' and 'read_from_file'
1991	!-- defined in init_grid.
1992	WRITE( message_string, * ) &
1993	'The value for "topography_grid_convention" ', &
1994	'is not set. Its default value is & only valid for ', &
1995	'"topography" = ''single_building'', ', &
1996	'''single_street_canyon'' & or ''read_from_file''.', &
1997	' & Choose ''cell_edge'' or ''cell_center''.'
1998	CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 )
1999	ELSE
2000	!-- The default value is applicable here.
2001	!-- Set convention according to topography.
2002	IF ( TRIM( topography ) == 'single_building' .OR. &
2003	TRIM( topography ) == 'single_street_canyon' ) THEN
2004	topography_grid_convention = 'cell_edge'
2005	ELSEIF ( TRIM( topography ) == 'read_from_file' .OR. &
2006	TRIM( topography ) == 'tunnel') THEN
2007	topography_grid_convention = 'cell_center'
2008	ENDIF
2009	ENDIF
2010	ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND. &
2011	TRIM( topography_grid_convention ) /= 'cell_center' ) THEN
2012	WRITE( message_string, * ) &
2013	'The value for "topography_grid_convention" is ', &
2014	'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
2015	CALL message( 'init_grid', 'PA0240', 1, 2, 0, 6, 0 )
2016	ENDIF
2017
2018
2019	IF ( topography_grid_convention == 'cell_edge' ) THEN
2020	!
2021	!-- The array nzb_local as defined using the 'cell_edge' convention
2022	!-- describes the actual total size of topography which is defined at the
2023	!-- cell edges where u=0 on the topography walls in x-direction and v=0
2024	!-- on the topography walls in y-direction. However, PALM uses individual
2025	!-- arrays nzb_u\|v\|w\|s_inner\|outer that are based on nzb_s_inner.
2026	!-- Therefore, the extent of topography in nzb_local is now reduced by
2027	!-- 1dx at the E topography walls and by 1dy at the N topography walls
2028	!-- to form the basis for nzb_s_inner.
2029	!-- Note, the reverse memory access (i-j instead of j-i) is absolutely
2030	!-- required at this point.
2031	DO j = nys+1, nyn+1
2032	DO i = nxl-1, nxr
2033	DO k = nzb, nzt+1
2034	IF ( BTEST( topo(k,j,i), 0 ) .OR. &
2035	BTEST( topo(k,j,i+1), 0 ) ) &
2036	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2037	ENDDO
2038	ENDDO
2039	ENDDO
2040	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2041
2042	DO i = nxl, nxr+1
2043	DO j = nys-1, nyn
2044	DO k = nzb, nzt+1
2045	IF ( BTEST( topo(k,j,i), 0 ) .OR. &
2046	BTEST( topo(k,j+1,i), 0 ) ) &
2047	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2048	ENDDO
2049	ENDDO
2050	ENDDO
2051	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2052
2053	ENDIF
2054	ENDIF
2055
2056
2057	END SUBROUTINE init_topo
2058
2059	SUBROUTINE set_topo_flags(topo)
2060
2061	USE control_parameters, &
2062	ONLY: bc_lr_cyc, bc_ns_cyc, constant_flux_layer, land_surface, &
2063	use_surface_fluxes, use_top_fluxes, urban_surface
2064
2065	USE indices, &
2066	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, &
2067	nzb, nzt, wall_flags_0
2068
2069	USE kinds
2070
2071	IMPLICIT NONE
2072
2073	INTEGER(iwp) :: i !< index variable along x
2074	INTEGER(iwp) :: j !< index variable along y
2075	INTEGER(iwp) :: k !< index variable along z
2076
2077	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo !< input array for 3D topography and dummy array for setting "outer"-flags
2078
2079	ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2080	wall_flags_0 = 0
2081	!
2082	!-- Set-up topography flags. First, set flags only for s, u, v and w-grid.
2083	!-- Further special flags will be set in following loops.
2084	DO i = nxl, nxr
2085	DO j = nys, nyn
2086	DO k = nzb, nzt+1
2087	!
2088	!-- scalar grid
2089	IF ( BTEST( topo(k,j,i), 0 ) ) &
2090	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
2091	!
2092	!-- u grid
2093	IF ( BTEST( topo(k,j,i), 0 ) .AND. &
2094	BTEST( topo(k,j,i-1), 0 ) ) &
2095	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
2096	!
2097	!-- v grid
2098	IF ( BTEST( topo(k,j,i), 0 ) .AND. &
2099	BTEST( topo(k,j-1,i), 0 ) ) &
2100	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
2101
2102	ENDDO
2103
2104	DO k = nzb, nzt
2105	!
2106	!-- w grid
2107	IF ( BTEST( topo(k,j,i), 0 ) .AND. &
2108	BTEST( topo(k+1,j,i), 0 ) ) &
2109	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
2110	ENDDO
2111	wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 3 )
2112
2113	ENDDO
2114	ENDDO
2115
2116	CALL exchange_horiz_int( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
2117	!
2118	!-- Set outer array for scalars to mask near-surface grid points in
2119	!-- production_e
2120	DO i = nxl, nxr
2121	DO j = nys, nyn
2122	DO k = nzb, nzt+1
2123	IF ( BTEST( wall_flags_0(k,j-1,i), 0 ) .AND. &
2124	BTEST( wall_flags_0(k,j+1,i), 0 ) .AND. &
2125	BTEST( wall_flags_0(k,j,i-1), 0 ) .AND. &
2126	BTEST( wall_flags_0(k,j-1,i-1), 0 ) .AND. &
2127	BTEST( wall_flags_0(k,j+1,i-1), 0 ) .AND. &
2128	BTEST( wall_flags_0(k,j-1,i+1), 0 ) .AND. &
2129	BTEST( wall_flags_0(k,j+1,i+1), 0 ) ) &
2130	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
2131	ENDDO
2132	ENDDO
2133	ENDDO
2134	!
2135	!-- Set further special flags
2136	DO i = nxl, nxr
2137	DO j = nys, nyn
2138	DO k = nzb, nzt+1
2139	!
2140	!-- scalar grid, former nzb_diff_s_inner.
2141	!-- Note, use this flag also to mask topography in diffusion_u and
2142	!-- diffusion_v along the vertical direction. In case of
2143	!-- use_surface_fluxes, fluxes are calculated via MOST, else, simple
2144	!-- gradient approach is applied. Please note, in case of u- and v-
2145	!-- diffuison, a small error is made at edges (on the east side for u,
2146	!-- at the north side for v), since topography on scalar grid point
2147	!-- is used instead of topography on u/v-grid. As number of topography grid
2148	!-- points on uv-grid is different than s-grid, different number of
2149	!-- surface elements would be required. In order to avoid this,
2150	!-- treat edges (u(k,j,i+1)) simply by a gradient approach, i.e. these
2151	!-- points are not masked within diffusion_u. Tests had shown that the
2152	!-- effect on the flow is negligible.
2153	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2154	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) &
2155	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
2156	ELSE
2157	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
2158	ENDIF
2159
2160	ENDDO
2161	!
2162	!-- Special flag to control vertical diffusion at model top - former
2163	!-- nzt_diff
2164	wall_flags_0(:,j,i) = IBSET( wall_flags_0(:,j,i), 9 )
2165	IF ( use_top_fluxes ) &
2166	wall_flags_0(nzt+1,j,i) = IBCLR( wall_flags_0(nzt+1,j,i), 9 )
2167
2168
2169	DO k = nzb+1, nzt
2170	!
2171	!-- Special flag on u grid, former nzb_u_inner + 1, required
2172	!-- for disturb_field and initialization. Do not disturb directly at
2173	!-- topography, as well as initialize u with zero one grid point outside
2174	!-- of topography.
2175	IF ( BTEST( wall_flags_0(k-1,j,i), 1 ) .AND. &
2176	BTEST( wall_flags_0(k,j,i), 1 ) .AND. &
2177	BTEST( wall_flags_0(k+1,j,i), 1 ) ) &
2178	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 )
2179	!
2180	!-- Special flag on v grid, former nzb_v_inner + 1, required
2181	!-- for disturb_field and initialization. Do not disturb directly at
2182	!-- topography, as well as initialize v with zero one grid point outside
2183	!-- of topography.
2184	IF ( BTEST( wall_flags_0(k-1,j,i), 2 ) .AND. &
2185	BTEST( wall_flags_0(k,j,i), 2 ) .AND. &
2186	BTEST( wall_flags_0(k+1,j,i), 2 ) ) &
2187	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
2188	!
2189	!-- Special flag on scalar grid, former nzb_s_inner+1. Used for
2190	!-- lpm_sgs_tke
2191	IF ( BTEST( wall_flags_0(k,j,i), 0 ) .AND. &
2192	BTEST( wall_flags_0(k-1,j,i), 0 ) .AND. &
2193	BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
2194	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 )
2195	!
2196	!-- Special flag on scalar grid, nzb_diff_s_outer - 1, required in
2197	!-- in production_e
2198	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2199	IF ( BTEST( wall_flags_0(k,j,i), 24 ) .AND. &
2200	BTEST( wall_flags_0(k-1,j,i), 24 ) .AND. &
2201	BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
2202	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 29 )
2203	ELSE
2204	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) &
2205	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 29 )
2206	ENDIF
2207	!
2208	!-- Special flag on scalar grid, nzb_diff_s_outer - 1, required in
2209	!-- in production_e
2210	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2211	IF ( BTEST( wall_flags_0(k,j,i), 0 ) .AND. &
2212	BTEST( wall_flags_0(k-1,j,i), 0 ) .AND. &
2213	BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
2214	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
2215	ELSE
2216	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) &
2217	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
2218	ENDIF
2219	ENDDO
2220	!
2221	!-- Flags indicating downward facing walls
2222	DO k = nzb+1, nzt
2223	!
2224	!-- Scalar grid
2225	IF ( BTEST( wall_flags_0(k-1,j,i), 0 ) .AND. &
2226	.NOT. BTEST( wall_flags_0(k,j,i), 0 ) ) &
2227	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
2228	!
2229	!-- Downward facing wall on u grid
2230	IF ( BTEST( wall_flags_0(k-1,j,i), 1 ) .AND. &
2231	.NOT. BTEST( wall_flags_0(k,j,i), 1 ) ) &
2232	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 )
2233	!
2234	!-- Downward facing wall on v grid
2235	IF ( BTEST( wall_flags_0(k-1,j,i), 2 ) .AND. &
2236	.NOT. BTEST( wall_flags_0(k,j,i), 2 ) ) &
2237	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 )
2238	!
2239	!-- Downward facing wall on w grid
2240	IF ( BTEST( wall_flags_0(k-1,j,i), 3 ) .AND. &
2241	.NOT. BTEST( wall_flags_0(k,j,i), 3 ) ) &
2242	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
2243	ENDDO
2244	!
2245	!-- Flags indicating upward facing walls
2246	DO k = nzb, nzt
2247	!
2248	!-- Upward facing wall on scalar grid
2249	IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND. &
2250	BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
2251	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
2252	!
2253	!-- Upward facing wall on u grid
2254	IF ( .NOT. BTEST( wall_flags_0(k,j,i), 1 ) .AND. &
2255	BTEST( wall_flags_0(k+1,j,i), 1 ) ) &
2256	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
2257
2258	!
2259	!-- Upward facing wall on v grid
2260	IF ( .NOT. BTEST( wall_flags_0(k,j,i), 2 ) .AND. &
2261	BTEST( wall_flags_0(k+1,j,i), 2 ) ) &
2262	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
2263
2264	!
2265	!-- Upward facing wall on w grid
2266	IF ( .NOT. BTEST( wall_flags_0(k,j,i), 3 ) .AND. &
2267	BTEST( wall_flags_0(k+1,j,i), 3 ) ) &
2268	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
2269	!
2270	!-- Special flag on scalar grid, former nzb_s_inner
2271	IF ( BTEST( wall_flags_0(k,j,i), 0 ) .OR. &
2272	BTEST( wall_flags_0(k,j,i), 12 ) .OR. &
2273	BTEST( wall_flags_0(k,j,i), 13 ) ) &
2274	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
2275	!
2276	!-- Special flag on scalar grid, nzb_diff_s_inner - 1, required for
2277	!-- flow_statistics
2278	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2279	IF ( BTEST( wall_flags_0(k,j,i), 0 ) .AND. &
2280	BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
2281	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
2282	ELSE
2283	IF ( BTEST( wall_flags_0(k,j,i), 22 ) ) &
2284	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
2285	ENDIF
2286
2287
2288	ENDDO
2289	wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 22 )
2290	wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 23 )
2291	ENDDO
2292	ENDDO
2293	!
2294	!-- Finally, set identification flags indicating natural terrain or buildings.
2295	!-- Natural terrain grid points.
2296	IF ( land_surface ) THEN
2297	DO i = nxl, nxr
2298	DO j = nys, nyn
2299	DO k = nzb, nzt+1
2300	!
2301	!-- Natural terrain grid point
2302	IF ( BTEST( topo(k,j,i), 1 ) ) &
2303	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
2304	ENDDO
2305	ENDDO
2306	ENDDO
2307	ENDIF
2308	!
2309	!-- Building grid points.
2310	IF ( urban_surface ) THEN
2311	DO i = nxl, nxr
2312	DO j = nys, nyn
2313	DO k = nzb, nzt+1
2314	IF ( BTEST( topo(k,j,i), 2 ) ) &
2315	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
2316	ENDDO
2317	ENDDO
2318	ENDDO
2319	ENDIF
2320	!
2321	!-- Exchange ghost points for wall flags
2322	CALL exchange_horiz_int( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
2323	!
2324	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
2325	!-- boundary conditions for topography.
2326	IF ( .NOT. bc_ns_cyc ) THEN
2327	IF ( nys == 0 ) THEN
2328	DO i = 1, nbgp
2329	wall_flags_0(:,nys-i,:) = wall_flags_0(:,nys,:)
2330	ENDDO
2331	ENDIF
2332	IF ( nyn == ny ) THEN
2333	DO i = 1, nbgp
2334	wall_flags_0(:,nyn+i,:) = wall_flags_0(:,nyn,:)
2335	ENDDO
2336	ENDIF
2337	ENDIF
2338	IF ( .NOT. bc_lr_cyc ) THEN
2339	IF ( nxl == 0 ) THEN
2340	DO i = 1, nbgp
2341	wall_flags_0(:,:,nxl-i) = wall_flags_0(:,:,nxl)
2342	ENDDO
2343	ENDIF
2344	IF ( nxr == nx ) THEN
2345	DO i = 1, nbgp
2346	wall_flags_0(:,:,nxr+i) = wall_flags_0(:,:,nxr)
2347	ENDDO
2348	ENDIF
2349	ENDIF
2350
2351
2352	END SUBROUTINE set_topo_flags
2353
2354
2355

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |