Home

Context Navigation

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

Last change on this file since 4626 was 4601, checked in by suehring, 4 years ago
Chemistry: Bugfix in variable name separation in profile-output initialization; Bugfix in counting the number of chemistry profiles; Surface-data output: Minor simplification in name creation for IO variables in restart files; init_grid: minor formatting adjustments
Property svn:keywords set to `Id`
File size: 118.1 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-2020 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: init_grid.f90 4601 2020-07-14 12:06:09Z raasch $
27	! Minor formatting adjustments
28	!
29	! 4564 2020-06-12 14:03:36Z raasch
30	! Vertical nesting method of Huq et al. (2019) removed
31	!
32	! 4543 2020-05-20 14:12:22Z gronemeier
33	! Remove non-required check for canyon height
34	!
35	! 4507 2020-04-22 18:21:45Z gronemeier
36	! update origin_z with shifting height of orography (oro_min)
37	!
38	! 4457 2020-03-11 14:20:43Z raasch
39	! use statement for exchange horiz added,
40	! bugfix for call of exchange horiz 2d
41	!
42	! 4444 2020-03-05 15:59:50Z raasch
43	! bugfix: cpp-directives for serial mode added
44	!
45	! 4414 2020-02-19 20:16:04Z suehring
46	! - Remove deprecated topography arrays nzb_s_inner, nzb_u_inner, etc.
47	! - Move initialization of boundary conditions and multigrid into an extra
48	! module interface.
49	!
50	! 4386 2020-01-27 15:07:30Z Giersch
51	! Allocation statements, comments, naming of variables revised and _wp added to
52	! real type values
53	!
54	! 4360 2020-01-07 11:25:50Z suehring
55	! Revise error messages for generic tunnel setup.
56	!
57	! 4346 2019-12-18 11:55:56Z motisi
58	! Introduction of wall_flags_total_0, which currently sets bits based on static
59	! topography information used in wall_flags_static_0
60	!
61	! 4340 2019-12-16 08:17:03Z Giersch
62	! Topography closed channel flow with symmetric boundaries implemented
63	!
64	! 4329 2019-12-10 15:46:36Z motisi
65	! Renamed wall_flags_0 to wall_flags_static_0
66	!
67	! 4328 2019-12-09 18:53:04Z suehring
68	! Minor change in nzb_max computation. Commentation added.
69	!
70	! 4314 2019-11-29 10:29:20Z suehring
71	! Set additional topography flag 4 to mark topography grid points emerged
72	! from the filtering process.
73	!
74	! 4294 2019-11-13 18:34:16Z suehring
75	! Bugfix, always set bit 5 and 6 of wall_flags, indicating terrain- and
76	! building surfaces in all cases, in order to enable terrain-following output
77	! also when no land- or urban-surface model is applied.
78	!
79	! 4265 2019-10-15 16:16:24Z suehring
80	! Bugfix for last commit, exchange oro_max variable only when it is allocated
81	! (not necessarily the case when topography is input from ASCII file).
82	!
83	! 4245 2019-09-30 08:40:37Z pavelkrc
84	! Store oro_max (building z-offset) in 2D for building surfaces
85	!
86	! 4189 2019-08-26 16:19:38Z suehring
87	! - Add check for proper setting of namelist parameter topography
88	! - Set flag to indicate land surfaces in case no topography is provided
89	!
90	! 4182 2019-08-22 15:20:23Z scharf
91	! Corrected "Former revisions" section
92	!
93	! 4168 2019-08-16 13:50:17Z suehring
94	! Pre-calculate topography top index and store it on an array (replaces former
95	! functions get_topography_top_index)
96	!
97	! 4159 2019-08-15 13:31:35Z suehring
98	! Revision of topography processing. This was not consistent between 2D and 3D
99	! buildings.
100	!
101	! 4144 2019-08-06 09:11:47Z raasch
102	! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
103	!
104	! 4115 2019-07-24 12:50:49Z suehring
105	! Bugfix in setting near-surface flag 24, inidicating wall-bounded grid points
106	!
107	! 4110 2019-07-22 17:05:21Z suehring
108	! - Separate initialization of advection flags for momentum and scalars.
109	! - Change subroutine interface for ws_init_flags_scalar to pass boundary flags
110	!
111	! 4109 2019-07-22 17:00:34Z suehring
112	! Fix bad commit
113	!
114	! 3926 2019-04-23 12:56:42Z suehring
115	! Minor bugfix in building mapping when all building IDs in the model domain
116	! are missing
117	!
118	! 3857 2019-04-03 13:00:16Z knoop
119	! In projection of non-building 3D objects onto numerical grid remove
120	! dependency on building_type
121	!
122	! 3763 2019-02-25 17:33:49Z suehring
123	! Replace work-around for ghost point exchange of 1-byte arrays with specific
124	! routine as already done in other routines
125	!
126	! 3761 2019-02-25 15:31:42Z raasch
127	! unused variables removed
128	!
129	! 3661 2019-01-08 18:22:50Z suehring
130	! Remove setting of nzb_max to nzt at non-cyclic boundary PEs, instead,
131	! order degradation of advection scheme is handeled directly in advec_ws
132	!
133	! 3655 2019-01-07 16:51:22Z knoop
134	! Comment added
135	!
136	! Revision 1.1 1997/08/11 06:17:45 raasch
137	! Initial revision (Testversion)
138	!
139	!
140	! Description:
141	! -----------------------------------------------------------------------------!
142	!> Creating grid depending constants
143	!> @todo: Rearrange topo flag list
144	!> @todo: reference 3D buildings on top of orography is not tested and may need
145	!> further improvement for steep slopes
146	!> @todo: Use more advanced setting of building type at filled holes
147	!------------------------------------------------------------------------------!
148	SUBROUTINE init_grid
149
150	USE arrays_3d, &
151	ONLY: dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, x, xu, y, yv, zu, zw
152
153	USE control_parameters, &
154	ONLY: constant_flux_layer, dz, dz_max, dz_stretch_factor, &
155	dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,&
156	dz_stretch_level_end_index, dz_stretch_level_start_index, &
157	dz_stretch_level_start, ibc_uv_b, message_string, &
158	number_stretch_level_end, &
159	number_stretch_level_start, &
160	ocean_mode, &
161	psolver, &
162	symmetry_flag, &
163	topography, &
164	use_surface_fluxes
165
166	USE grid_variables, &
167	ONLY: ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner
168
169	USE indices, &
170	ONLY: nbgp, &
171	nx, &
172	nxl, &
173	nxlg, &
174	nxr, &
175	nxrg, &
176	ny, &
177	nyn, &
178	nyng, &
179	nys, &
180	nysg, &
181	nz, &
182	nzb, &
183	nzb_diff, &
184	nzb_max, &
185	nzt, &
186	topo_top_ind, &
187	topo_min_level
188
189	USE kinds
190
191	USE pegrid
192
193	IMPLICIT NONE
194
195	INTEGER(iwp) :: i !< index variable along x
196	INTEGER(iwp) :: j !< index variable along y
197	INTEGER(iwp) :: k !< index variable along z
198	INTEGER(iwp) :: k_top !< topography top index on local PE
199	INTEGER(iwp) :: n !< loop variable for stretching
200	INTEGER(iwp) :: number_dz !< number of user-specified dz values
201	INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height
202	INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height
203
204	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo !< input array for 3D topography and dummy array for setting "outer"-flags
205
206	REAL(wp) :: dz_level_end !< distance between calculated height level for u/v-grid and user-specified end level for stretching
207	REAL(wp) :: dz_stretched !< stretched vertical grid spacing
208
209	REAL(wp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end !< Array that contains all minimum heights where the stretching can end
210
211
212	!
213	!-- Calculation of horizontal array bounds including ghost layers
214	nxlg = nxl - nbgp
215	nxrg = nxr + nbgp
216	nysg = nys - nbgp
217	nyng = nyn + nbgp
218
219	!
220	!-- Allocate grid arrays
221	ALLOCATE( x(0:nx) )
222	ALLOCATE( xu(0:nx) )
223
224	DO i = 0, nx
225	xu(i) = i * dx
226	x(i) = i * dx + 0.5_wp * dx
227	ENDDO
228
229	ALLOCATE( y(0:ny) )
230	ALLOCATE( yv(0:ny) )
231
232	DO j = 0, ny
233	yv(j) = j * dy
234	y(j) = j * dy + 0.5_wp * dy
235	ENDDO
236
237	ALLOCATE( ddzu(1:nzt+1) )
238	ALLOCATE( ddzw(1:nzt+1) )
239	ALLOCATE( dd2zu(1:nzt) )
240	ALLOCATE( dzu(1:nzt+1) )
241	ALLOCATE( dzw(1:nzt+1) )
242	ALLOCATE( zu(nzb:nzt+1) )
243	ALLOCATE( zw(nzb:nzt+1) )
244
245	!
246	!-- For constructing an appropriate grid, the vertical grid spacing dz has to
247	!-- be specified with a non-negative value in the parameter file
248	IF ( dz(1) == -1.0_wp ) THEN
249	message_string = 'missing dz'
250	CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
251	ELSEIF ( dz(1) <= 0.0_wp ) THEN
252	WRITE( message_string, * ) 'dz=',dz(1),' <= 0.0'
253	CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
254	ENDIF
255
256	!
257	!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
258	!-- if it was set by the user
259	IF ( dz_stretch_level /= -9999999.9_wp ) THEN
260	dz_stretch_level_start(1) = dz_stretch_level
261	ENDIF
262
263	!
264	!-- Determine number of dz values and stretching levels specified by the
265	!-- user to allow right controlling of the stretching mechanism and to
266	!-- perform error checks. The additional requirement that dz /= dz_max
267	!-- for counting number of user-specified dz values is necessary. Otherwise
268	!-- restarts would abort if the old stretching mechanism with dz_stretch_level
269	!-- is used (Attention: The user is not allowed to specify a dz value equal
270	!-- to the default of dz_max = 999.0).
271	number_dz = COUNT( dz /= -1.0_wp .AND. dz /= dz_max)
272	number_stretch_level_start = COUNT( dz_stretch_level_start /= &
273	-9999999.9_wp )
274	number_stretch_level_end = COUNT( dz_stretch_level_end /= &
275	9999999.9_wp )
276
277	!
278	!-- The number of specified end levels +1 has to be the same as the number
279	!-- of specified dz values
280	IF ( number_dz /= number_stretch_level_end + 1 ) THEN
281	WRITE( message_string, * ) 'The number of values for dz = ', &
282	number_dz, 'has to be the same as& ', &
283	'the number of values for ', &
284	'dz_stretch_level_end + 1 = ', &
285	number_stretch_level_end+1
286	CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 )
287	ENDIF
288
289	!
290	!-- The number of specified start levels has to be the same or one less than
291	!-- the number of specified dz values
292	IF ( number_dz /= number_stretch_level_start + 1 .AND. &
293	number_dz /= number_stretch_level_start ) THEN
294	WRITE( message_string, * ) 'The number of values for dz = ', &
295	number_dz, 'has to be the same as or one ', &
296	'more than& the number of values for ', &
297	'dz_stretch_level_start = ', &
298	number_stretch_level_start
299	CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
300	ENDIF
301
302	!-- The number of specified start levels has to be the same or one more than
303	!-- the number of specified end levels
304	IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. &
305	number_stretch_level_start /= number_stretch_level_end ) THEN
306	WRITE( message_string, * ) 'The number of values for ', &
307	'dz_stretch_level_start = ', &
308	dz_stretch_level_start, 'has to be the ', &
309	'same or one more than& the number of ', &
310	'values for dz_stretch_level_end = ', &
311	number_stretch_level_end
312	CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 )
313	ENDIF
314
315	!
316	!-- Initialize dz for the free atmosphere with the value of dz_max
317	IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND. &
318	number_stretch_level_start /= 0 ) THEN
319	dz(number_stretch_level_start+1) = dz_max
320	ENDIF
321
322	!
323	!-- Initialize the stretching factor if (infinitely) stretching in the free
324	!-- atmosphere is desired (dz_stretch_level_end was not specified for the
325	!-- free atmosphere)
326	IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN
327	dz_stretch_factor_array(number_stretch_level_start) = &
328	dz_stretch_factor
329	ENDIF
330
331	!
332	!-- Allocation of arrays for stretching
333	ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
334
335	!
336	!-- Define the vertical grid levels. Start with atmosphere branch
337	IF ( .NOT. ocean_mode ) THEN
338
339	!
340	!-- The stretching region has to be large enough to allow for a smooth
341	!-- transition between two different grid spacings. The number 4 is an
342	!-- empirical value
343	DO n = 1, number_stretch_level_start
344	min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + &
345	4 * MAX( dz(n),dz(n+1) )
346	ENDDO
347
348	IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > &
349	dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
350	message_string= 'Each dz_stretch_level_end has to be larger ' // &
351	'than its corresponding value for &' // &
352	'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '// &
353	'to allow for smooth grid stretching'
354	CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
355	ENDIF
356
357	!
358	!-- Stretching must not be applied within the surface layer
359	!-- (first two grid points). For the default case dz_stretch_level_start
360	!-- is negative. Therefore the absolut value is checked here.
361	IF ( ANY( ABS( dz_stretch_level_start ) <= dz(1) * 1.5_wp ) ) THEN
362	WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
363	'larger than ', dz(1) * 1.5
364	CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
365	ENDIF
366
367	!
368	!-- The stretching has to start and end on a grid level. Therefore
369	!-- user-specified values are mapped to the next lowest level. The
370	!-- calculation of the first level is realized differently just because of
371	!-- historical reasons (the advanced/new stretching mechanism was realized
372	!-- in a way that results don't change if the old parameters
373	!-- dz_stretch_level, dz_stretch_factor and dz_max are used)
374	IF ( number_stretch_level_start /= 0 ) THEN
375	dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - &
376	dz(1)/2.0) / dz(1) ) &
377	* dz(1) + dz(1)/2.0
378	ENDIF
379
380	IF ( number_stretch_level_start > 1 ) THEN
381	DO n = 2, number_stretch_level_start
382	dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / &
383	dz(n) ) * dz(n)
384	ENDDO
385	ENDIF
386
387	IF ( number_stretch_level_end /= 0 ) THEN
388	DO n = 1, number_stretch_level_end
389	dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / &
390	dz(n+1) ) * dz(n+1)
391	ENDDO
392	ENDIF
393
394	!
395	!-- Determine stretching factor if necessary
396	IF ( number_stretch_level_end >= 1 ) THEN
397	CALL calculate_stretching_factor( number_stretch_level_end )
398	ENDIF
399
400	!
401	!-- Grid for atmosphere with surface at z=0 (k=0, w-grid).
402	!-- First compute the u- and v-levels. In case of dirichlet bc for u and v
403	!-- the first u/v- and w-level (k=0) are defined at same height (z=0).
404	!-- The second u-level (k=1) corresponds to the top of the
405	!-- surface layer. In case of symmetric boundaries (closed channel flow),
406	!-- the first grid point is always at z=0.
407	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR. &
408	topography == 'closed_channel' ) THEN
409	zu(0) = 0.0_wp
410	ELSE
411	zu(0) = - dz(1) * 0.5_wp
412	ENDIF
413
414	zu(1) = dz(1) * 0.5_wp
415
416	!
417	!-- Determine u and v height levels considering the possibility of grid
418	!-- stretching in several heights.
419	n = 1
420	dz_stretch_level_start_index = nzt+1
421	dz_stretch_level_end_index = nzt+1
422	dz_stretched = dz(1)
423
424	!-- The default value of dz_stretch_level_start is negative, thus the first
425	!-- condition is true even if no stretching shall be applied. Hence, the
426	!-- second condition is also necessary.
427	DO k = 2, nzt+1-symmetry_flag
428	IF ( dz_stretch_level_start(n) <= zu(k-1) .AND. &
429	dz_stretch_level_start(n) /= -9999999.9_wp ) THEN
430	dz_stretched = dz_stretched * dz_stretch_factor_array(n)
431
432	IF ( dz(n) > dz(n+1) ) THEN
433	dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
434	ELSE
435	dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
436	ENDIF
437
438	IF ( dz_stretch_level_start_index(n) == nzt+1 ) &
439	dz_stretch_level_start_index(n) = k-1
440
441	ENDIF
442
443	zu(k) = zu(k-1) + dz_stretched
444
445	!
446	!-- Make sure that the stretching ends exactly at dz_stretch_level_end
447	dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) )
448
449	IF ( dz_level_end < dz(n+1)/3.0 ) THEN
450	zu(k) = dz_stretch_level_end(n)
451	dz_stretched = dz(n+1)
452	dz_stretch_level_end_index(n) = k
453	n = n + 1
454	ENDIF
455	ENDDO
456
457	!
458	!-- If a closed channel flow is simulated, make sure that grid structure is
459	!-- the same for both bottom and top boundary. (Hint: Using a different dz
460	!-- at the bottom and at the top makes no sense due to symmetric boundaries
461	!-- where dz should be equal. Therefore, different dz at the bottom and top
462	!-- causes an abort (see check_parameters).)
463	IF ( topography == 'closed_channel' ) THEN
464	zu(nzt+1) = zu(nzt) + dz(1) * 0.5_wp
465	ENDIF
466
467	!
468	!-- Compute the w-levels. They are always staggered half-way between the
469	!-- corresponding u-levels. In case of dirichlet bc for u and v at the
470	!-- ground the first u- and w-level (k=0) are defined at same height (z=0).
471	!-- Per default, the top w-level is extrapolated linearly. In case of
472	!-- a closed channel flow, zu(nzt+1) and zw(nzt) must be set explicitely.
473	!-- (Hint: Using a different dz at the bottom and at the top makes no sense
474	!-- due to symmetric boundaries where dz should be equal. Therefore,
475	!-- different dz at the bottom and top causes an abort (see
476	!-- check_parameters).)
477	zw(0) = 0.0_wp
478	DO k = 1, nzt-symmetry_flag
479	zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
480	ENDDO
481	IF ( topography == 'closed_channel' ) THEN
482	zw(nzt) = zw(nzt-1) + dz(1)
483	zw(nzt+1) = zw(nzt) + dz(1)
484	ELSE
485	zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
486	ENDIF
487
488	ELSE !ocean branch
489
490	!
491	!-- The stretching region has to be large enough to allow for a smooth
492	!-- transition between two different grid spacings. The number 4 is an
493	!-- empirical value
494	DO n = 1, number_stretch_level_start
495	min_dz_stretch_level_end(n) = dz_stretch_level_start(n) - &
496	4 * MAX( dz(n),dz(n+1) )
497	ENDDO
498
499	IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) < &
500	dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
501	message_string= 'Each dz_stretch_level_end has to be less ' // &
502	'than its corresponding value for &' // &
503	'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//&
504	'to allow for smooth grid stretching'
505	CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
506	ENDIF
507
508	!
509	!-- Stretching must not be applied close to the surface (last two grid
510	!-- points). For the default case dz_stretch_level_start is negative.
511	IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) ) THEN
512	WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
513	'less than ', -dz(1) * 1.5
514	CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
515	ENDIF
516
517	!
518	!-- The stretching has to start and end on a grid level. Therefore
519	!-- user-specified values are mapped to the next highest level. The
520	!-- calculation of the first level is realized differently just because of
521	!-- historical reasons (the advanced/new stretching mechanism was realized
522	!-- in a way that results don't change if the old parameters
523	!-- dz_stretch_level, dz_stretch_factor and dz_max are used)
524	IF ( number_stretch_level_start /= 0 ) THEN
525	dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) + &
526	dz(1)/2.0) / dz(1) ) &
527	* dz(1) - dz(1)/2.0
528	ENDIF
529
530	IF ( number_stretch_level_start > 1 ) THEN
531	DO n = 2, number_stretch_level_start
532	dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / &
533	dz(n) ) * dz(n)
534	ENDDO
535	ENDIF
536
537	IF ( number_stretch_level_end /= 0 ) THEN
538	DO n = 1, number_stretch_level_end
539	dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / &
540	dz(n+1) ) * dz(n+1)
541	ENDDO
542	ENDIF
543
544	!
545	!-- Determine stretching factor if necessary
546	IF ( number_stretch_level_end >= 1 ) THEN
547	CALL calculate_stretching_factor( number_stretch_level_end )
548	ENDIF
549
550	!
551	!-- Grid for ocean with free water surface is at k=nzt (w-grid).
552	!-- In case of neumann bc at the ground the first first u-level (k=0) lies
553	!-- below the first w-level (k=0). In case of dirichlet bc the first u- and
554	!-- w-level are defined at same height, but staggered from the second level.
555	!-- The second u-level (k=1) corresponds to the top of the surface layer.
556	!-- z values are negative starting from z=0 (surface)
557	zu(nzt+1) = dz(1) * 0.5_wp
558	zu(nzt) = - dz(1) * 0.5_wp
559
560	!
561	!-- Determine u and v height levels considering the possibility of grid
562	!-- stretching in several heights.
563	n = 1
564	dz_stretch_level_start_index = 0
565	dz_stretch_level_end_index = 0
566	dz_stretched = dz(1)
567
568	DO k = nzt-1, 0, -1
569
570	IF ( dz_stretch_level_start(n) >= zu(k+1) ) THEN
571	dz_stretched = dz_stretched * dz_stretch_factor_array(n)
572
573	IF ( dz(n) > dz(n+1) ) THEN
574	dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
575	ELSE
576	dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
577	ENDIF
578
579	IF ( dz_stretch_level_start_index(n) == 0 ) &
580	dz_stretch_level_start_index(n) = k+1
581
582	ENDIF
583
584	zu(k) = zu(k+1) - dz_stretched
585
586	!
587	!-- Make sure that the stretching ends exactly at dz_stretch_level_end
588	dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) )
589
590	IF ( dz_level_end < dz(n+1)/3.0 ) THEN
591	zu(k) = dz_stretch_level_end(n)
592	dz_stretched = dz(n+1)
593	dz_stretch_level_end_index(n) = k
594	n = n + 1
595	ENDIF
596	ENDDO
597
598	!
599	!-- Compute the w-levels. They are always staggered half-way between the
600	!-- corresponding u-levels, except in case of dirichlet bc for u and v
601	!-- at the ground. In this case the first u- and w-level are defined at
602	!-- same height. The top w-level (nzt+1) is not used but set for
603	!-- consistency, since w and all scalar variables are defined up tp nzt+1.
604	zw(nzt+1) = dz(1)
605	zw(nzt) = 0.0_wp
606	DO k = 0, nzt
607	zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
608	ENDDO
609
610	!
611	!-- In case of dirichlet bc for u and v the first u- and w-level are defined
612	!-- at same height.
613	IF ( ibc_uv_b == 0 ) THEN
614	zu(0) = zw(0)
615	ENDIF
616
617	ENDIF !End of defining the vertical grid levels
618
619	!
620	!-- Compute grid lengths.
621	DO k = 1, nzt+1
622	dzu(k) = zu(k) - zu(k-1)
623	ddzu(k) = 1.0_wp / dzu(k)
624	dzw(k) = zw(k) - zw(k-1)
625	ddzw(k) = 1.0_wp / dzw(k)
626	ENDDO
627
628	DO k = 1, nzt
629	dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
630	ENDDO
631
632	!
633	!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
634	!-- everywhere. For the actual grid, the grid spacing at the lowest level
635	!-- is only dz/2, but should be dz. Therefore, an additional array
636	!-- containing with appropriate grid information is created for these
637	!-- solvers.
638	IF ( psolver(1:9) /= 'multigrid' ) THEN
639	ALLOCATE( ddzu_pres(1:nzt+1) )
640	ddzu_pres = ddzu
641	ddzu_pres(1) = ddzu_pres(2) ! change for lowest level
642	ENDIF
643
644	!
645	!-- Compute the reciprocal values of the horizontal grid lengths.
646	ddx = 1.0_wp / dx
647	ddy = 1.0_wp / dy
648	dx2 = dx * dx
649	dy2 = dy * dy
650	ddx2 = 1.0_wp / dx2
651	ddy2 = 1.0_wp / dy2
652
653	!
654	!-- Allocate 3D array to set topography
655	ALLOCATE( topo(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
656	topo = 0
657	!
658	!-- Initialize topography by generic topography or read topography from file.
659	CALL init_topo( topo )
660	!
661	!-- Set flags to mask topography on the grid.
662	CALL set_topo_flags( topo )
663
664	!
665	!-- Determine the maximum level of topography. It is used for
666	!-- steering the degradation of order of the applied advection scheme,
667	!-- as well in the lpm.
668	k_top = 0
669	DO i = nxl, nxr
670	DO j = nys, nyn
671	DO k = nzb, nzt + 1
672	k_top = MAX( k_top, MERGE( k, 0, .NOT. BTEST( topo(k,j,i), 0 ) ) )
673	ENDDO
674	ENDDO
675	ENDDO
676	#if defined( __parallel )
677	CALL MPI_ALLREDUCE( k_top, nzb_max, 1, MPI_INTEGER, &
678	MPI_MAX, comm2d, ierr )
679	#else
680	nzb_max = k_top
681	#endif
682	!
683	!-- Increment nzb_max by 1 in order to allow for proper diverengence correction.
684	!-- Further, in case topography extents up to the model top, limit to nzt.
685	nzb_max = MIN( nzb_max+1, nzt )
686	!
687	!-- Determine minimum index of topography. Usually, this will be nzb. In case
688	!-- there is elevated topography, however, the lowest topography will be higher.
689	!-- This index is e.g. used to calculate mean first-grid point atmosphere
690	!-- temperature, surface pressure and density, etc. .
691	topo_min_level = 0
692	#if defined( __parallel )
693	CALL MPI_ALLREDUCE( MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ), &
694	topo_min_level, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
695	#else
696	topo_min_level = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
697	#endif
698
699	!
700	!-- Check topography for consistency with model domain. Therefore, use
701	!-- maximum and minium topography-top indices. Note, minimum topography top
702	!-- index is already calculated.
703	IF ( TRIM( topography ) /= 'flat' ) THEN
704	#if defined( __parallel )
705	CALL MPI_ALLREDUCE( MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ), &
706	nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
707	#else
708	nzb_local_max = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
709	#endif
710	nzb_local_min = topo_min_level
711	!
712	!-- Consistency checks
713	IF ( nzb_local_min < 0 .OR. nzb_local_max > nz + 1 ) THEN
714	WRITE( message_string, * ) 'nzb_local values are outside the', &
715	' model domain', &
716	'&MINVAL( nzb_local ) = ', nzb_local_min, &
717	'&MAXVAL( nzb_local ) = ', nzb_local_max
718	CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
719	ENDIF
720	ENDIF
721	!
722	!-- Define vertical gridpoint from (or to) which on the usual finite difference
723	!-- form (which does not use surface fluxes) is applied
724	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
725	nzb_diff = nzb + 2
726	ELSE
727	nzb_diff = nzb + 1
728	ENDIF
729
730	IF ( TRIM( topography ) /= 'flat' ) THEN
731	!
732	!-- Allocate and set the arrays containing the topography height (for output
733	!-- reasons only).
734	IF ( nxr == nx .AND. nyn /= ny ) THEN
735	ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn), &
736	zw_w_inner(nxl:nxr+1,nys:nyn) )
737	ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
738	ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1), &
739	zw_w_inner(nxl:nxr,nys:nyn+1) )
740	ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
741	ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1), &
742	zw_w_inner(nxl:nxr+1,nys:nyn+1) )
743	ELSE
744	ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn), &
745	zw_w_inner(nxl:nxr,nys:nyn) )
746	ENDIF
747
748	zu_s_inner = 0.0_wp
749	zw_w_inner = 0.0_wp
750	!
751	!-- Determine local topography height on scalar and w-grid. Note, setting
752	!-- lateral boundary values is not necessary, realized via wall_flags_static_0
753	!-- array. Further, please note that loop bounds are different from
754	!-- nxl to nxr and nys to nyn on south and right model boundary, hence,
755	!-- use intrinsic lbound and ubound functions to infer array bounds.
756	DO i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
757	DO j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
758	!
759	!-- Topography height on scalar grid. Therefore, determine index of
760	!-- upward-facing surface element on scalar grid.
761	zu_s_inner(i,j) = zu(topo_top_ind(j,i,0))
762	!
763	!-- Topography height on w grid. Therefore, determine index of
764	!-- upward-facing surface element on w grid.
765	zw_w_inner(i,j) = zw(topo_top_ind(j,i,3))
766	ENDDO
767	ENDDO
768	ENDIF
769
770	END SUBROUTINE init_grid
771
772
773	! Description:
774	! -----------------------------------------------------------------------------!
775	!> Calculation of the stretching factor through an iterative method. Ideas were
776	!> taken from the paper "Regional stretched grid generation and its application
777	!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
778	!> system of equations has two variables (r,l) but four requirements
779	!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
780	!> results into an overdetermined system.
781	!------------------------------------------------------------------------------!
782	SUBROUTINE calculate_stretching_factor( number_end )
783
784	USE control_parameters, &
785	ONLY: dz, dz_stretch_factor_array, &
786	dz_stretch_level_end, dz_stretch_level_start, message_string
787
788	USE kinds
789
790	IMPLICIT NONE
791
792	INTEGER(iwp) :: iterations !< number of iterations until stretch_factor_lower/upper_limit is reached
793	INTEGER(iwp) :: l_rounded !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
794	INTEGER(iwp) :: n !< loop variable for stretching
795
796	INTEGER(iwp), INTENT(IN) :: number_end !< number of user-specified end levels for stretching
797
798	REAL(wp) :: delta_l !< absolute difference between l and l_rounded
799	REAL(wp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_2
800	REAL(wp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
801	REAL(wp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration
802	REAL(wp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
803	REAL(wp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
804	REAL(wp) :: numerator !< numerator of the quotient
805	REAL(wp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly
806	REAL(wp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
807
808	REAL(wp) :: dz_stretch_factor_array_2(9) = 1.08_wp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
809
810	REAL(wp), PARAMETER :: stretch_factor_interval = 1.0E-06_wp !< interval for sampling possible stretching factors
811	REAL(wp), PARAMETER :: stretch_factor_lower_limit = 0.88_wp !< lowest possible stretching factor
812	REAL(wp), PARAMETER :: stretch_factor_upper_limit = 1.12_wp !< highest possible stretching factor
813
814
815	l = 0
816	DO n = 1, number_end
817
818	iterations = 1
819	stretch_factor_1 = 1.0_wp
820	stretch_factor_2 = 1.0_wp
821	delta_total_old = 1.0_wp
822
823	!
824	!-- First branch for stretching from rough to fine
825	IF ( dz(n) > dz(n+1) ) THEN
826	DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
827
828	stretch_factor_1 = 1.0_wp - iterations * stretch_factor_interval
829	distance = ABS( dz_stretch_level_end(n) - &
830	dz_stretch_level_start(n) )
831	numerator = distance*stretch_factor_1/dz(n) + &
832	stretch_factor_1 - distance/dz(n)
833
834	IF ( numerator > 0.0_wp ) THEN
835	l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
836	l_rounded = NINT( l )
837	delta_l = ABS( l_rounded - l ) / l
838	ENDIF
839
840	stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
841
842	delta_stretch_factor = ABS( stretch_factor_1 - &
843	stretch_factor_2 ) / &
844	stretch_factor_2
845
846	delta_total_new = delta_l + delta_stretch_factor
847
848	!
849	!-- stretch_factor_1 is taken to guarantee that the stretching
850	!-- procedure ends as close as possible to dz_stretch_level_end.
851	!-- stretch_factor_2 would guarantee that the stretched dz(n) is
852	!-- equal to dz(n+1) after l_rounded grid levels.
853	IF (delta_total_new < delta_total_old) THEN
854	dz_stretch_factor_array(n) = stretch_factor_1
855	dz_stretch_factor_array_2(n) = stretch_factor_2
856	delta_total_old = delta_total_new
857	ENDIF
858
859	iterations = iterations + 1
860
861	ENDDO
862
863	!
864	!-- Second branch for stretching from fine to rough
865	ELSEIF ( dz(n) < dz(n+1) ) THEN
866	DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
867
868	stretch_factor_1 = 1.0_wp + iterations * stretch_factor_interval
869	distance = ABS( dz_stretch_level_end(n) - &
870	dz_stretch_level_start(n) )
871	numerator = distance*stretch_factor_1/dz(n) + &
872	stretch_factor_1 - distance/dz(n)
873
874	l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
875	l_rounded = NINT( l )
876	delta_l = ABS( l_rounded - l ) / l
877
878	stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
879
880	delta_stretch_factor = ABS( stretch_factor_1 - &
881	stretch_factor_2 ) / &
882	stretch_factor_2
883
884	delta_total_new = delta_l + delta_stretch_factor
885
886	!
887	!-- stretch_factor_1 is taken to guarantee that the stretching
888	!-- procedure ends as close as possible to dz_stretch_level_end.
889	!-- stretch_factor_2 would guarantee that the stretched dz(n) is
890	!-- equal to dz(n+1) after l_rounded grid levels.
891	IF (delta_total_new < delta_total_old) THEN
892	dz_stretch_factor_array(n) = stretch_factor_1
893	dz_stretch_factor_array_2(n) = stretch_factor_2
894	delta_total_old = delta_total_new
895	ENDIF
896
897	iterations = iterations + 1
898	ENDDO
899
900	ELSE
901	message_string= 'Two adjacent values of dz must be different'
902	CALL message( 'init_grid', 'PA0228', 1, 2, 0, 6, 0 )
903
904	ENDIF
905
906	!
907	!-- Check if also the second stretching factor fits into the allowed
908	!-- interval. If not, print a warning for the user.
909	IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR. &
910	dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
911	WRITE( message_string, * ) 'stretch_factor_2 = ', &
912	dz_stretch_factor_array_2(n), ' which is',&
913	' responsible for exactly reaching& dz =',&
914	dz(n+1), 'after a specific amount of', &
915	' grid levels& exceeds the upper', &
916	' limit =', stretch_factor_upper_limit, &
917	' &or lower limit = ', &
918	stretch_factor_lower_limit
919	CALL message( 'init_grid', 'PA0499', 0, 1, 0, 6, 0 )
920
921	ENDIF
922	ENDDO
923
924	END SUBROUTINE calculate_stretching_factor
925
926
927	! Description:
928	! -----------------------------------------------------------------------------!
929	!> Set temporary topography flags and reference buildings on top of underlying
930	!> orography.
931	!------------------------------------------------------------------------------!
932	SUBROUTINE process_topography( topo_3d )
933
934	USE arrays_3d, &
935	ONLY: zu, zw
936
937	USE control_parameters, &
938	ONLY: bc_lr_cyc, bc_ns_cyc, ocean_mode
939
940	USE exchange_horiz_mod, &
941	ONLY: exchange_horiz_int, exchange_horiz_2d
942
943	USE indices, &
944	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, &
945	nzt
946
947	USE netcdf_data_input_mod, &
948	ONLY: buildings_f, building_id_f, building_type_f, &
949	init_model, &
950	input_pids_static, &
951	terrain_height_f
952
953	USE kinds
954
955	USE pegrid
956
957	IMPLICIT NONE
958
959	INTEGER(iwp) :: i !< running index along x-direction
960	INTEGER(iwp) :: j !< running index along y-direction
961	INTEGER(iwp) :: k !< running index along z-direction with respect to numeric grid
962	INTEGER(iwp) :: k2 !< running index along z-direction with respect to netcdf grid
963	INTEGER(iwp) :: nr !< index variable indication maximum terrain height for respective building ID
964	INTEGER(iwp) :: num_build !< counter for number of buildings
965	INTEGER(iwp) :: topo_top_index !< orography top index, used to map 3D buildings onto terrain
966
967	#if defined( __parallel )
968	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displace_dum !< displacements of start addresses, used for MPI_ALLGATHERV
969	#endif
970	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids !< building IDs on entire model domain
971	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final !< building IDs on entire model domain, multiple occurences are sorted out
972	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final_tmp !< temporary array used for resizing
973	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l !< building IDs on local subdomain
974	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l_tmp !< temporary array used to resize array of building IDs
975
976	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings !< number of buildings with different ID on entire model domain
977	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings_l !< number of buildings with different ID on local subdomain
978
979	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags
980
981	REAL(wp) :: ocean_offset !< offset to consider inverse vertical coordinate at topography definition
982	REAL(wp) :: oro_min = 0.0_wp !< minimum terrain height in entire model domain, used to reference terrain to zero
983	REAL(wp), DIMENSION(:), ALLOCATABLE :: oro_max !< maximum terrain height occupied by an building with certain id
984	REAL(wp), DIMENSION(:), ALLOCATABLE :: oro_max_l !< maximum terrain height occupied by an building with certain id, on local subdomain
985
986	!
987	!-- Reference lowest terrain height to zero. This ensures that first,
988	!-- non-required gird levels (those which lie entirely below the minimum
989	!-- orography) are avoided, and second, that also negative orography can be used
990	!-- within the input file.
991	!-- Please note, in case of a nested run, the global minimum from all parent and
992	!-- childs need to be remove to avoid steep edges at the child-domain boundaries.
993	IF ( input_pids_static ) THEN
994
995	#if defined( __parallel )
996	CALL MPI_ALLREDUCE( MINVAL( terrain_height_f%var ), oro_min, 1, &
997	MPI_REAL, MPI_MIN, MPI_COMM_WORLD, ierr )
998	#else
999	oro_min = MINVAL( terrain_height_f%var )
1000	#endif
1001	terrain_height_f%var = terrain_height_f%var - oro_min
1002	!
1003	!-- Update reference height used within output files
1004	init_model%origin_z = init_model%origin_z + oro_min
1005
1006	ENDIF
1007
1008	!
1009	!-- In the following, buildings and orography are further preprocessed
1010	!-- before they are mapped on the LES grid.
1011	!-- Buildings are mapped on top of the orography by maintaining the roof
1012	!-- shape of the building. This can be achieved by referencing building on
1013	!-- top of the maximum terrain height within the area occupied by the
1014	!-- respective building. As buildings and terrain height are defined PE-wise,
1015	!-- parallelization of this referencing is required (a building can be
1016	!-- distributed between different PEs).
1017	!-- In a first step, determine the number of buildings with different
1018	!-- building id on each PE. In a next step, all building ids are gathered
1019	!-- into one array which is present to all PEs. For each building ID,
1020	!-- the maximum terrain height occupied by the respective building is
1021	!-- computed and distributed to each PE.
1022	!-- Finally, for each building id and its respective reference orography,
1023	!-- builidings are mapped on top.
1024	!--
1025	!-- First, pre-set topography flags, bit 1 indicates orography, bit 2
1026	!-- buildings
1027	!-- classify the respective surfaces.
1028	topo_3d = IBSET( topo_3d, 0 )
1029	topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
1030	!
1031	!-- In order to map topography on PALM grid also in case of ocean simulations,
1032	!-- pre-calculate an offset value.
1033	ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode )
1034	!
1035	!-- Reference buildings on top of orography. This is not necessary
1036	!-- if topography is read from ASCII file as no distinction between buildings
1037	!-- and terrain height can be made. Moreover, this is also not necessary if
1038	!-- urban-surface and land-surface model are used at the same time.
1039	IF ( input_pids_static ) THEN
1040
1041	IF ( buildings_f%from_file ) THEN
1042	num_buildings_l = 0
1043	num_buildings = 0
1044	!
1045	!-- Allocate at least one element for building ids and give it an inital
1046	!-- negative value that will be overwritten later. This, however, is
1047	!-- necessary in case there all IDs in the model domain are fill values.
1048	ALLOCATE( build_ids_l(1) )
1049	build_ids_l = -1
1050	DO i = nxl, nxr
1051	DO j = nys, nyn
1052	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1053	IF ( num_buildings_l(myid) > 0 ) THEN
1054	IF ( ANY( building_id_f%var(j,i) == build_ids_l ) ) &
1055	THEN
1056	CYCLE
1057	ELSE
1058	num_buildings_l(myid) = num_buildings_l(myid) + 1
1059	!
1060	!-- Resize array with different local building ids
1061	ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
1062	build_ids_l_tmp = build_ids_l
1063	DEALLOCATE( build_ids_l )
1064	ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
1065	build_ids_l(1:num_buildings_l(myid)-1) = &
1066	build_ids_l_tmp(1:num_buildings_l(myid)-1)
1067	build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
1068	DEALLOCATE( build_ids_l_tmp )
1069	ENDIF
1070	!
1071	!-- First occuring building id on PE
1072	ELSE
1073	num_buildings_l(myid) = num_buildings_l(myid) + 1
1074	build_ids_l(1) = building_id_f%var(j,i)
1075	ENDIF
1076	ENDIF
1077	ENDDO
1078	ENDDO
1079	!
1080	!-- Determine number of different building ids for the entire domain
1081	#if defined( __parallel )
1082	CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs, &
1083	MPI_INTEGER, MPI_SUM, comm2d, ierr )
1084	#else
1085	num_buildings = num_buildings_l
1086	#endif
1087	!
1088	!-- Gather all buildings ids on each PEs.
1089	!-- First, allocate array encompassing all building ids in model domain.
1090	ALLOCATE( build_ids(1:SUM(num_buildings)) )
1091	#if defined( __parallel )
1092	!
1093	!-- Allocate array for displacements.
1094	!-- As each PE may has a different number of buildings, so that
1095	!-- the block sizes send by each PE may not be equal. Hence,
1096	!-- information about the respective displacement is required, indicating
1097	!-- the respective adress where each MPI-task writes into the receive
1098	!-- buffer array
1099	ALLOCATE( displace_dum(0:numprocs-1) )
1100	displace_dum(0) = 0
1101	DO i = 1, numprocs-1
1102	displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
1103	ENDDO
1104
1105	CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)), &
1106	num_buildings(myid), &
1107	MPI_INTEGER, &
1108	build_ids, &
1109	num_buildings, &
1110	displace_dum, &
1111	MPI_INTEGER, &
1112	comm2d, ierr )
1113
1114	DEALLOCATE( displace_dum )
1115
1116	#else
1117	build_ids = build_ids_l
1118	#endif
1119
1120	!
1121	!-- Note, in parallel mode building ids can occure mutliple times, as
1122	!-- each PE has send its own ids. Therefore, sort out building ids which
1123	!-- appear more than one time.
1124	num_build = 0
1125	DO nr = 1, SIZE(build_ids)
1126
1127	IF ( ALLOCATED(build_ids_final) ) THEN
1128	IF ( ANY( build_ids(nr) == build_ids_final ) ) THEN
1129	CYCLE
1130	ELSE
1131	num_build = num_build + 1
1132	!
1133	!-- Resize
1134	ALLOCATE( build_ids_final_tmp(1:num_build) )
1135	build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
1136	DEALLOCATE( build_ids_final )
1137	ALLOCATE( build_ids_final(1:num_build) )
1138	build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
1139	build_ids_final(num_build) = build_ids(nr)
1140	DEALLOCATE( build_ids_final_tmp )
1141	ENDIF
1142	ELSE
1143	num_build = num_build + 1
1144	ALLOCATE( build_ids_final(1:num_build) )
1145	build_ids_final(num_build) = build_ids(nr)
1146	ENDIF
1147	ENDDO
1148
1149	!
1150	!-- Determine maximumum terrain height occupied by the respective
1151	!-- building and temporalily store on oro_max
1152	ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
1153	ALLOCATE( oro_max(1:SIZE(build_ids_final)) )
1154	oro_max_l = 0.0_wp
1155
1156	DO nr = 1, SIZE(build_ids_final)
1157	oro_max_l(nr) = MAXVAL( &
1158	MERGE( terrain_height_f%var(nys:nyn,nxl:nxr), &
1159	0.0_wp, &
1160	building_id_f%var(nys:nyn,nxl:nxr) == &
1161	build_ids_final(nr) ) )
1162	ENDDO
1163
1164	#if defined( __parallel )
1165	IF ( SIZE(build_ids_final) >= 1 ) THEN
1166	CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL,&
1167	MPI_MAX, comm2d, ierr )
1168	ENDIF
1169	#else
1170	oro_max = oro_max_l
1171	#endif
1172	!
1173	!-- Finally, determine discrete grid height of maximum orography occupied
1174	!-- by a building. Use all-or-nothing approach, i.e. if terrain
1175	!-- exceeds the scalar level the grid box is fully terrain and the
1176	!-- maximum terrain is set to the zw level.
1177	!-- terrain or
1178	oro_max_l = 0.0
1179	DO nr = 1, SIZE(build_ids_final)
1180	DO k = nzb, nzt
1181	IF ( zu(k) - ocean_offset <= oro_max(nr) ) &
1182	oro_max_l(nr) = zw(k) - ocean_offset
1183	ENDDO
1184	oro_max(nr) = oro_max_l(nr)
1185	ENDDO
1186	ENDIF
1187	!
1188	!-- Allocate array for storing terrain height under buildings
1189	IF ( buildings_f%from_file ) THEN
1190	ALLOCATE( buildings_f%oro_max(nysg:nyng,nxlg:nxrg) )
1191	buildings_f%oro_max = buildings_f%fill1
1192	END IF
1193	!
1194	!-- Map orography as well as buildings onto grid.
1195	DO i = nxl, nxr
1196	DO j = nys, nyn
1197	topo_top_index = 0
1198	!
1199	!-- Obtain index in global building_id array
1200	IF ( buildings_f%from_file ) THEN
1201	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1202	!
1203	!-- Determine index where maximum terrain height occupied by
1204	!-- the respective building height is stored.
1205	nr = MINLOC( ABS( build_ids_final - &
1206	building_id_f%var(j,i) ), DIM = 1 )
1207	!
1208	!-- Save grid-indexed oro_max
1209	buildings_f%oro_max(j,i) = oro_max(nr)
1210	ENDIF
1211	ENDIF
1212	DO k = nzb, nzt
1213	!
1214	!-- In a first step, if grid point is below or equal the given
1215	!-- terrain height, grid point is flagged to be of type natural.
1216	!-- Please note, in case there is also a building which is lower
1217	!-- than the vertical grid spacing, initialization of surface
1218	!-- attributes will not be correct as given surface information
1219	!-- will not be in accordance to the classified grid points.
1220	!-- Hence, in this case, also a building flag.
1221	IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) ) THEN
1222	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1223	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1224	topo_top_index = k ! topo_top_index + 1
1225	ENDIF
1226	!
1227	!-- Set building grid points. Here, only consider 2D buildings.
1228	!-- 3D buildings require separate treatment.
1229	IF ( buildings_f%from_file .AND. buildings_f%lod == 1 ) THEN
1230	!
1231	!-- Fill-up the terrain to the level of maximum orography
1232	!-- within the building-covered area.
1233	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1234	!
1235	!-- Note, oro_max is always on zw level
1236	IF ( zu(k) - ocean_offset < oro_max(nr) ) THEN
1237	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1238	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1239	ELSEIF ( zu(k) - ocean_offset <= &
1240	oro_max(nr) + buildings_f%var_2d(j,i) ) THEN
1241	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1242	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1243	ENDIF
1244	ENDIF
1245	ENDIF
1246	ENDDO
1247	!
1248	!-- Special treatment for non grid-resolved buildings. This case,
1249	!-- the uppermost terrain grid point is flagged as building as well
1250	!-- well, even though no building exists at all. However, the
1251	!-- surface element will be identified as urban-surface and the
1252	!-- input data provided by the drivers is consistent to the surface
1253	!-- classification. Else, all non grid-resolved buildings would vanish
1254	!-- and identified as terrain grid points, which, however, won't be
1255	!-- consistent with the input data.
1256	IF ( buildings_f%from_file .AND. buildings_f%lod == 1 ) THEN
1257	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1258	DO k = nzb, nzt
1259	IF( zw(k) - ocean_offset == oro_max(nr) ) THEN
1260	IF ( buildings_f%var_2d(j,i) <= zu(k+1) - zw(k) ) THEN
1261	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1262	ENDIF
1263	ENDIF
1264	ENDDO
1265	ENDIF
1266	ENDIF
1267	!
1268	!-- Map 3D buildings onto terrain height.
1269	!-- In case of any slopes, map building on top of maximum terrain
1270	!-- height covered by the building. In other words, extend
1271	!-- building down to the respective local terrain-surface height.
1272	IF ( buildings_f%from_file .AND. buildings_f%lod == 2 ) THEN
1273	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1274	!
1275	!-- Extend building down to the terrain surface, i.e. fill-up
1276	!-- surface irregularities below a building. Note, oro_max
1277	!-- is already a discrete height according to the all-or-nothing
1278	!-- approach, i.e. grid box is either topography or atmosphere,
1279	!-- terrain top is defined at upper bound of the grid box.
1280	!-- Hence, check for zw in this case.
1281	!-- Note, do this only for buildings which are surface mounted,
1282	!-- i.e. building types 1-6. Below bridges, which are represented
1283	!-- exclusively by building type 7, terrain shape should be
1284	!-- maintained.
1285	IF ( building_type_f%from_file ) THEN
1286	IF ( building_type_f%var(j,i) /= 7 ) THEN
1287	DO k = topo_top_index + 1, nzt + 1
1288	IF ( zu(k) - ocean_offset <= oro_max(nr) ) THEN
1289	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1290	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1291	ENDIF
1292	ENDDO
1293	!
1294	!-- After surface irregularities are smoothen, determine
1295	!-- lower start index where building starts.
1296	DO k = nzb, nzt
1297	IF ( zu(k) - ocean_offset <= oro_max(nr) ) &
1298	topo_top_index = k
1299	ENDDO
1300	ENDIF
1301	ENDIF
1302	!
1303	!-- Finally, map building on top.
1304	k2 = 0
1305	DO k = topo_top_index, nzt + 1
1306	IF ( k2 <= buildings_f%nz-1 ) THEN
1307	IF ( buildings_f%var_3d(k2,j,i) == 1 ) THEN
1308	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1309	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1310	ENDIF
1311	ENDIF
1312	k2 = k2 + 1
1313	ENDDO
1314	ENDIF
1315	ENDIF
1316	ENDDO
1317	ENDDO
1318	!
1319	!-- Horizontal exchange the oro_max array, which is required to for
1320	!-- initialization of building-surface properties.
1321	IF ( ALLOCATED( buildings_f%oro_max ) ) THEN
1322	CALL exchange_horiz_2d( buildings_f%oro_max(:,:) )
1323	ENDIF
1324	!
1325	!-- Deallocate temporary arrays required for processing and reading data
1326	IF ( ALLOCATED( oro_max ) ) DEALLOCATE( oro_max )
1327	IF ( ALLOCATED( oro_max_l ) ) DEALLOCATE( oro_max_l )
1328	IF ( ALLOCATED( build_ids_final ) ) DEALLOCATE( build_ids_final )
1329	!
1330	!-- Topography input via ASCII format.
1331	ELSE
1332	ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode )
1333	!
1334	!-- Initialize topography bit 0 (indicates obstacle) everywhere to zero
1335	!-- and clear all grid points at nzb, where alway a surface is defined.
1336	!-- Further, set also bit 1 (indicates terrain) at nzb, which is further
1337	!-- used for masked data output and further processing. Note, in the
1338	!-- ASCII case no distinction is made between buildings and terrain,
1339	!-- so that setting of bit 1 and 2 at the same time has no effect.
1340	topo_3d = IBSET( topo_3d, 0 )
1341	topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
1342	topo_3d(nzb,:,:) = IBSET( topo_3d(nzb,:,:), 1 )
1343	DO i = nxl, nxr
1344	DO j = nys, nyn
1345	DO k = nzb, nzt
1346	!
1347	!-- Flag topography for all grid points which are below
1348	!-- the local topography height.
1349	!-- Note, each topography is flagged as building.
1350	IF ( zu(k) - ocean_offset <= buildings_f%var_2d(j,i) ) THEN
1351	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1352	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates building
1353	ENDIF
1354	ENDDO
1355	ENDDO
1356	ENDDO
1357	ENDIF
1358
1359	CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1360
1361	IF ( .NOT. bc_ns_cyc ) THEN
1362	IF ( nys == 0 ) topo_3d(:,-1,:) = topo_3d(:,0,:)
1363	IF ( nyn == ny ) topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
1364	ENDIF
1365
1366	IF ( .NOT. bc_lr_cyc ) THEN
1367	IF ( nxl == 0 ) topo_3d(:,:,-1) = topo_3d(:,:,0)
1368	IF ( nxr == nx ) topo_3d(:,:,nx+1) = topo_3d(:,:,nx)
1369	ENDIF
1370
1371	END SUBROUTINE process_topography
1372
1373
1374	! Description:
1375	! -----------------------------------------------------------------------------!
1376	!> Filter topography, i.e. fill holes resolved by only one grid point.
1377	!> Such holes are suspected to lead to velocity blow-ups as continuity
1378	!> equation on discrete grid cannot be fulfilled in such case.
1379	!------------------------------------------------------------------------------!
1380	SUBROUTINE filter_topography( topo_3d )
1381
1382	USE control_parameters, &
1383	ONLY: bc_lr_cyc, bc_ns_cyc, message_string
1384
1385	USE exchange_horiz_mod, &
1386	ONLY: exchange_horiz_int, exchange_horiz_2d_byte, exchange_horiz_2d_int
1387
1388	USE indices, &
1389	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
1390
1391	USE netcdf_data_input_mod, &
1392	ONLY: building_id_f, building_type_f
1393
1394	USE pegrid
1395
1396	IMPLICIT NONE
1397
1398	LOGICAL :: filled = .FALSE. !< flag indicating if holes were filled
1399
1400	INTEGER(iwp) :: i !< running index along x-direction
1401	INTEGER(iwp) :: j !< running index along y-direction
1402	INTEGER(iwp) :: k !< running index along z-direction
1403	INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point
1404	INTEGER(iwp) :: num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE
1405	INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point
1406
1407	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo_tmp !< temporary 3D-topography used to fill holes
1408	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo_3d !< 3D-topography array merging buildings and orography
1409	!
1410	!-- Before checking for holes, set lateral boundary conditions for
1411	!-- topography. After hole-filling, boundary conditions must be set again.
1412	!-- Several iterations are performed, in order to fill holes which might
1413	!-- emerge by the filling-algorithm itself.
1414	ALLOCATE( topo_tmp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1415	topo_tmp = 0
1416
1417	num_hole = 99999
1418	DO WHILE ( num_hole > 0 )
1419
1420	num_hole = 0
1421	CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1422	!
1423	!-- Exchange also building ID and type. Note, building_type is an one-byte
1424	!-- variable.
1425	IF ( building_id_f%from_file ) &
1426	CALL exchange_horiz_2d_int( building_id_f%var, nys, nyn, nxl, nxr, nbgp )
1427	IF ( building_type_f%from_file ) &
1428	CALL exchange_horiz_2d_byte( building_type_f%var, nys, nyn, nxl, nxr, nbgp )
1429
1430	topo_tmp = topo_3d
1431	!
1432	!-- In case of non-cyclic lateral boundaries, assume lateral boundary to be
1433	!-- a solid wall. Thus, intermediate spaces of one grid point between
1434	!-- boundary and some topographic structure will be filled.
1435	IF ( .NOT. bc_ns_cyc ) THEN
1436	IF ( nys == 0 ) topo_tmp(:,-1,:) = IBCLR( topo_tmp(:,0,:), 0 )
1437	IF ( nyn == ny ) topo_tmp(:,ny+1,:) = IBCLR( topo_tmp(:,ny,:), 0 )
1438	ENDIF
1439
1440	IF ( .NOT. bc_lr_cyc ) THEN
1441	IF ( nxl == 0 ) topo_tmp(:,:,-1) = IBCLR( topo_tmp(:,:,0), 0 )
1442	IF ( nxr == nx ) topo_tmp(:,:,nx+1) = IBCLR( topo_tmp(:,:,nx), 0 )
1443	ENDIF
1444
1445	num_hole_l = 0
1446	DO i = nxl, nxr
1447	DO j = nys, nyn
1448	DO k = nzb+1, nzt
1449	IF ( BTEST( topo_tmp(k,j,i), 0 ) ) THEN
1450	num_wall = 0
1451	IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) ) &
1452	num_wall = num_wall + 1
1453	IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) ) &
1454	num_wall = num_wall + 1
1455	IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) ) &
1456	num_wall = num_wall + 1
1457	IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) ) &
1458	num_wall = num_wall + 1
1459	IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) ) &
1460	num_wall = num_wall + 1
1461	IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) ) &
1462	num_wall = num_wall + 1
1463
1464	IF ( num_wall >= 4 ) THEN
1465	num_hole_l = num_hole_l + 1
1466	!
1467	!-- Clear flag 0 and set special flag ( bit 4) to indicate
1468	!-- that new topography point is a result of filtering process.
1469	topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
1470	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 4 )
1471	!
1472	!-- If filled grid point is occupied by a building, classify
1473	!-- it as building grid point.
1474	IF ( building_type_f%from_file ) THEN
1475	IF ( building_type_f%var(j,i) /= &
1476	building_type_f%fill .OR. &
1477	building_type_f%var(j+1,i) /= &
1478	building_type_f%fill .OR. &
1479	building_type_f%var(j-1,i) /= &
1480	building_type_f%fill .OR. &
1481	building_type_f%var(j,i+1) /= &
1482	building_type_f%fill .OR. &
1483	building_type_f%var(j,i-1) /= &
1484	building_type_f%fill ) THEN
1485	!
1486	!-- Set flag indicating building surfaces
1487	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
1488	!
1489	!-- Set building_type and ID at this position if not
1490	!-- already set. This is required for proper
1491	!-- initialization of urban-surface energy balance
1492	!-- solver.
1493	IF ( building_type_f%var(j,i) == &
1494	building_type_f%fill ) THEN
1495
1496	IF ( building_type_f%var(j+1,i) /= &
1497	building_type_f%fill ) THEN
1498	building_type_f%var(j,i) = &
1499	building_type_f%var(j+1,i)
1500	building_id_f%var(j,i) = &
1501	building_id_f%var(j+1,i)
1502	ELSEIF ( building_type_f%var(j-1,i) /= &
1503	building_type_f%fill ) THEN
1504	building_type_f%var(j,i) = &
1505	building_type_f%var(j-1,i)
1506	building_id_f%var(j,i) = &
1507	building_id_f%var(j-1,i)
1508	ELSEIF ( building_type_f%var(j,i+1) /= &
1509	building_type_f%fill ) THEN
1510	building_type_f%var(j,i) = &
1511	building_type_f%var(j,i+1)
1512	building_id_f%var(j,i) = &
1513	building_id_f%var(j,i+1)
1514	ELSEIF ( building_type_f%var(j,i-1) /= &
1515	building_type_f%fill ) THEN
1516	building_type_f%var(j,i) = &
1517	building_type_f%var(j,i-1)
1518	building_id_f%var(j,i) = &
1519	building_id_f%var(j,i-1)
1520	ENDIF
1521	ENDIF
1522	ENDIF
1523	ENDIF
1524	!
1525	!-- If filled grid point is already classified as building
1526	!-- everything is fine, else classify this grid point as
1527	!-- natural type grid point. This case, values for the
1528	!-- surface type are already set.
1529	IF ( .NOT. BTEST( topo_3d(k,j,i), 2 ) ) THEN
1530	topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
1531	ENDIF
1532	ENDIF
1533	ENDIF
1534	ENDDO
1535	ENDDO
1536	ENDDO
1537	!
1538	!-- Count the total number of holes, required for informative message.
1539	#if defined( __parallel )
1540	CALL MPI_ALLREDUCE( num_hole_l, num_hole, 1, MPI_INTEGER, MPI_SUM, &
1541	comm2d, ierr )
1542	#else
1543	num_hole = num_hole_l
1544	#endif
1545	IF ( num_hole > 0 .AND. .NOT. filled ) filled = .TRUE.
1546
1547	ENDDO
1548	!
1549	!-- Create an informative message if any holes were filled.
1550	IF ( filled ) THEN
1551	WRITE( message_string, * ) 'Topography was filtered, i.e. holes ' // &
1552	'resolved by only one grid point ' // &
1553	'were filled during initialization.'
1554	CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 )
1555	ENDIF
1556
1557	DEALLOCATE( topo_tmp )
1558	!
1559	!-- Finally, exchange topo_3d array again and if necessary set Neumann boundary
1560	!-- condition in case of non-cyclic lateral boundaries.
1561	CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
1562
1563	IF ( .NOT. bc_ns_cyc ) THEN
1564	IF ( nys == 0 ) topo_3d(:,-1,:) = topo_3d(:,0,:)
1565	IF ( nyn == ny ) topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
1566	ENDIF
1567
1568	IF ( .NOT. bc_lr_cyc ) THEN
1569	IF ( nxl == 0 ) topo_3d(:,:,-1) = topo_3d(:,:,0)
1570	IF ( nxr == nx ) topo_3d(:,:,nx+1) = topo_3d(:,:,nx)
1571	ENDIF
1572	!
1573	!-- Exchange building ID and type. Note, building_type is an one-byte variable.
1574	IF ( building_id_f%from_file ) &
1575	CALL exchange_horiz_2d_int( building_id_f%var, nys, nyn, nxl, nxr, nbgp )
1576	IF ( building_type_f%from_file ) &
1577	CALL exchange_horiz_2d_byte( building_type_f%var, nys, nyn, nxl, nxr, nbgp )
1578
1579	END SUBROUTINE filter_topography
1580
1581
1582	! Description:
1583	! -----------------------------------------------------------------------------!
1584	!> Reads topography information from file or sets generic topography. Moreover,
1585	!> all topography-relevant topography arrays are initialized, and grid flags
1586	!> are set.
1587	!------------------------------------------------------------------------------!
1588	SUBROUTINE init_topo( topo )
1589
1590	USE arrays_3d, &
1591	ONLY: zw
1592
1593	USE control_parameters, &
1594	ONLY: bc_lr_cyc, bc_ns_cyc, building_height, building_length_x, &
1595	building_length_y, building_wall_left, building_wall_south, &
1596	canyon_height, canyon_wall_left, canyon_wall_south, &
1597	canyon_width_x, canyon_width_y, dp_level_ind_b, dz, &
1598	message_string, topography, topography_grid_convention, &
1599	tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y, &
1600	tunnel_wall_depth
1601
1602	USE exchange_horiz_mod, &
1603	ONLY: exchange_horiz_int
1604
1605	USE grid_variables, &
1606	ONLY: dx, dy
1607
1608	USE indices, &
1609	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, &
1610	nzb, nzt
1611
1612	USE kinds
1613
1614	USE netcdf_data_input_mod, &
1615	ONLY: buildings_f, terrain_height_f
1616
1617	USE pegrid
1618
1619	IMPLICIT NONE
1620
1621	INTEGER(iwp) :: bh !< temporary vertical index of building height
1622	INTEGER(iwp) :: ngp_bx !< grid point number of building size along x
1623	INTEGER(iwp) :: ngp_by !< grid point number of building size along y
1624	INTEGER(iwp) :: index_left_bwall !< index for left building wall
1625	INTEGER(iwp) :: index_right_bwall !< index for right building wall
1626	INTEGER(iwp) :: index_north_bwall !< index for north building wall
1627	INTEGER(iwp) :: index_south_bwall !< index for south building wall
1628	INTEGER(iwp) :: ch !< temporary vertical index for canyon height
1629	INTEGER(iwp) :: ngp_cx !< grid point number of canyon size along x
1630	INTEGER(iwp) :: ngp_cy !< grid point number of canyon size along y
1631	INTEGER(iwp) :: index_left_cwall !< index for left canyon wall
1632	INTEGER(iwp) :: index_right_cwall !< index for right canyon wall
1633	INTEGER(iwp) :: index_north_cwall !< index for north canyon wall
1634	INTEGER(iwp) :: index_south_cwall !< index for south canyon wall
1635	INTEGER(iwp) :: i !< index variable along x
1636	INTEGER(iwp) :: j !< index variable along y
1637	INTEGER(iwp) :: k !< index variable along z
1638	INTEGER(iwp) :: hv_in !< heavyside function to model inner tunnel surface
1639	INTEGER(iwp) :: hv_out !< heavyside function to model outer tunnel surface
1640	INTEGER(iwp) :: txe_out !< end position of outer tunnel wall in x
1641	INTEGER(iwp) :: txs_out !< start position of outer tunnel wall in x
1642	INTEGER(iwp) :: tye_out !< end position of outer tunnel wall in y
1643	INTEGER(iwp) :: tys_out !< start position of outer tunnel wall in y
1644	INTEGER(iwp) :: txe_in !< end position of inner tunnel wall in x
1645	INTEGER(iwp) :: txs_in !< start position of inner tunnel wall in x
1646	INTEGER(iwp) :: tye_in !< end position of inner tunnel wall in y
1647	INTEGER(iwp) :: tys_in !< start position of inner tunnel wall in y
1648	INTEGER(iwp) :: td !< tunnel wall depth
1649	INTEGER(iwp) :: th !< height of outer tunnel wall
1650
1651	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center
1652	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo !< input array for 3D topography and dummy array for setting "outer"-flags
1653	!
1654	!-- Check for correct setting of the namelist parameter topography. If
1655	!-- topography information is read from file but topography = 'flat',
1656	!-- initialization does not work properly.
1657	IF ( ( buildings_f%from_file .OR. terrain_height_f%from_file ) .AND. &
1658	TRIM( topography ) /= 'read_from_file' ) THEN
1659	message_string = 'If topography information is provided (via ' // &
1660	'Netcdf or ASCII input) topography = ' // &
1661	'"read_from_file" is required.'
1662	CALL message( 'init_grid', 'PA0437', 1, 2, 0, 6, 0 )
1663	ENDIF
1664	!
1665	!-- Set outer and inner index arrays for non-flat topography.
1666	!-- Here consistency checks concerning domain size and periodicity are
1667	!-- necessary.
1668	!-- Within this SELECT CASE structure only nzb_local is initialized
1669	!-- individually depending on the chosen topography type, all other index
1670	!-- arrays are initialized further below.
1671	SELECT CASE ( TRIM( topography ) )
1672
1673	CASE ( 'flat' )
1674	!
1675	!-- Initialilize 3D topography array, used later for initializing flags
1676	topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
1677
1678	CASE ( 'closed_channel' )
1679	!
1680	!-- Initialilize 3D topography array, used later for initializing flags
1681	topo(nzb+1:nzt,:,:) = IBSET( topo(nzb+1:nzt,:,:), 0 )
1682
1683	CASE ( 'single_building' )
1684	!
1685	!-- Single rectangular building, by default centered in the middle of the
1686	!-- total domain
1687	ngp_bx = NINT( building_length_x / dx )
1688	ngp_by = NINT( building_length_y / dy )
1689	bh = MINLOC( ABS( zw - building_height ), 1 ) - 1
1690	IF ( ABS( zw(bh) - building_height ) == &
1691	ABS( zw(bh+1) - building_height ) ) bh = bh + 1
1692	IF ( building_wall_left == 9999999.9_wp ) THEN
1693	building_wall_left = ( nx + 1 - ngp_bx ) / 2 * dx
1694	ENDIF
1695	index_left_bwall = NINT( building_wall_left / dx )
1696	index_right_bwall = index_left_bwall + ngp_bx
1697
1698	IF ( building_wall_south == 9999999.9_wp ) THEN
1699	building_wall_south = ( ny + 1 - ngp_by ) / 2 * dy
1700	ENDIF
1701	index_south_bwall = NINT( building_wall_south / dy )
1702	index_north_bwall = index_south_bwall + ngp_by
1703
1704	!
1705	!-- Building size has to meet some requirements
1706	IF ( ( index_left_bwall < 1 ) .OR. ( index_right_bwall > nx-1 ) .OR. &
1707	( index_right_bwall < index_left_bwall+3 ) .OR. &
1708	( index_south_bwall < 1 ) .OR. ( index_north_bwall > ny-1 ) .OR.&
1709	( index_north_bwall < index_south_bwall+3 ) ) THEN
1710	WRITE( message_string, * ) 'inconsistent building parameters:', &
1711	'&index_left_bwall=', index_left_bwall, &
1712	'index_right_bwall=', index_right_bwall, &
1713	'index_south_bwall=', index_south_bwall, &
1714	'index_north_bwall=', index_north_bwall, &
1715	'nx=', nx, 'ny=', ny
1716	CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
1717	ENDIF
1718
1719	ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
1720	nzb_local = 0
1721	!
1722	!-- Define the building.
1723	IF ( index_left_bwall <= nxr .AND. index_right_bwall >= nxl .AND. &
1724	index_south_bwall <= nyn .AND. index_north_bwall >= nys ) &
1725	nzb_local(MAX(nys,index_south_bwall):MIN(nyn,index_north_bwall), &
1726	MAX(nxl,index_left_bwall):MIN(nxr,index_right_bwall)) = bh
1727	!
1728	!-- Set bit array on basis of nzb_local
1729	DO i = nxl, nxr
1730	DO j = nys, nyn
1731	topo(nzb_local(j,i)+1:nzt+1,j,i) = &
1732	IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 )
1733	ENDDO
1734	ENDDO
1735
1736	DEALLOCATE( nzb_local )
1737
1738	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1739	!
1740	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
1741	!-- boundary conditions for topography.
1742	IF ( .NOT. bc_ns_cyc ) THEN
1743	IF ( nys == 0 ) THEN
1744	DO i = 1, nbgp
1745	topo(:,nys-i,:) = topo(:,nys,:)
1746	ENDDO
1747	ENDIF
1748	IF ( nyn == ny ) THEN
1749	DO i = 1, nbgp
1750	topo(:,nyn+i,:) = topo(:,nyn,:)
1751	ENDDO
1752	ENDIF
1753	ENDIF
1754	IF ( .NOT. bc_lr_cyc ) THEN
1755	IF ( nxl == 0 ) THEN
1756	DO i = 1, nbgp
1757	topo(:,:,nxl-i) = topo(:,:,nxl)
1758	ENDDO
1759	ENDIF
1760	IF ( nxr == nx ) THEN
1761	DO i = 1, nbgp
1762	topo(:,:,nxr+i) = topo(:,:,nxr)
1763	ENDDO
1764	ENDIF
1765	ENDIF
1766
1767	CASE ( 'single_street_canyon' )
1768	!
1769	!-- Single quasi-2D street canyon of infinite length in x or y direction.
1770	!-- The canyon is centered in the other direction by default.
1771	IF ( canyon_width_x /= 9999999.9_wp ) THEN
1772	!
1773	!-- Street canyon in y direction
1774	ngp_cx = NINT( canyon_width_x / dx )
1775	IF ( canyon_wall_left == 9999999.9_wp ) THEN
1776	canyon_wall_left = ( nx + 1 - ngp_cx ) / 2 * dx
1777	ENDIF
1778	index_left_cwall= NINT( canyon_wall_left / dx )
1779	index_right_cwall= index_left_cwall+ ngp_cx
1780	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
1781	!
1782	!-- Street canyon in x direction
1783	ngp_cy = NINT( canyon_width_y / dy )
1784	IF ( canyon_wall_south == 9999999.9_wp ) THEN
1785	canyon_wall_south = ( ny + 1 - ngp_cy ) / 2 * dy
1786	ENDIF
1787	index_south_cwall = NINT( canyon_wall_south / dy )
1788	index_north_cwall = index_south_cwall + ngp_cy
1789
1790	ELSE
1791
1792	message_string = 'no street canyon width given'
1793	CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
1794
1795	ENDIF
1796
1797	ch = MINLOC( ABS( zw - canyon_height ), 1 ) - 1
1798	IF ( ABS( zw(ch) - canyon_height ) == &
1799	ABS( zw(ch+1) - canyon_height ) ) ch = ch + 1
1800	dp_level_ind_b = ch
1801	!
1802	!-- Street canyon size has to meet some requirements
1803	IF ( canyon_width_x /= 9999999.9_wp ) THEN
1804	IF ( ( index_left_cwall< 1 ) .OR. ( index_right_cwall> nx-1 ) .OR.&
1805	( ngp_cx < 3 ) ) THEN
1806	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
1807	'&index_left_cwall=', index_left_cwall, &
1808	' index_right_cwall=', index_right_cwall, &
1809	' ngp_cx=', ngp_cx, &
1810	' ch=', ch, ' nx=', nx, ' ny=', ny
1811	CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
1812	ENDIF
1813	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
1814	IF ( ( index_south_cwall < 1 ) .OR. &
1815	( index_north_cwall > ny-1 ) .OR. ( ngp_cy < 3 ) ) THEN
1816	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
1817	'&index_south_cwall=', index_south_cwall, &
1818	' index_north_cwall=', index_north_cwall, &
1819	' ngp_cy=', ngp_cy, &
1820	' ch=', ch, ' nx=', nx, ' ny=', ny
1821	CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
1822	ENDIF
1823	ENDIF
1824	IF ( canyon_width_x /= 9999999.9_wp .AND. &
1825	canyon_width_y /= 9999999.9_wp ) THEN
1826	message_string = 'inconsistent canyon parameters:' // &
1827	'&street canyon can only be oriented' // &
1828	' either in x- or in y-direction'
1829	CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
1830	ENDIF
1831
1832	ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
1833	nzb_local = ch
1834	IF ( canyon_width_x /= 9999999.9_wp ) THEN
1835	IF ( index_left_cwall<= nxr .AND. index_right_cwall>= nxl ) &
1836	nzb_local(:,MAX(nxl,index_left_cwall+1):MIN(nxr,index_right_cwall-1)) = 0
1837	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
1838	IF ( index_south_cwall <= nyn .AND. index_north_cwall >= nys ) &
1839	nzb_local(MAX(nys,index_south_cwall+1):MIN(nyn,index_north_cwall-1),:) = 0
1840	ENDIF
1841	!
1842	!-- Set bit array on basis of nzb_local
1843	DO i = nxl, nxr
1844	DO j = nys, nyn
1845	topo(nzb_local(j,i)+1:nzt+1,j,i) = &
1846	IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 )
1847	ENDDO
1848	ENDDO
1849	DEALLOCATE( nzb_local )
1850
1851	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
1852	!
1853	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
1854	!-- boundary conditions for topography.
1855	IF ( .NOT. bc_ns_cyc ) THEN
1856	IF ( nys == 0 ) THEN
1857	DO i = 1, nbgp
1858	topo(:,nys-i,:) = topo(:,nys,:)
1859	ENDDO
1860	ENDIF
1861	IF ( nyn == ny ) THEN
1862	DO i = 1, nbgp
1863	topo(:,nyn+i,:) = topo(:,nyn,:)
1864	ENDDO
1865	ENDIF
1866	ENDIF
1867	IF ( .NOT. bc_lr_cyc ) THEN
1868	IF ( nxl == 0 ) THEN
1869	DO i = 1, nbgp
1870	topo(:,:,nxl-i) = topo(:,:,nxl)
1871	ENDDO
1872	ENDIF
1873	IF ( nxr == nx ) THEN
1874	DO i = 1, nbgp
1875	topo(:,:,nxr+i) = topo(:,:,nxr)
1876	ENDDO
1877	ENDIF
1878	ENDIF
1879
1880	CASE ( 'tunnel' )
1881
1882	!
1883	!-- Tunnel height
1884	IF ( tunnel_height == 9999999.9_wp ) THEN
1885	th = zw( INT( 0.2 * nz) )
1886	ELSE
1887	th = tunnel_height
1888	ENDIF
1889	!
1890	!-- Tunnel-wall depth
1891	IF ( tunnel_wall_depth == 9999999.9_wp ) THEN
1892	td = MAX ( dx, dy, dz(1) )
1893	ELSE
1894	td = tunnel_wall_depth
1895	ENDIF
1896	!
1897	!-- Check for tunnel width
1898	IF ( tunnel_width_x == 9999999.9_wp .AND. &
1899	tunnel_width_y == 9999999.9_wp ) THEN
1900	message_string = 'No tunnel width is given. '
1901	CALL message( 'init_grid', 'PA0280', 1, 2, 0, 6, 0 )
1902	ENDIF
1903	IF ( tunnel_width_x /= 9999999.9_wp .AND. &
1904	tunnel_width_y /= 9999999.9_wp ) THEN
1905	message_string = 'Inconsistent tunnel parameters:' // &
1906	'tunnel can only be oriented' // &
1907	'either in x- or in y-direction.'
1908	CALL message( 'init_grid', 'PA0281', 1, 2, 0, 6, 0 )
1909	ENDIF
1910	!
1911	!-- Check for too small tunnel width in x- and y-direction
1912	IF ( tunnel_width_x /= 9999999.9_wp .AND. &
1913	tunnel_width_x - 2.0_wp * td <= 2.0_wp * dx ) THEN
1914	message_string = 'tunnel_width_x too small'
1915	CALL message( 'init_grid', 'PA0175', 1, 2, 0, 6, 0 )
1916	ENDIF
1917	IF ( tunnel_width_y /= 9999999.9_wp .AND. &
1918	tunnel_width_y - 2.0_wp * td <= 2.0_wp * dy ) THEN
1919	message_string = 'tunnel_width_y too small'
1920	CALL message( 'init_grid', 'PA0455', 1, 2, 0, 6, 0 )
1921	ENDIF
1922	!
1923	!-- Check for too large tunnel width.
1924	!-- Tunnel axis along y.
1925	IF ( tunnel_width_x /= 9999999.9_wp ) THEN
1926	IF ( tunnel_width_x > ( nx + 1 ) * dx ) THEN
1927	message_string = 'tunnel_width_x too large'
1928	CALL message( 'init_grid', 'PA0282', 1, 2, 0, 6, 0 )
1929	ENDIF
1930
1931	txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_width_x * 0.5_wp )
1932	txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_width_x * 0.5_wp )
1933	txs_in = INT( ( nx + 1 ) * 0.5_wp * dx - &
1934	( tunnel_width_x * 0.5_wp - td ) )
1935	txe_in = INT( ( nx + 1 ) * 0.5_wp * dx + &
1936	( tunnel_width_x * 0.5_wp - td ) )
1937
1938	tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_length * 0.5_wp )
1939	tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_length * 0.5_wp )
1940	tys_in = tys_out
1941	tye_in = tye_out
1942	ENDIF
1943	!
1944	!-- Tunnel axis along x.
1945	IF ( tunnel_width_y /= 9999999.9_wp ) THEN
1946	IF ( tunnel_width_y > ( ny + 1 ) * dy ) THEN
1947	message_string = 'tunnel_width_y too large'
1948	CALL message( 'init_grid', 'PA0456', 1, 2, 0, 6, 0 )
1949	ENDIF
1950
1951	txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_length * 0.5_wp )
1952	txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_length * 0.5_wp )
1953	txs_in = txs_out
1954	txe_in = txe_out
1955
1956	tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_width_y * 0.5_wp )
1957	tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_width_y * 0.5_wp )
1958	tys_in = INT( ( ny + 1 ) * 0.5_wp * dy - &
1959	( tunnel_width_y * 0.5_wp - td ) )
1960	tye_in = INT( ( ny + 1 ) * 0.5_wp * dy + &
1961	( tunnel_width_y * 0.5_wp - td ) )
1962	ENDIF
1963
1964	topo = 0
1965	DO i = nxl, nxr
1966	DO j = nys, nyn
1967	!
1968	!-- Use heaviside function to model outer tunnel surface
1969	hv_out = th * 0.5_wp * &
1970	( ( SIGN( 1.0_wp, i * dx - txs_out ) + 1.0_wp ) &
1971	- ( SIGN( 1.0_wp, i * dx - txe_out ) + 1.0_wp ) )
1972
1973	hv_out = hv_out * 0.5_wp * &
1974	( ( SIGN( 1.0_wp, j * dy - tys_out ) + 1.0_wp ) &
1975	- ( SIGN( 1.0_wp, j * dy - tye_out ) + 1.0_wp ) )
1976	!
1977	!-- Use heaviside function to model inner tunnel surface
1978	hv_in = ( th - td ) * 0.5_wp * &
1979	( ( SIGN( 1.0_wp, i * dx - txs_in ) + 1.0_wp ) &
1980	- ( SIGN( 1.0_wp, i * dx - txe_in ) + 1.0_wp ) )
1981
1982	hv_in = hv_in * 0.5_wp * &
1983	( ( SIGN( 1.0_wp, j * dy - tys_in ) + 1.0_wp ) &
1984	- ( SIGN( 1.0_wp, j * dy - tye_in ) + 1.0_wp ) )
1985	!
1986	!-- Set flags at x-y-positions without any tunnel surface
1987	IF ( hv_out - hv_in == 0.0_wp ) THEN
1988	topo(nzb+1:nzt+1,j,i) = IBSET( topo(nzb+1:nzt+1,j,i), 0 )
1989	!
1990	!-- Set flags at x-y-positions with tunnel surfaces
1991	ELSE
1992	DO k = nzb + 1, nzt + 1
1993	!
1994	!-- Inner tunnel
1995	IF ( hv_out - hv_in == th ) THEN
1996	IF ( zw(k) <= hv_out ) THEN
1997	topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
1998	ELSE
1999	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2000	ENDIF
2001	ENDIF
2002	!
2003	!-- Lateral tunnel walls
2004	IF ( hv_out - hv_in == td ) THEN
2005	IF ( zw(k) <= hv_in ) THEN
2006	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2007	ELSEIF ( zw(k) > hv_in .AND. zw(k) <= hv_out ) THEN
2008	topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
2009	ELSEIF ( zw(k) > hv_out ) THEN
2010	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2011	ENDIF
2012	ENDIF
2013	ENDDO
2014	ENDIF
2015	ENDDO
2016	ENDDO
2017
2018	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2019	!
2020	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
2021	!-- boundary conditions for topography.
2022	IF ( .NOT. bc_ns_cyc ) THEN
2023	IF ( nys == 0 ) THEN
2024	DO i = 1, nbgp
2025	topo(:,nys-i,:) = topo(:,nys,:)
2026	ENDDO
2027	ENDIF
2028	IF ( nyn == ny ) THEN
2029	DO i = 1, nbgp
2030	topo(:,nyn+i,:) = topo(:,nyn,:)
2031	ENDDO
2032	ENDIF
2033	ENDIF
2034	IF ( .NOT. bc_lr_cyc ) THEN
2035	IF ( nxl == 0 ) THEN
2036	DO i = 1, nbgp
2037	topo(:,:,nxl-i) = topo(:,:,nxl)
2038	ENDDO
2039	ENDIF
2040	IF ( nxr == nx ) THEN
2041	DO i = 1, nbgp
2042	topo(:,:,nxr+i) = topo(:,:,nxr)
2043	ENDDO
2044	ENDIF
2045	ENDIF
2046
2047	CASE ( 'read_from_file' )
2048	!
2049	!-- Note, topography information have been already read.
2050	!-- If required, further process topography, i.e. reference buildings on
2051	!-- top of orography and set temporary 3D topography array, which is
2052	!-- used later to set grid flags. Calling of this rouinte is also
2053	!-- required in case of ASCII input, even though no distinction between
2054	!-- terrain- and building height is made in this case.
2055	CALL process_topography( topo )
2056	!
2057	!-- Filter holes resolved by only one grid-point
2058	CALL filter_topography( topo )
2059	!
2060	!-- Exchange ghost-points, as well as add cyclic or Neumann boundary
2061	!-- conditions.
2062	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2063	!
2064	!-- Set lateral boundary conditions for topography on all ghost layers
2065	IF ( .NOT. bc_ns_cyc ) THEN
2066	IF ( nys == 0 ) THEN
2067	DO i = 1, nbgp
2068	topo(:,nys-i,:) = topo(:,nys,:)
2069	ENDDO
2070	ENDIF
2071	IF ( nyn == ny ) THEN
2072	DO i = 1, nbgp
2073	topo(:,nyn+i,:) = topo(:,nyn,:)
2074	ENDDO
2075	ENDIF
2076	ENDIF
2077
2078	IF ( .NOT. bc_lr_cyc ) THEN
2079	IF ( nxl == 0 ) THEN
2080	DO i = 1, nbgp
2081	topo(:,:,nxl-i) = topo(:,:,nxl)
2082	ENDDO
2083	ENDIF
2084	IF ( nxr == nx ) THEN
2085	DO i = 1, nbgp
2086	topo(:,:,nxr+i) = topo(:,:,nxr)
2087	ENDDO
2088	ENDIF
2089	ENDIF
2090
2091
2092	CASE DEFAULT
2093	!
2094	!-- The DEFAULT case is reached either if the parameter topography
2095	!-- contains a wrong character string or if the user has defined a special
2096	!-- case in the user interface. There, the subroutine user_init_grid
2097	!-- checks which of these two conditions applies.
2098	CALL user_init_grid( topo )
2099	CALL filter_topography( topo )
2100
2101	END SELECT
2102	!
2103	!-- Consistency checks and index array initialization are only required for
2104	!-- non-flat topography.
2105	IF ( TRIM( topography ) /= 'flat' ) THEN
2106	!
2107	!-- In case of non-flat topography, check whether the convention how to
2108	!-- define the topography grid has been set correctly, or whether the default
2109	!-- is applicable. If this is not possible, abort.
2110	IF ( TRIM( topography_grid_convention ) == ' ' ) THEN
2111	IF ( TRIM( topography ) /= 'closed_channel' .AND. &
2112	TRIM( topography ) /= 'single_building' .AND. &
2113	TRIM( topography ) /= 'single_street_canyon' .AND. &
2114	TRIM( topography ) /= 'tunnel' .AND. &
2115	TRIM( topography ) /= 'read_from_file') THEN
2116	!-- The default value is not applicable here, because it is only valid
2117	!-- for the four standard cases 'single_building',
2118	!-- 'single_street_canyon', 'tunnel' and 'read_from_file'
2119	!-- defined in init_grid.
2120	WRITE( message_string, * ) &
2121	'The value for "topography_grid_convention" ', &
2122	'is not set. Its default value is & only valid for ', &
2123	'"topography" = ''single_building'', ''tunnel'' ', &
2124	'''single_street_canyon'', ''closed_channel'' & or ', &
2125	'''read_from_file''.', &
2126	'& Choose ''cell_edge'' or ''cell_center''.'
2127	CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 )
2128	ELSE
2129	!-- The default value is applicable here.
2130	!-- Set convention according to topography.
2131	IF ( TRIM( topography ) == 'single_building' .OR. &
2132	TRIM( topography ) == 'single_street_canyon' ) THEN
2133	topography_grid_convention = 'cell_edge'
2134	ELSEIF ( TRIM( topography ) == 'read_from_file' .OR. &
2135	TRIM( topography ) == 'tunnel') THEN
2136	topography_grid_convention = 'cell_center'
2137	ENDIF
2138	ENDIF
2139	ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND. &
2140	TRIM( topography_grid_convention ) /= 'cell_center' ) THEN
2141	WRITE( message_string, * ) &
2142	'The value for "topography_grid_convention" is ', &
2143	'not recognized.& Choose ''cell_edge'' or ''cell_center''.'
2144	CALL message( 'init_grid', 'PA0240', 1, 2, 0, 6, 0 )
2145	ENDIF
2146
2147
2148	IF ( topography_grid_convention == 'cell_edge' ) THEN
2149	!
2150	!-- The array nzb_local as defined using the 'cell_edge' convention
2151	!-- describes the actual total size of topography which is defined at the
2152	!-- cell edges where u=0 on the topography walls in x-direction and v=0
2153	!-- on the topography walls in y-direction. However, PALM uses individual
2154	!-- arrays nzb_u\|v\|w\|s_inner\|outer that are based on nzb_s_inner.
2155	!-- Therefore, the extent of topography in nzb_local is now reduced by
2156	!-- 1dx at the E topography walls and by 1dy at the N topography walls
2157	!-- to form the basis for nzb_s_inner.
2158	!-- Note, the reverse memory access (i-j instead of j-i) is absolutely
2159	!-- required at this point.
2160	DO j = nys+1, nyn+1
2161	DO i = nxl-1, nxr
2162	DO k = nzb, nzt+1
2163	IF ( BTEST( topo(k,j,i), 0 ) .OR. &
2164	BTEST( topo(k,j,i+1), 0 ) ) &
2165	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2166	ENDDO
2167	ENDDO
2168	ENDDO
2169	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2170
2171	DO i = nxl, nxr+1
2172	DO j = nys-1, nyn
2173	DO k = nzb, nzt+1
2174	IF ( BTEST( topo(k,j,i), 0 ) .OR. &
2175	BTEST( topo(k,j+1,i), 0 ) ) &
2176	topo(k,j,i) = IBSET( topo(k,j,i), 0 )
2177	ENDDO
2178	ENDDO
2179	ENDDO
2180	CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
2181
2182	ENDIF
2183	ENDIF
2184
2185
2186	END SUBROUTINE init_topo
2187
2188	SUBROUTINE set_topo_flags(topo)
2189
2190	USE control_parameters, &
2191	ONLY: bc_lr_cyc, bc_ns_cyc, constant_flux_layer, &
2192	scalar_advec, topography, use_surface_fluxes, use_top_fluxes
2193
2194	USE exchange_horiz_mod, &
2195	ONLY: exchange_horiz_int
2196
2197	USE indices, &
2198	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, &
2199	nzt, topo_top_ind, wall_flags_static_0, wall_flags_total_0
2200
2201	USE kinds
2202
2203	IMPLICIT NONE
2204
2205	INTEGER(iwp) :: i !< index variable along x
2206	INTEGER(iwp) :: ibit !< integer bit position of topgraphy masking array
2207	INTEGER(iwp) :: j !< index variable along y
2208	INTEGER(iwp) :: k !< index variable along z
2209
2210	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo !< input array for 3D topography and dummy array for setting "outer"-flags
2211
2212	ALLOCATE( wall_flags_static_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2213	wall_flags_static_0 = 0
2214	!
2215	!-- Set-up topography flags. First, set flags only for s, u, v and w-grid.
2216	!-- Further special flags will be set in following loops.
2217	DO i = nxl, nxr
2218	DO j = nys, nyn
2219	DO k = nzb, nzt+1
2220	!
2221	!-- scalar grid
2222	IF ( BTEST( topo(k,j,i), 0 ) ) &
2223	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 0 )
2224	!
2225	!-- u grid
2226	IF ( BTEST( topo(k,j,i), 0 ) .AND. &
2227	BTEST( topo(k,j,i-1), 0 ) ) &
2228	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 1 )
2229	!
2230	!-- v grid
2231	IF ( BTEST( topo(k,j,i), 0 ) .AND. &
2232	BTEST( topo(k,j-1,i), 0 ) ) &
2233	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 2 )
2234
2235	ENDDO
2236
2237	DO k = nzb, nzt
2238	!
2239	!-- w grid
2240	IF ( BTEST( topo(k,j,i), 0 ) .AND. &
2241	BTEST( topo(k+1,j,i), 0 ) ) &
2242	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 3 )
2243	ENDDO
2244
2245	IF ( topography /= 'closed_channel' ) THEN
2246	wall_flags_static_0(nzt+1,j,i) = IBSET( wall_flags_static_0(nzt+1,j,i), 3 )
2247	ENDIF
2248
2249	ENDDO
2250	ENDDO
2251
2252	CALL exchange_horiz_int( wall_flags_static_0, nys, nyn, nxl, nxr, nzt, nbgp )
2253
2254	!
2255	!-- Set outer array for scalars to mask near-surface grid points. Note, on
2256	!-- basis of flag 24 futher flags will be derived which are used to control
2257	!-- production of subgrid TKE production near walls.
2258
2259	ALLOCATE( wall_flags_total_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2260	wall_flags_total_0 = 0
2261
2262	DO i = nxl, nxr
2263	DO j = nys, nyn
2264	DO k = nzb, nzt+1
2265	wall_flags_total_0(k,j,i) = IOR( wall_flags_total_0(k,j,i), wall_flags_static_0(k,j,i) )
2266	ENDDO
2267	ENDDO
2268	ENDDO
2269
2270	CALL exchange_horiz_int( wall_flags_total_0, nys, nyn, nxl, nxr, nzt, nbgp )
2271
2272	DO i = nxl, nxr
2273	DO j = nys, nyn
2274	DO k = nzb, nzt+1
2275	IF ( BTEST( wall_flags_total_0(k,j-1,i), 0 ) .AND. &
2276	BTEST( wall_flags_total_0(k,j+1,i), 0 ) .AND. &
2277	BTEST( wall_flags_total_0(k,j,i-1), 0 ) .AND. &
2278	BTEST( wall_flags_total_0(k,j,i+1), 0 ) .AND. &
2279	BTEST( wall_flags_total_0(k,j-1,i-1), 0 ) .AND. &
2280	BTEST( wall_flags_total_0(k,j+1,i-1), 0 ) .AND. &
2281	BTEST( wall_flags_total_0(k,j-1,i+1), 0 ) .AND. &
2282	BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) ) &
2283	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 24 )
2284	ENDDO
2285	ENDDO
2286	ENDDO
2287	!
2288	!-- Set further special flags
2289	DO i = nxl, nxr
2290	DO j = nys, nyn
2291	DO k = nzb, nzt+1
2292	!
2293	!-- scalar grid, former nzb_diff_s_inner.
2294	!-- Note, use this flag also to mask topography in diffusion_u and
2295	!-- diffusion_v along the vertical direction. In case of
2296	!-- use_surface_fluxes, fluxes are calculated via MOST, else, simple
2297	!-- gradient approach is applied. Please note, in case of u- and v-
2298	!-- diffuison, a small error is made at edges (on the east side for u,
2299	!-- at the north side for v), since topography on scalar grid point
2300	!-- is used instead of topography on u/v-grid. As number of topography grid
2301	!-- points on uv-grid is different than s-grid, different number of
2302	!-- surface elements would be required. In order to avoid this,
2303	!-- treat edges (u(k,j,i+1)) simply by a gradient approach, i.e. these
2304	!-- points are not masked within diffusion_u. Tests had shown that the
2305	!-- effect on the flow is negligible.
2306	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2307	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) ) &
2308	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 8 )
2309	ELSE
2310	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 8 )
2311	ENDIF
2312
2313	ENDDO
2314	!
2315	!-- Special flag to control vertical diffusion at model top - former
2316	!-- nzt_diff
2317	wall_flags_total_0(:,j,i) = IBSET( wall_flags_total_0(:,j,i), 9 )
2318	IF ( use_top_fluxes ) &
2319	wall_flags_total_0(nzt+1,j,i) = IBCLR( wall_flags_total_0(nzt+1,j,i), 9 )
2320
2321
2322	DO k = nzb+1, nzt
2323	!
2324	!-- Special flag on u grid, former nzb_u_inner + 1, required
2325	!-- for disturb_field and initialization. Do not disturb directly at
2326	!-- topography, as well as initialize u with zero one grid point outside
2327	!-- of topography.
2328	IF ( BTEST( wall_flags_total_0(k-1,j,i), 1 ) .AND. &
2329	BTEST( wall_flags_total_0(k,j,i), 1 ) .AND. &
2330	BTEST( wall_flags_total_0(k+1,j,i), 1 ) ) &
2331	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 20 )
2332	!
2333	!-- Special flag on v grid, former nzb_v_inner + 1, required
2334	!-- for disturb_field and initialization. Do not disturb directly at
2335	!-- topography, as well as initialize v with zero one grid point outside
2336	!-- of topography.
2337	IF ( BTEST( wall_flags_total_0(k-1,j,i), 2 ) .AND. &
2338	BTEST( wall_flags_total_0(k,j,i), 2 ) .AND. &
2339	BTEST( wall_flags_total_0(k+1,j,i), 2 ) ) &
2340	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 21 )
2341	!
2342	!-- Special flag on scalar grid, former nzb_s_inner+1. Used for
2343	!-- lpm_sgs_tke
2344	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
2345	BTEST( wall_flags_total_0(k-1,j,i), 0 ) .AND. &
2346	BTEST( wall_flags_total_0(k+1,j,i), 0 ) ) &
2347	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 25 )
2348	!
2349	!-- Special flag on scalar grid, nzb_diff_s_outer - 1, required in
2350	!-- in production_e
2351	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2352	IF ( BTEST( wall_flags_total_0(k,j,i), 24 ) .AND. &
2353	BTEST( wall_flags_total_0(k-1,j,i), 24 ) .AND. &
2354	BTEST( wall_flags_total_0(k+1,j,i), 0 ) ) &
2355	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 29 )
2356	ELSE
2357	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) ) &
2358	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 29 )
2359	ENDIF
2360	!
2361	!-- Special flag on scalar grid, nzb_diff_s_outer - 1, required in
2362	!-- in production_e
2363	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2364	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
2365	BTEST( wall_flags_total_0(k-1,j,i), 0 ) .AND. &
2366	BTEST( wall_flags_total_0(k+1,j,i), 0 ) ) &
2367	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 30 )
2368	ELSE
2369	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) ) &
2370	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 30 )
2371	ENDIF
2372	ENDDO
2373	!
2374	!-- Flags indicating downward facing walls
2375	DO k = nzb+1, nzt+1
2376	!
2377	!-- Scalar grid
2378	IF ( BTEST( wall_flags_total_0(k-1,j,i), 0 ) .AND. &
2379	.NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) ) &
2380	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 13 )
2381	!
2382	!-- Downward facing wall on u grid
2383	IF ( BTEST( wall_flags_total_0(k-1,j,i), 1 ) .AND. &
2384	.NOT. BTEST( wall_flags_total_0(k,j,i), 1 ) ) &
2385	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 15 )
2386	!
2387	!-- Downward facing wall on v grid
2388	IF ( BTEST( wall_flags_total_0(k-1,j,i), 2 ) .AND. &
2389	.NOT. BTEST( wall_flags_total_0(k,j,i), 2 ) ) &
2390	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 17 )
2391	!
2392	!-- Downward facing wall on w grid
2393	IF ( BTEST( wall_flags_total_0(k-1,j,i), 3 ) .AND. &
2394	.NOT. BTEST( wall_flags_total_0(k,j,i), 3 ) ) &
2395	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 19 )
2396	ENDDO
2397	!
2398	!-- Flags indicating upward facing walls
2399	DO k = nzb, nzt
2400	!
2401	!-- Upward facing wall on scalar grid
2402	IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
2403	BTEST( wall_flags_total_0(k+1,j,i), 0 ) ) &
2404	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 12 )
2405	!
2406	!-- Upward facing wall on u grid
2407	IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 1 ) .AND. &
2408	BTEST( wall_flags_total_0(k+1,j,i), 1 ) ) &
2409	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 14 )
2410
2411	!
2412	!-- Upward facing wall on v grid
2413	IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 2 ) .AND. &
2414	BTEST( wall_flags_total_0(k+1,j,i), 2 ) ) &
2415	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 16 )
2416
2417	!
2418	!-- Upward facing wall on w grid
2419	IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 3 ) .AND. &
2420	BTEST( wall_flags_total_0(k+1,j,i), 3 ) ) &
2421	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 18 )
2422	!
2423	!-- Special flag on scalar grid, former nzb_s_inner
2424	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) .OR. &
2425	BTEST( wall_flags_total_0(k,j,i), 12 ) .OR. &
2426	BTEST( wall_flags_total_0(k,j,i), 13 ) ) &
2427	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 22 )
2428	!
2429	!-- Special flag on scalar grid, nzb_diff_s_inner - 1, required for
2430	!-- flow_statistics
2431	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
2432	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
2433	BTEST( wall_flags_total_0(k+1,j,i), 0 ) ) &
2434	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 23 )
2435	ELSE
2436	IF ( BTEST( wall_flags_total_0(k,j,i), 22 ) ) &
2437	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 23 )
2438	ENDIF
2439
2440
2441	ENDDO
2442	wall_flags_total_0(nzt+1,j,i) = IBSET( wall_flags_total_0(nzt+1,j,i), 22 )
2443	wall_flags_total_0(nzt+1,j,i) = IBSET( wall_flags_total_0(nzt+1,j,i), 23 )
2444	!
2445	!-- Set flags indicating that topography is close by in horizontal
2446	!-- direction, i.e. flags that infold the topography. These will be used
2447	!-- to set advection flags for passive scalars, where due to large
2448	!-- gradients near buildings stationary numerical oscillations can produce
2449	!-- unrealistically high concentrations. This is only necessary if
2450	!-- WS-scheme is applied for scalar advection. Note, these flags will be
2451	!-- only used for passive scalars such as chemical species or aerosols.
2452	IF ( scalar_advec == 'ws-scheme' ) THEN
2453	DO k = nzb, nzt
2454	IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. ( &
2455	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-1), 0 ) ) .OR.&
2456	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-2), 0 ) ) .OR.&
2457	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-3), 0 ) ) .OR.&
2458	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i+1), 0 ) ) .OR.&
2459	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i+2), 0 ) ) .OR.&
2460	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i+3), 0 ) ) .OR.&
2461	ANY( .NOT. BTEST( wall_flags_total_0(k,j-1,i-3:i+3), 0 ) ) .OR.&
2462	ANY( .NOT. BTEST( wall_flags_total_0(k,j-2,i-3:i+3), 0 ) ) .OR.&
2463	ANY( .NOT. BTEST( wall_flags_total_0(k,j-3,i-3:i+3), 0 ) ) .OR.&
2464	ANY( .NOT. BTEST( wall_flags_total_0(k,j+1,i-3:i+3), 0 ) ) .OR.&
2465	ANY( .NOT. BTEST( wall_flags_total_0(k,j+2,i-3:i+3), 0 ) ) .OR.&
2466	ANY( .NOT. BTEST( wall_flags_total_0(k,j+3,i-3:i+3), 0 ) ) &
2467	) ) &
2468	wall_flags_total_0(k,j,i) = IBSET( wall_flags_total_0(k,j,i), 31 )
2469
2470	ENDDO
2471	ENDIF
2472	ENDDO
2473	ENDDO
2474	!
2475	!-- Finally, set identification flags indicating natural terrain or buildings.
2476	!-- Natural terrain grid points. Information on the type of the surface is
2477	!-- stored in bit 1 of 3D Integer array topo. However, this bit is only set
2478	!-- when topography is read from file. In order to run the land-surface model
2479	!-- also without topography information, set bit 1 explicitely in this case.
2480	!--
2481	!-- Natural terrain grid points
2482	!-- If no topography is initialized, the land-surface is at k = nzb.
2483	IF ( TRIM( topography ) /= 'read_from_file' ) THEN
2484	wall_flags_static_0(nzb,:,:) = IBSET( wall_flags_static_0(nzb,:,:), 5 )
2485	ELSE
2486	DO i = nxl, nxr
2487	DO j = nys, nyn
2488	DO k = nzb, nzt+1
2489	!
2490	!-- Natural terrain grid point
2491	IF ( BTEST( topo(k,j,i), 1 ) ) &
2492	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 5 )
2493	ENDDO
2494	ENDDO
2495	ENDDO
2496	ENDIF
2497	!
2498	!-- Building grid points.
2499	DO i = nxl, nxr
2500	DO j = nys, nyn
2501	DO k = nzb, nzt+1
2502	IF ( BTEST( topo(k,j,i), 2 ) ) &
2503	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 6 )
2504	ENDDO
2505	ENDDO
2506	ENDDO
2507	!
2508	!-- Set flag 4, indicating new topography grid points due to filtering.
2509	DO i = nxl, nxr
2510	DO j = nys, nyn
2511	DO k = nzb, nzt+1
2512	IF ( BTEST( topo(k,j,i), 4 ) ) &
2513	wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 4 )
2514	ENDDO
2515	ENDDO
2516	ENDDO
2517
2518	CALL exchange_horiz_int( wall_flags_static_0, nys, nyn, nxl, nxr, nzt, nbgp )
2519
2520	DO i = nxl, nxr
2521	DO j = nys, nyn
2522	DO k = nzb, nzt+1
2523	wall_flags_total_0(k,j,i) = IOR( wall_flags_total_0(k,j,i), wall_flags_static_0(k,j,i) )
2524	ENDDO
2525	ENDDO
2526	ENDDO
2527	!
2528	!-- Exchange ghost points for wall flags
2529	CALL exchange_horiz_int( wall_flags_total_0, nys, nyn, nxl, nxr, nzt, nbgp )
2530	!
2531	!-- Set boundary conditions also for flags. Can be interpreted as Neumann
2532	!-- boundary conditions for topography.
2533	IF ( .NOT. bc_ns_cyc ) THEN
2534	IF ( nys == 0 ) THEN
2535	DO i = 1, nbgp
2536	wall_flags_total_0(:,nys-i,:) = wall_flags_total_0(:,nys,:)
2537	ENDDO
2538	ENDIF
2539	IF ( nyn == ny ) THEN
2540	DO i = 1, nbgp
2541	wall_flags_total_0(:,nyn+i,:) = wall_flags_total_0(:,nyn,:)
2542	ENDDO
2543	ENDIF
2544	ENDIF
2545	IF ( .NOT. bc_lr_cyc ) THEN
2546	IF ( nxl == 0 ) THEN
2547	DO i = 1, nbgp
2548	wall_flags_total_0(:,:,nxl-i) = wall_flags_total_0(:,:,nxl)
2549	ENDDO
2550	ENDIF
2551	IF ( nxr == nx ) THEN
2552	DO i = 1, nbgp
2553	wall_flags_total_0(:,:,nxr+i) = wall_flags_total_0(:,:,nxr)
2554	ENDDO
2555	ENDIF
2556	ENDIF
2557	!
2558	!-- Pre-calculate topography top indices (former get_topography_top_index
2559	!-- function)
2560	ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:4) )
2561	!
2562	!-- Uppermost topography index on scalar grid
2563	ibit = 12
2564	topo_top_ind(:,:,0) = MAXLOC( &
2565	MERGE( 1, 0, &
2566	BTEST( wall_flags_total_0(:,:,:), ibit ) &
2567	), DIM = 1 &
2568	) - 1
2569	!
2570	!-- Uppermost topography index on u grid
2571	ibit = 14
2572	topo_top_ind(:,:,1) = MAXLOC( &
2573	MERGE( 1, 0, &
2574	BTEST( wall_flags_total_0(:,:,:), ibit ) &
2575	), DIM = 1 &
2576	) - 1
2577	!
2578	!-- Uppermost topography index on v grid
2579	ibit = 16
2580	topo_top_ind(:,:,2) = MAXLOC( &
2581	MERGE( 1, 0, &
2582	BTEST( wall_flags_total_0(:,:,:), ibit ) &
2583	), DIM = 1 &
2584	) - 1
2585	!
2586	!-- Uppermost topography index on w grid
2587	ibit = 18
2588	topo_top_ind(:,:,3) = MAXLOC( &
2589	MERGE( 1, 0, &
2590	BTEST( wall_flags_total_0(:,:,:), ibit ) &
2591	), DIM = 1 &
2592	) - 1
2593	!
2594	!-- Uppermost topography index on scalar outer grid
2595	ibit = 24
2596	topo_top_ind(:,:,4) = MAXLOC( &
2597	MERGE( 1, 0, &
2598	BTEST( wall_flags_total_0(:,:,:), ibit ) &
2599	), DIM = 1 &
2600	) - 1
2601
2602	END SUBROUTINE set_topo_flags

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |