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

Last change on this file since 2893 was 2893, checked in by suehring, 3 years ago

Revise informative message concerning filtered topography.

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