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

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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