Ignore:
Timestamp:
Dec 14, 2017 5:12:51 PM (6 years ago)
Author:
kanani
Message:

Merge of branch palm4u into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/init_grid.f90

    r2550 r2696  
    11!> @file init_grid.f90
    22!------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
     3! This file is part of the PALM model system.
    44!
    55! PALM is free software: you can redistribute it and/or modify it under the
     
    2525! -----------------
    2626! $Id$
     27! Revised topography input
     28! Set nzb_max not for the entire nest domain, only for boundary PEs
     29! Re-organize routine, split-up into several subroutines
     30! Modularize poismg_noopt
     31! Remove setting bit 26, 27, 28 in wall_flags_0, indicating former '_outer'
     32! arrays (not required any more). 
     33! Bugfix in generic tunnel setup (MS)
     34!
     35! 2550 2017-10-16 17:12:01Z boeske
    2736! Set lateral boundary conditions for topography on all three ghost layers
    2837!
     
    239248!
    240249! Description:
    241 ! ------------
     250! -----------------------------------------------------------------------------!
    242251!> Creating grid depending constants
    243 !> To Do: Setting topo flags only based on topo_3d array - flags for former
    244 !> nzb_outer arrays are still not set properly.
    245 !> To Do: Rearrange topo flag list
     252!> @todo: Move initialization mixing length
     253!> @todo: Rearrange topo flag list
     254!> @todo: reference 3D buildings on top of orography is not tested and may need
     255!>        further improvement for steep slopes
     256!> @todo: Use more advanced setting of building type at filled holes
    246257!------------------------------------------------------------------------------!
    247258 SUBROUTINE init_grid
     
    251262
    252263    USE arrays_3d,                                                             &
    253         ONLY:  dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzu_mg, dzw, dzw_mg, f1_mg,  &
    254                f2_mg, f3_mg, l_grid, l_wall, zu, zw
     264        ONLY:  dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, zu, zw
    255265       
    256266    USE control_parameters,                                                    &
     
    259269               canyon_height, canyon_wall_left, canyon_wall_south,             &
    260270               canyon_width_x, canyon_width_y, constant_flux_layer,            &
    261                coupling_char, coupling_mode,                                   &
    262271               dp_level_ind_b, dz, dz_max, dz_stretch_factor,                  &
    263                dz_stretch_level, dz_stretch_level_index, grid_level, ibc_uv_b, &
    264                io_blocks, io_group, inflow_l, inflow_n, inflow_r, inflow_s,    &
    265                lod, masking_method, maximum_grid_level, message_string,        &
     272               dz_stretch_level, dz_stretch_level_index, grid_level,           &
     273               force_bound_l, force_bound_r, force_bound_n, force_bound_s,     &
     274               ibc_uv_b, inflow_l, inflow_n, inflow_r, inflow_s,               &
     275               masking_method, maximum_grid_level, message_string,             &
    266276               momentum_advec, nest_domain, nest_bound_l, nest_bound_n,        &
    267277               nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n,        &
     
    275285       
    276286    USE indices,                                                               &
    277         ONLY:  advc_flags_1, advc_flags_2, flags, nbgp, nx, nxl, nxlg, nxl_mg, &
    278                nxr, nxrg, nxr_mg, ny, nyn, nyng, nyn_mg, nys, nys_mg, nysg, nz,&
     287        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
    279288               nzb, nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer,              &
    280289               nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,                 &
    281290               nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner,             &
    282                nzb_w_outer, nzt, nzt_mg, wall_flags_0, wall_flags_1,           &
    283                wall_flags_10, wall_flags_2, wall_flags_3,  wall_flags_4,       &
    284                wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,         &
    285                wall_flags_9
     291               nzb_w_outer, nzt
    286292   
    287293    USE kinds
    288 #if defined ( __netcdf )
    289     USE netcdf_interface,                                                      &
    290         ONLY:  netcdf_close_file, netcdf_open_read_file, netcdf_get_attribute, &
    291                netcdf_get_variable
    292 #endif
     294
    293295    USE pegrid
     296
     297    USE poismg_noopt_mod
    294298
    295299    USE surface_mod,                                                           &
     
    298302    USE vertical_nesting_mod,                                                  &
    299303        ONLY:  vnested, vnest_init_grid
     304
     305    IMPLICIT NONE
     306
     307    INTEGER(iwp) ::  i             !< index variable along x
     308    INTEGER(iwp) ::  j             !< index variable along y
     309    INTEGER(iwp) ::  k             !< index variable along z
     310    INTEGER(iwp) ::  k_top         !< topography top index on local PE
     311    INTEGER(iwp) ::  l             !< loop variable
     312    INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
     313    INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
     314                                     
     315    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
     316    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp        !< dummy to calculate topography indices on u- and v-grid
     317
     318    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
     319
     320    REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
     321
     322
     323!
     324!-- Calculation of horizontal array bounds including ghost layers
     325    nxlg = nxl - nbgp
     326    nxrg = nxr + nbgp
     327    nysg = nys - nbgp
     328    nyng = nyn + nbgp
     329
     330!
     331!-- Allocate grid arrays
     332    ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1),        &
     333              dzw(1:nzt+1), zu(nzb:nzt+1), zw(nzb:nzt+1) )
     334
     335!
     336!-- Compute height of u-levels from constant grid length and dz stretch factors
     337    IF ( dz == -1.0_wp )  THEN
     338       message_string = 'missing dz'
     339       CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
     340    ELSEIF ( dz <= 0.0_wp )  THEN
     341       WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
     342       CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
     343    ENDIF
     344
     345!
     346!-- Define the vertical grid levels
     347    IF ( .NOT. ocean )  THEN
     348!
     349!--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
     350!--    The second u-level (k=1) corresponds to the top of the
     351!--    Prandtl-layer.
     352
     353       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
     354          zu(0) = 0.0_wp
     355      !    zu(0) = - dz * 0.5_wp
     356       ELSE
     357          zu(0) = - dz * 0.5_wp
     358       ENDIF
     359       zu(1) =   dz * 0.5_wp
     360
     361       dz_stretch_level_index = nzt+1
     362       dz_stretched = dz
     363       DO  k = 2, nzt+1
     364          IF ( dz_stretch_level <= zu(k-1)  .AND.  dz_stretched < dz_max )  THEN
     365             dz_stretched = dz_stretched * dz_stretch_factor
     366             dz_stretched = MIN( dz_stretched, dz_max )
     367             IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
     368          ENDIF
     369          zu(k) = zu(k-1) + dz_stretched
     370       ENDDO
     371
     372!
     373!--    Compute the w-levels. They are always staggered half-way between the
     374!--    corresponding u-levels. In case of dirichlet bc for u and v at the
     375!--    ground the first u- and w-level (k=0) are defined at same height (z=0).
     376!--    The top w-level is extrapolated linearly.
     377       zw(0) = 0.0_wp
     378       DO  k = 1, nzt
     379          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
     380       ENDDO
     381       zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
     382
     383    ELSE
     384!
     385!--    Grid for ocean with free water surface is at k=nzt (w-grid).
     386!--    In case of neumann bc at the ground the first first u-level (k=0) lies
     387!--    below the first w-level (k=0). In case of dirichlet bc the first u- and
     388!--    w-level are defined at same height, but staggered from the second level.
     389!--    The second u-level (k=1) corresponds to the top of the Prandtl-layer.
     390       zu(nzt+1) =   dz * 0.5_wp
     391       zu(nzt)   = - dz * 0.5_wp
     392
     393       dz_stretch_level_index = 0
     394       dz_stretched = dz
     395       DO  k = nzt-1, 0, -1
     396!
     397!--       The default value of dz_stretch_level is positive, thus the first
     398!--       condition is always true. Hence, the second condition is necessary.
     399          IF ( dz_stretch_level >= zu(k+1)  .AND.  dz_stretch_level <= 0.0  &
     400               .AND.  dz_stretched < dz_max )  THEN
     401             dz_stretched = dz_stretched * dz_stretch_factor
     402             dz_stretched = MIN( dz_stretched, dz_max )
     403             IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
     404          ENDIF
     405          zu(k) = zu(k+1) - dz_stretched
     406       ENDDO
     407
     408!
     409!--    Compute the w-levels. They are always staggered half-way between the
     410!--    corresponding u-levels, except in case of dirichlet bc for u and v
     411!--    at the ground. In this case the first u- and w-level are defined at
     412!--    same height. The top w-level (nzt+1) is not used but set for
     413!--    consistency, since w and all scalar variables are defined up tp nzt+1.
     414       zw(nzt+1) = dz
     415       zw(nzt)   = 0.0_wp
     416       DO  k = 0, nzt
     417          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
     418       ENDDO
     419
     420!
     421!--    In case of dirichlet bc for u and v the first u- and w-level are defined
     422!--    at same height.
     423       IF ( ibc_uv_b == 0 ) THEN
     424          zu(0) = zw(0)
     425       ENDIF
     426
     427    ENDIF
     428
     429!
     430!-- Compute grid lengths.
     431    DO  k = 1, nzt+1
     432       dzu(k)  = zu(k) - zu(k-1)
     433       ddzu(k) = 1.0_wp / dzu(k)
     434       dzw(k)  = zw(k) - zw(k-1)
     435       ddzw(k) = 1.0_wp / dzw(k)
     436    ENDDO
     437
     438    DO  k = 1, nzt
     439       dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
     440    ENDDO
     441   
     442!   
     443!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
     444!-- everywhere. For the actual grid, the grid spacing at the lowest level
     445!-- is only dz/2, but should be dz. Therefore, an additional array
     446!-- containing with appropriate grid information is created for these
     447!-- solvers.
     448    IF ( psolver(1:9) /= 'multigrid' )  THEN
     449       ALLOCATE( ddzu_pres(1:nzt+1) )
     450       ddzu_pres = ddzu
     451       ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
     452    ENDIF
     453
     454!
     455!-- Compute the reciprocal values of the horizontal grid lengths.
     456    ddx = 1.0_wp / dx
     457    ddy = 1.0_wp / dy
     458    dx2 = dx * dx
     459    dy2 = dy * dy
     460    ddx2 = 1.0_wp / dx2
     461    ddy2 = 1.0_wp / dy2
     462
     463!
     464!-- Allocate 3D array to set topography
     465    ALLOCATE( topo(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     466    topo = 0
     467!
     468!-- Initialize topography by generic topography or read from topography from file. 
     469    CALL init_topo( topo )
     470!
     471!-- Set flags to mask topography on the grid.
     472    CALL set_topo_flags( topo )   
     473!
     474!-- Calculate wall flag arrays for the multigrid method.
     475!-- Please note, wall flags are only applied in the non-optimized version.
     476    IF ( psolver == 'multigrid_noopt' )  CALL poismg_noopt_init
     477
     478!
     479!-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e.
     480!-- to decrease the numerical stencil appropriately.
     481    IF ( momentum_advec == 'ws-scheme'  .OR.  scalar_advec == 'ws-scheme' )    &
     482       CALL ws_init_flags
     483
     484!
     485!-- Calculated grid-dependent as well as near-wall mixing length
     486    CALL init_mixing_length
     487
     488!
     489!-- Determine the maximum level of topography. It is used for
     490!-- steering the degradation of order of the applied advection scheme,
     491!-- as well in the lpm.
     492!-- In case of non-cyclic lateral boundaries, the order of the advection
     493!-- scheme has to be reduced up to nzt (required at the lateral boundaries).
     494    k_top = 0
     495    DO  i = nxl, nxr
     496       DO  j = nys, nyn
     497          DO  k = nzb, nzt + 1
     498             k_top = MAX( k_top, MERGE( k, 0,                                  &
     499                                        .NOT. BTEST( topo(k,j,i), 0 ) ) )
     500          ENDDO
     501       ENDDO
     502    ENDDO
     503#if defined( __parallel )
     504    CALL MPI_ALLREDUCE( k_top + 1, nzb_max, 1, MPI_INTEGER,                    & !is +1 really necessary here?
     505                        MPI_MAX, comm2d, ierr )
     506#else
     507    nzb_max = k_top + 1
     508#endif
     509    IF ( inflow_l  .OR.  outflow_l  .OR.  force_bound_l  .OR.  nest_bound_l  .OR.&
     510         inflow_r  .OR.  outflow_r  .OR.  force_bound_r  .OR.  nest_bound_r  .OR.&
     511         inflow_n  .OR.  outflow_n  .OR.  force_bound_n  .OR.  nest_bound_n  .OR.&
     512         inflow_s  .OR.  outflow_s  .OR.  force_bound_s  .OR.  nest_bound_s )    &
     513         nzb_max = nzt
     514!   
     515!-- Finally, if topography extents up to the model top, limit nzb_max to nzt.
     516    nzb_max = MIN( nzb_max, nzt )
     517
     518!
     519!-- Initialize boundary conditions via surface type
     520    CALL init_bc
     521
     522!
     523!-- Allocate and set topography height arrays required for data output
     524    IF ( TRIM( topography ) /= 'flat' )  THEN
     525!
     526!--    Allocate and set the arrays containing the topography height
     527       IF ( nxr == nx  .AND.  nyn /= ny )  THEN
     528          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
     529                    zw_w_inner(nxl:nxr+1,nys:nyn) )
     530       ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
     531          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
     532                    zw_w_inner(nxl:nxr,nys:nyn+1) )
     533       ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
     534          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
     535                    zw_w_inner(nxl:nxr+1,nys:nyn+1) )
     536       ELSE
     537          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
     538                    zw_w_inner(nxl:nxr,nys:nyn) )
     539       ENDIF
     540
     541       zu_s_inner   = 0.0_wp
     542       zw_w_inner   = 0.0_wp
     543!
     544!--    Determine local topography height on scalar and w-grid. Note, setting
     545!--    lateral boundary values is not necessary, realized via wall_flags_0
     546!--    array. Further, please note that loop bounds are different from
     547!--    nxl to nxr and nys to nyn on south and right model boundary, hence,
     548!--    use intrinsic lbound and ubound functions to infer array bounds.
     549       DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
     550          DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
     551!
     552!--          Topography height on scalar grid. Therefore, determine index of
     553!--          upward-facing surface element on scalar grid.
     554             zu_s_inner(i,j) = zu( get_topography_top_index( j, i, 's' ) )
     555!
     556!--          Topography height on w grid. Therefore, determine index of
     557!--          upward-facing surface element on w grid.
     558             zw_w_inner(i,j) = zw( get_topography_top_index( j, i, 's' ) )
     559          ENDDO
     560       ENDDO
     561    ENDIF
     562
     563!
     564!-- In the following, calculate 2D index arrays. Note, these will be removed
     565!-- soon.
     566!-- Allocate outer and inner index arrays for topography and set
     567!-- defaults.                   
     568    ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg),                                &
     569              nzb_s_outer(nysg:nyng,nxlg:nxrg),                                &
     570              nzb_u_inner(nysg:nyng,nxlg:nxrg),                                &
     571              nzb_u_outer(nysg:nyng,nxlg:nxrg),                                &
     572              nzb_v_inner(nysg:nyng,nxlg:nxrg),                                &
     573              nzb_v_outer(nysg:nyng,nxlg:nxrg),                                &
     574              nzb_w_inner(nysg:nyng,nxlg:nxrg),                                &
     575              nzb_w_outer(nysg:nyng,nxlg:nxrg),                                &
     576              nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                           &
     577              nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                           &
     578              nzb_local(nysg:nyng,nxlg:nxrg),                                  &
     579              nzb_tmp(nysg:nyng,nxlg:nxrg) )
     580!
     581!-- Initialize 2D-index arrays. Note, these will be removed soon!
     582    nzb_local(nys:nyn,nxl:nxr) = get_topography_top_index( 's' )
     583    CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
     584
     585    IF ( TRIM( topography ) /= 'flat' )  THEN
     586#if defined( __parallel )
     587       CALL MPI_ALLREDUCE( MAXVAL( get_topography_top_index( 's' ) ),          &
     588                           nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
     589       CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ),          &
     590                           nzb_local_min, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )                   
     591#else
     592       nzb_local_max = MAXVAL( get_topography_top_index( 's' ) )
     593       nzb_local_min = MINVAL( get_topography_top_index( 's' ) )
     594#endif
     595!
     596!--    Consistency checks
     597       IF ( nzb_local_min < 0  .OR.  nzb_local_max  > nz + 1 )  THEN
     598          WRITE( message_string, * ) 'nzb_local values are outside the',       &
     599                                'model domain',                                &
     600                                '&MINVAL( nzb_local ) = ', nzb_local_min,      &
     601                                '&MAXVAL( nzb_local ) = ', nzb_local_max
     602          CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
     603       ENDIF
     604    ENDIF
     605
     606    nzb_s_inner = nzb;  nzb_s_outer = nzb
     607    nzb_u_inner = nzb;  nzb_u_outer = nzb
     608    nzb_v_inner = nzb;  nzb_v_outer = nzb
     609    nzb_w_inner = nzb;  nzb_w_outer = nzb
     610
     611!
     612!-- Define vertical gridpoint from (or to) which on the usual finite difference
     613!-- form (which does not use surface fluxes) is applied
     614    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
     615       nzb_diff = nzb + 2
     616    ELSE
     617       nzb_diff = nzb + 1
     618    ENDIF
     619
     620    nzb_diff_s_inner = nzb_diff;  nzb_diff_s_outer = nzb_diff
     621!
     622!-- Set Neumann conditions for topography. Will be removed soon.
     623    IF ( .NOT. bc_ns_cyc )  THEN
     624       IF ( nys == 0  )  THEN
     625          nzb_local(-1,:)   = nzb_local(0,:)
     626       ELSEIF ( nyn == ny )  THEN
     627          nzb_local(ny+1,:) = nzb_local(ny,:)
     628       ENDIF
     629    ENDIF
     630
     631    IF ( .NOT. bc_lr_cyc )  THEN
     632       IF ( nxl == 0  )  THEN
     633          nzb_local(:,-1)   = nzb_local(:,0)
     634          nzb_local(:,-2)   = nzb_local(:,0)
     635          nzb_local(:,-3)   = nzb_local(:,0)
     636       ELSEIF ( nxr == nx )  THEN
     637          nzb_local(:,nx+1) = nzb_local(:,nx)
     638          nzb_local(:,nx+2) = nzb_local(:,nx)         
     639          nzb_local(:,nx+3) = nzb_local(:,nx)         
     640       ENDIF         
     641    ENDIF
     642!
     643!-- Initialization of 2D index arrays, will be removed soon!
     644!-- Initialize nzb_s_inner and nzb_w_inner
     645    nzb_s_inner = nzb_local
     646    nzb_w_inner = nzb_local
     647
     648!
     649!-- Initialize remaining index arrays:
     650!-- first pre-initialize them with nzb_s_inner...
     651    nzb_u_inner = nzb_s_inner
     652    nzb_u_outer = nzb_s_inner
     653    nzb_v_inner = nzb_s_inner
     654    nzb_v_outer = nzb_s_inner
     655    nzb_w_outer = nzb_s_inner
     656    nzb_s_outer = nzb_s_inner
     657
     658!
     659!-- nzb_s_outer:
     660!-- extend nzb_local east-/westwards first, then north-/southwards
     661    nzb_tmp = nzb_local
     662    DO  j = nys, nyn
     663       DO  i = nxl, nxr
     664          nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i),             &
     665                              nzb_local(j,i+1) )
     666       ENDDO
     667    ENDDO
     668       
     669    CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
     670     
     671    DO  i = nxl, nxr
     672       DO  j = nys, nyn
     673          nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
     674                                  nzb_tmp(j+1,i) )
     675       ENDDO
     676!
     677!--    non-cyclic boundary conditions (overwritten by call of
     678!--    exchange_horiz_2d_int below in case of cyclic boundary conditions)
     679       IF ( nys == 0 )  THEN
     680          j = -1
     681          nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
     682       ENDIF
     683       IF ( nyn == ny )  THEN
     684          j = ny + 1
     685          nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
     686       ENDIF
     687    ENDDO
     688!
     689!-- nzb_w_outer:
     690!-- identical to nzb_s_outer
     691    nzb_w_outer = nzb_s_outer
     692!
     693!-- nzb_u_inner:
     694!-- extend nzb_local rightwards only
     695    nzb_tmp = nzb_local
     696    DO  j = nys, nyn
     697       DO  i = nxl, nxr
     698          nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
     699       ENDDO
     700    ENDDO
     701       
     702    CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
     703       
     704    nzb_u_inner = nzb_tmp
     705!
     706!-- nzb_u_outer:
     707!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
     708    DO  i = nxl, nxr
     709       DO  j = nys, nyn
     710          nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
     711                                  nzb_tmp(j+1,i) )
     712       ENDDO
     713!
     714!--    non-cyclic boundary conditions (overwritten by call of
     715!--    exchange_horiz_2d_int below in case of cyclic boundary conditions)
     716       IF ( nys == 0 )  THEN
     717          j = -1
     718          nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
     719       ENDIF
     720       IF ( nyn == ny )  THEN
     721          j = ny + 1
     722          nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
     723       ENDIF
     724    ENDDO
     725
     726!
     727!-- nzb_v_inner:
     728!-- extend nzb_local northwards only
     729    nzb_tmp = nzb_local
     730    DO  i = nxl, nxr
     731       DO  j = nys, nyn
     732          nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
     733       ENDDO
     734    ENDDO
     735       
     736    CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )     
     737    nzb_v_inner = nzb_tmp
     738
     739!
     740!-- nzb_v_outer:
     741!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
     742    DO  j = nys, nyn
     743       DO  i = nxl, nxr
     744          nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i),                &
     745                                  nzb_tmp(j,i+1) )
     746       ENDDO
     747!
     748!--    non-cyclic boundary conditions (overwritten by call of
     749!--    exchange_horiz_2d_int below in case of cyclic boundary conditions)
     750       IF ( nxl == 0 )  THEN
     751          i = -1
     752          nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
     753       ENDIF
     754       IF ( nxr == nx )  THEN
     755          i = nx + 1
     756          nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
     757       ENDIF
     758    ENDDO
     759
     760!
     761!-- Exchange of lateral boundary values (parallel computers) and cyclic
     762!-- boundary conditions, if applicable.
     763!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
     764!-- they do not require exchange and are not included here.
     765    CALL exchange_horiz_2d_int( nzb_u_inner, nys, nyn, nxl, nxr, nbgp )
     766    CALL exchange_horiz_2d_int( nzb_u_outer, nys, nyn, nxl, nxr, nbgp )
     767    CALL exchange_horiz_2d_int( nzb_v_inner, nys, nyn, nxl, nxr, nbgp )
     768    CALL exchange_horiz_2d_int( nzb_v_outer, nys, nyn, nxl, nxr, nbgp )
     769    CALL exchange_horiz_2d_int( nzb_w_outer, nys, nyn, nxl, nxr, nbgp )
     770    CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp )
     771
     772!
     773!-- Set the individual index arrays which define the k index from which on
     774!-- the usual finite difference form (which does not use surface fluxes) is
     775!-- applied
     776    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
     777       nzb_diff_s_inner   = nzb_s_inner + 2
     778       nzb_diff_s_outer   = nzb_s_outer + 2
     779    ELSE
     780       nzb_diff_s_inner   = nzb_s_inner + 1
     781       nzb_diff_s_outer   = nzb_s_outer + 1
     782    ENDIF
     783!
     784!-- Vertical nesting: communicate vertical grid level arrays between fine and
     785!-- coarse grid
     786    IF ( vnested )  CALL vnest_init_grid
     787
     788 END SUBROUTINE init_grid
     789
     790! Description:
     791! -----------------------------------------------------------------------------!
     792!> Set temporary topography flags and reference buildings on top of underlying
     793!> orography.
     794!------------------------------------------------------------------------------!
     795 SUBROUTINE process_topography( topo_3d )
     796
     797    USE arrays_3d,                                                             &
     798        ONLY:  zw
     799
     800    USE control_parameters,                                                    &
     801        ONLY:  bc_lr_cyc, bc_ns_cyc, land_surface, ocean, urban_surface
     802
     803    USE indices,                                                               &
     804        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
     805               nzt
     806
     807    USE netcdf_data_input_mod,                                                 &
     808        ONLY:  buildings_f, building_id_f, input_pids_static,                  &
     809               terrain_height_f
     810
     811    USE kinds
     812
     813    USE pegrid
     814
     815    IMPLICIT NONE
     816
     817    INTEGER(iwp) ::  i         !< running index along x-direction
     818    INTEGER(iwp) ::  j         !< running index along y-direction
     819    INTEGER(iwp) ::  k         !< running index along z-direction with respect to numeric grid
     820    INTEGER(iwp) ::  k2        !< running index along z-direction with respect to netcdf grid
     821    INTEGER(iwp) ::  k_surf    !< orography top index, used to map 3D buildings onto terrain
     822    INTEGER(iwp) ::  nr        !< index variable indication maximum terrain height for respective building ID
     823    INTEGER(iwp) ::  num_build !< counter for number of buildings
     824
     825    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displace_dum        !< displacements of start addresses, used for MPI_ALLGATHERV
     826    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids           !< building IDs on entire model domain
     827    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final     !< building IDs on entire model domain, multiple occurences are sorted out
     828    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final_tmp !< temporary array used for resizing
     829    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l         !< building IDs on local subdomain
     830    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l_tmp     !< temporary array used to resize array of building IDs
     831
     832    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings     !< number of buildings with different ID on entire model domain
     833    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings_l   !< number of buildings with different ID on local subdomain
     834
     835    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags
     836
     837    REAL(wp)                            ::  ocean_offset        !< offset to consider inverse vertical coordinate at topography definition
     838    REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max             !< maximum terrain height occupied by an building with certain id
     839    REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max_l           !< maximum terrain height occupied by an building with certain id, on local subdomain
     840
     841!
     842!-- In the following, buildings and orography are further preprocessed
     843!-- before they are mapped on the LES grid.
     844!-- Buildings are mapped on top of the orography by maintaining the roof
     845!-- shape of the building. This can be achieved by referencing building on
     846!-- top of the maximum terrain height within the area occupied by the
     847!-- respective building. As buildings and terrain height are defined PE-wise,
     848!-- parallelization of this referencing is required (a building can be
     849!-- distributed between different PEs). 
     850!-- In a first step, determine the number of buildings with different
     851!-- building id on each PE. In a next step, all building ids are gathered
     852!-- into one array which is present to all PEs. For each building ID,
     853!-- the maximum terrain height occupied by the respective building is
     854!-- computed and distributed to each PE. 
     855!-- Finally, for each building id and its respective reference orography,
     856!-- builidings are mapped on top.   
     857!--
     858!-- First, pre-set topography flags, bit 1 indicates orography, bit 2
     859!-- buildings
     860!-- classify the respective surfaces.
     861    topo_3d          = IBSET( topo_3d, 0 )
     862    topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
     863!
     864!-- Reference buildings on top of orography. This is not necessary
     865!-- if topography is read from ASCII file as no distinction between buildings
     866!-- and terrain height can be made. Moreover, this is also not necessary if
     867!-- urban-surface and land-surface model are used at the same time.
     868    IF ( input_pids_static )  THEN
     869       num_buildings_l = 0
     870       num_buildings   = 0
     871!
     872!--    Allocate at least one element for building ids,
     873       ALLOCATE( build_ids_l(1) )
     874       DO  i = nxl, nxr
     875          DO  j = nys, nyn
     876             IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
     877                IF ( num_buildings_l(myid) > 0 )  THEN
     878                   IF ( ANY( building_id_f%var(j,i) .EQ.  build_ids_l ) )  THEN
     879                      CYCLE
     880                   ELSE
     881                      num_buildings_l(myid) = num_buildings_l(myid) + 1
     882!
     883!--                   Resize array with different local building ids
     884                      ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
     885                      build_ids_l_tmp = build_ids_l
     886                      DEALLOCATE( build_ids_l )
     887                      ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
     888                      build_ids_l(1:num_buildings_l(myid)-1) =                 &
     889                                  build_ids_l_tmp(1:num_buildings_l(myid)-1)
     890                      build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
     891                      DEALLOCATE( build_ids_l_tmp )
     892                   ENDIF
     893!
     894!--             First occuring building id on PE
     895                ELSE
     896                   num_buildings_l(myid) = num_buildings_l(myid) + 1
     897                   build_ids_l(1) = building_id_f%var(j,i)
     898                ENDIF
     899             ENDIF
     900          ENDDO
     901       ENDDO
     902!
     903!--    Determine number of different building ids for the entire domain
     904#if defined( __parallel )
     905       CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs,              &
     906                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
     907#else
     908       num_buildings = num_buildings_l
     909#endif
     910!
     911!--    Gather all buildings ids on each PEs.
     912!--    First, allocate array encompassing all building ids in model domain. 
     913       ALLOCATE( build_ids(1:SUM(num_buildings)) )
     914#if defined( __parallel )
     915!
     916!--    Allocate array for displacements.
     917!--    As each PE may has a different number of buildings, so that
     918!--    the block sizes send by each PE may not be equal. Hence,
     919!--    information about the respective displacement is required, indicating
     920!--    the respective adress where each MPI-task writes into the receive
     921!--    buffer array 
     922       ALLOCATE( displace_dum(0:numprocs-1) )
     923       displace_dum(0) = 0
     924       DO i = 1, numprocs-1
     925          displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
     926       ENDDO
     927
     928       CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)),                 &
     929                            num_buildings(myid),                                  &
     930                            MPI_INTEGER,                                          &
     931                            build_ids,                                            &
     932                            num_buildings,                                        &
     933                            displace_dum,                                         &
     934                            MPI_INTEGER,                                          &
     935                            comm2d, ierr )   
     936
     937       DEALLOCATE( displace_dum )
     938
     939#else
     940       build_ids = build_ids_l
     941#endif
     942
     943!
     944!--    Note, in parallel mode building ids can occure mutliple times, as
     945!--    each PE has send its own ids. Therefore, sort out building ids which
     946!--    appear more than one time.
     947       num_build = 0
     948       DO  nr = 1, SIZE(build_ids)
     949
     950          IF ( ALLOCATED(build_ids_final) )  THEN
     951             IF ( ANY( build_ids(nr) .EQ. build_ids_final ) )  THEN
     952                CYCLE
     953             ELSE
     954                num_build = num_build + 1
     955!
     956!--             Resize
     957                ALLOCATE( build_ids_final_tmp(1:num_build) )
     958                build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
     959                DEALLOCATE( build_ids_final )
     960                ALLOCATE( build_ids_final(1:num_build) )
     961                build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
     962                build_ids_final(num_build) = build_ids(nr)
     963                DEALLOCATE( build_ids_final_tmp )
     964             ENDIF             
     965          ELSE
     966             num_build = num_build + 1
     967             ALLOCATE( build_ids_final(1:num_build) )
     968             build_ids_final(num_build) = build_ids(nr)
     969          ENDIF
     970       ENDDO
     971
     972!
     973!--    Finally, determine maximumum terrain height occupied by the
     974!--    respective building. First, on each PE locally.
     975       ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
     976       ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
     977       oro_max_l = 0.0_wp
     978
     979       DO  nr = 1, SIZE(build_ids_final)
     980          oro_max_l(nr) = MAXVAL(                                              &
     981                           MERGE( terrain_height_f%var, 0.0_wp,                &
     982                                  building_id_f%var(nys:nyn,nxl:nxr) .EQ.      &
     983                                  build_ids_final(nr) ) )
     984       ENDDO
     985   
     986#if defined( __parallel )   
     987       IF ( SIZE(build_ids_final) >= 1 ) THEN
     988          CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL,   &
     989                              MPI_MAX, comm2d, ierr )
     990       ENDIF
     991#else
     992       oro_max = oro_max_l
     993#endif
     994!
     995!--    Map orography as well as buildings on grid.
     996!--    In case of ocean simulations, add an offset. 
     997       ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
     998       DO  i = nxl, nxr
     999          DO  j = nys, nyn
     1000             DO  k = nzb, nzt
     1001!
     1002!--             In a first step, if grid point is below or equal the given
     1003!--             terrain height, grid point is flagged to be of type natural.
     1004!--             Please note, in case there is also a building which is lower
     1005!--             than the vertical grid spacing, initialization of surface
     1006!--             attributes will not be correct as given surface information
     1007!--             will not be in accordance to the classified grid points.
     1008!--             Hence, in this case, de-flag the grid point and give it
     1009!--             urban type instead.
     1010                IF ( zw(k) - ocean_offset <= terrain_height_f%var(j,i) )  THEN
     1011                    topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1012                    topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
     1013                ENDIF
     1014!
     1015!--             Set building grid points. Here, only consider 2D buildings.
     1016!--             3D buildings require separate treatment.
     1017                IF ( buildings_f%lod == 1 )  THEN
     1018                   IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
     1019!
     1020!--                   Determine index where maximum terrain height occupied by
     1021!--                   the respective building height is stored.
     1022                      nr = MINLOC( ABS( build_ids_final -                      &
     1023                                        building_id_f%var(j,i) ), DIM = 1 )
     1024       
     1025                      IF ( zw(k) - ocean_offset <=                             &
     1026                           oro_max(nr) + buildings_f%var_2d(j,i) )  THEN
     1027                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1028                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     1029!
     1030!--                      De-flag grid point of type natural. See comment above.
     1031                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
     1032                      ENDIF
     1033                   ENDIF
     1034                ENDIF
     1035             ENDDO
     1036!
     1037!--          Map 3D buildings onto terrain height. 
     1038             IF ( buildings_f%lod == 2 )  THEN
     1039                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
     1040!
     1041!--                Determine topo index occupied by terrain
     1042                   k_surf = MAXLOC( MERGE( 1, 0,                               &
     1043                                           .NOT. BTEST( topo_3d(:,j,i), 0 ) ), &
     1044                                    DIM = 1 ) - 1
     1045                   k2 = nzb+1
     1046                   DO k = k_surf + 1, nzt + 1
     1047                      IF ( k2 <= buildings_f%nz )  THEN
     1048                         IF ( buildings_f%var_3d(k2,j,i) == 1 )  THEN
     1049                            topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1050                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     1051                         ENDIF
     1052                      ENDIF
     1053                      k2 = k2 + 1
     1054                   ENDDO
     1055                ENDIF
     1056             ENDIF
     1057          ENDDO
     1058       ENDDO
     1059!
     1060!--    Deallocate temporary arrays required for processing and reading data
     1061       IF ( ALLOCATED( oro_max         ) )  DEALLOCATE( oro_max         )
     1062       IF ( ALLOCATED( oro_max_l       ) )  DEALLOCATE( oro_max_l       )
     1063       IF ( ALLOCATED( build_ids_final ) )  DEALLOCATE( build_ids_final )
     1064!
     1065!-- Topography input via ASCII format.
     1066    ELSE
     1067       ocean_offset     = MERGE( zw(0), 0.0_wp, ocean )
     1068       topo_3d          = IBSET( topo_3d, 0 )
     1069       topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
     1070       DO  i = nxl, nxr
     1071          DO  j = nys, nyn
     1072             DO  k = nzb, nzt
     1073                IF ( zw(k) - ocean_offset <= buildings_f%var_2d(j,i) )  THEN
     1074                    topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1075                    topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates terrain
     1076                ENDIF
     1077             ENDDO
     1078          ENDDO
     1079       ENDDO
     1080    ENDIF
     1081
     1082    CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
     1083
     1084    IF ( .NOT. bc_ns_cyc )  THEN
     1085       IF ( nys == 0  )  topo_3d(:,-1,:)   = topo_3d(:,0,:)
     1086       IF ( nyn == ny )  topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
     1087    ENDIF
     1088
     1089    IF ( .NOT. bc_lr_cyc )  THEN
     1090       IF ( nxl == 0  )  topo_3d(:,:,-1)   = topo_3d(:,:,0)
     1091       IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)         
     1092    ENDIF
     1093
     1094 END SUBROUTINE process_topography
     1095
     1096
     1097! Description:
     1098! -----------------------------------------------------------------------------!
     1099!> Filter topography, i.e. fill holes resolved by only one grid point. 
     1100!> Such holes are suspected to lead to velocity blow-ups as continuity
     1101!> equation on discrete grid cannot be fulfilled in such case.
     1102!------------------------------------------------------------------------------!
     1103 SUBROUTINE filter_topography( topo_3d )
     1104
     1105    USE control_parameters,                                                    &
     1106        ONLY:  bc_lr_cyc, bc_ns_cyc, message_string
     1107
     1108    USE indices,                                                               &
     1109        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
     1110
     1111    USE netcdf_data_input_mod,                                                 &
     1112        ONLY:  building_id_f, building_type_f
     1113
     1114    USE  pegrid
     1115
     1116    IMPLICIT NONE
     1117
     1118    INTEGER(iwp) ::  i          !< running index along x-direction
     1119    INTEGER(iwp) ::  it         !< counter for number of iterations
     1120    INTEGER(iwp) ::  j          !< running index along y-direction
     1121    INTEGER(iwp) ::  k          !< running index along z-direction
     1122    INTEGER(iwp) ::  num_hole   !< number of holes (in topography) resolved by only one grid point
     1123    INTEGER(iwp) ::  num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE     
     1124    INTEGER(iwp) ::  num_wall   !< number of surrounding vertical walls for a single grid point
     1125
     1126    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo_tmp      !< temporary 3D-topography used to fill holes
     1127    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d       !< 3D-topography array merging buildings and orography
     1128!
     1129!-- Before checking for holes, set lateral boundary conditions for
     1130!-- topography. After hole-filling, boundary conditions must be set again.
     1131!-- Several iterations are performed, in order to fill holes which might
     1132!-- emerge by the filling-algorithm itself.
     1133    ALLOCATE( topo_tmp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1134    topo_tmp = 0
     1135
     1136    num_hole = 99999
     1137    it = 0
     1138    DO WHILE ( num_hole > 0 )       
     1139
     1140       num_hole = 0   
     1141       CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
     1142
     1143       topo_tmp = topo_3d
     1144!
     1145!--    In case of non-cyclic lateral boundaries, assume lateral boundary to be
     1146!--    a solid wall. Thus, intermediate spaces of one grid point between
     1147!--    boundary and some topographic structure will be filled.           
     1148       IF ( .NOT. bc_ns_cyc )  THEN
     1149          IF ( nys == 0  )  topo_tmp(:,-1,:)   = IBCLR( topo_tmp(:,0,:),  0 )
     1150          IF ( nyn == ny )  topo_tmp(:,ny+1,:) = IBCLR( topo_tmp(:,ny,:), 0 )
     1151       ENDIF
     1152
     1153       IF ( .NOT. bc_lr_cyc )  THEN
     1154          IF ( nxl == 0  )  topo_tmp(:,:,-1)   = IBCLR( topo_tmp(:,:,0),  0 )
     1155          IF ( nxr == nx )  topo_tmp(:,:,nx+1) = IBCLR( topo_tmp(:,:,nx), 0 )         
     1156       ENDIF
     1157
     1158       it = it + 1
     1159       num_hole_l = 0
     1160       DO i = nxl, nxr
     1161          DO j = nys, nyn
     1162             DO  k = nzb+1, nzt
     1163                IF ( BTEST( topo_tmp(k,j,i), 0 ) )  THEN
     1164                   num_wall = 0
     1165                   IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) )                  &
     1166                      num_wall = num_wall + 1
     1167                   IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) )                  &
     1168                      num_wall = num_wall + 1
     1169                   IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) )                  &
     1170                      num_wall = num_wall + 1
     1171                   IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) )                  &
     1172                      num_wall = num_wall + 1
     1173                   IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) )                  &
     1174                      num_wall = num_wall + 1   
     1175                   IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) )                  &
     1176                      num_wall = num_wall + 1
     1177
     1178                   IF ( num_wall >= 4 )  THEN
     1179                      num_hole_l     = num_hole_l + 1
     1180!
     1181!--                   Clear flag 0 and set special flag ( bit 3) to indicate
     1182!--                   that new topography point is a result of filtering process.
     1183                      topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1184                      topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 3 )
     1185!
     1186!--                   If filled grid point is occupied by a building, classify
     1187!--                   it as building grid point.
     1188                      IF ( building_type_f%from_file )  THEN
     1189                         IF ( building_type_f%var(j,i)   /=                    & 
     1190                              building_type_f%fill            .OR.             &       
     1191                              building_type_f%var(j+1,i) /=                    & 
     1192                              building_type_f%fill            .OR.             &               
     1193                              building_type_f%var(j-1,i) /=                    &               
     1194                              building_type_f%fill            .OR.             &               
     1195                              building_type_f%var(j,i+1) /=                    &               
     1196                              building_type_f%fill            .OR.             &               
     1197                              building_type_f%var(j,i-1) /=                    &               
     1198                              building_type_f%fill )  THEN
     1199!
     1200!--                         Set flag indicating building surfaces
     1201                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     1202!
     1203!--                         Set building_type and ID at this position if not
     1204!--                         already set. This is required for proper
     1205!--                         initialization of urban-surface energy balance
     1206!--                         solver.
     1207                            IF ( building_type_f%var(j,i) ==                   &
     1208                                 building_type_f%fill )  THEN
     1209
     1210                               IF ( building_type_f%var(j+1,i) /=              &
     1211                                    building_type_f%fill )  THEN
     1212                                  building_type_f%var(j,i) =                   &
     1213                                                    building_type_f%var(j+1,i)
     1214                                  building_id_f%var(j,i) =                     &
     1215                                                    building_id_f%var(j+1,i)
     1216                               ELSEIF ( building_type_f%var(j-1,i) /=          &
     1217                                        building_type_f%fill )  THEN
     1218                                  building_type_f%var(j,i) =                   &
     1219                                                    building_type_f%var(j-1,i)
     1220                                  building_id_f%var(j,i) =                     &
     1221                                                    building_id_f%var(j-1,i)
     1222                               ELSEIF ( building_type_f%var(j,i+1) /=          &
     1223                                        building_type_f%fill )  THEN
     1224                                  building_type_f%var(j,i) =                   &
     1225                                                    building_type_f%var(j,i+1)
     1226                                  building_id_f%var(j,i) =                     &
     1227                                                    building_id_f%var(j,i+1)
     1228                               ELSEIF ( building_type_f%var(j,i-1) /=          &
     1229                                        building_type_f%fill )  THEN
     1230                                  building_type_f%var(j,i) =                   &
     1231                                                    building_type_f%var(j,i-1)
     1232                                  building_id_f%var(j,i) =                     &
     1233                                                    building_id_f%var(j,i-1)
     1234                               ENDIF
     1235                            ENDIF
     1236                         ENDIF
     1237                      ENDIF
     1238!
     1239!--                   If filled grid point is already classified as building
     1240!--                   everything is fine, else classify this grid point as
     1241!--                   natural type grid point. This case, values for the
     1242!--                   surface type are already set.
     1243                      IF ( .NOT. BTEST( topo_3d(k,j,i), 2 ) )  THEN
     1244                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
     1245                      ENDIF
     1246                   ENDIF
     1247                ENDIF
     1248             ENDDO
     1249          ENDDO
     1250       ENDDO
     1251!
     1252!--    Count the total number of holes, required for informative message.
     1253#if defined( __parallel )
     1254       CALL MPI_ALLREDUCE( num_hole_l, num_hole, 1, MPI_INTEGER, MPI_SUM,      &
     1255                           comm2d, ierr )
     1256#else
     1257       num_hole = num_hole_l
     1258#endif   
     1259
     1260!
     1261!--    Create an informative message if any holes were filled.
     1262       IF ( num_hole > 0 )  THEN
     1263          WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//   &
     1264                                     'one grid point were filled at ', it,     &
     1265                                     'th iteration '
     1266          CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 )
     1267       ENDIF
     1268
     1269
     1270    ENDDO
     1271
     1272    DEALLOCATE( topo_tmp )
     1273!
     1274!-- Finally, exchange topo_3d array again and if necessary set Neumann boundary
     1275!-- condition in case of non-cyclic lateral boundaries.
     1276    CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
     1277
     1278    IF ( .NOT. bc_ns_cyc )  THEN
     1279       IF ( nys == 0  )  topo_3d(:,-1,:)   = topo_3d(:,0,:)
     1280       IF ( nyn == ny )  topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
     1281    ENDIF
     1282
     1283    IF ( .NOT. bc_lr_cyc )  THEN
     1284       IF ( nxl == 0  )  topo_3d(:,:,-1)   = topo_3d(:,:,0)
     1285       IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)         
     1286    ENDIF
     1287
     1288 END SUBROUTINE filter_topography
     1289
     1290! Description:
     1291! -----------------------------------------------------------------------------!
     1292!> Pre-computation of grid-dependent and near-wall mixing length.
     1293!> @todo: move subroutine to turbulence module developed by M1 (Tobi).
     1294!------------------------------------------------------------------------------!
     1295 SUBROUTINE init_mixing_length
     1296
     1297    USE arrays_3d,                                                             &
     1298        ONLY:  dzw, l_grid, l_wall, zu, zw
     1299
     1300    USE control_parameters,                                                    &
     1301        ONLY:  message_string, wall_adjustment_factor
     1302
     1303    USE grid_variables,                                                        &
     1304        ONLY:  dx, dy
     1305
     1306    USE indices,                                                               &
     1307        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
     1308               wall_flags_0
     1309   
     1310    USE kinds
     1311
     1312    IMPLICIT NONE
     1313
     1314    INTEGER(iwp) ::  i             !< index variable along x
     1315    INTEGER(iwp) ::  j             !< index variable along y
     1316    INTEGER(iwp) ::  k             !< index variable along z
     1317
     1318    ALLOCATE( l_grid(1:nzt) )
     1319    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1320!
     1321!-- Compute the grid-dependent mixing length.
     1322    DO  k = 1, nzt
     1323       l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
     1324    ENDDO
     1325!
     1326!-- Initialize near-wall mixing length l_wall only in the vertical direction
     1327!-- for the moment, multiplication with wall_adjustment_factor further below
     1328    l_wall(nzb,:,:)   = l_grid(1)
     1329    DO  k = nzb+1, nzt
     1330       l_wall(k,:,:)  = l_grid(k)
     1331    ENDDO
     1332    l_wall(nzt+1,:,:) = l_grid(nzt)
     1333
     1334    DO  k = 1, nzt
     1335       IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.  &
     1336            l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
     1337          WRITE( message_string, * ) 'grid anisotropy exceeds ', &
     1338                                     'threshold given by only local', &
     1339                                     ' &horizontal reduction of near_wall ', &
     1340                                     'mixing length l_wall', &
     1341                                     ' &starting from height level k = ', k, '.'
     1342          CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
     1343          EXIT
     1344       ENDIF
     1345    ENDDO
     1346!
     1347!-- In case of topography: limit near-wall mixing length l_wall further:
     1348!-- Go through all points of the subdomain one by one and look for the closest
     1349!-- surface.
     1350!-- Is this correct in the ocean case?
     1351    DO  i = nxl, nxr
     1352       DO  j = nys, nyn
     1353          DO  k = nzb+1, nzt
     1354!
     1355!--          Check if current gridpoint belongs to the atmosphere
     1356             IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
     1357!
     1358!--             Check for neighbouring grid-points.
     1359!--             Vertical distance, down
     1360                IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )                 &
     1361                   l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
     1362!
     1363!--             Vertical distance, up
     1364                IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )                 &
     1365                   l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
     1366!
     1367!--             y-distance
     1368                IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .OR.             &
     1369                     .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )                 &
     1370                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
     1371!
     1372!--             x-distance
     1373                IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .OR.             &
     1374                     .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )                 &
     1375                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
     1376!
     1377!--              yz-distance (vertical edges, down)
     1378                 IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 )  .OR.          &
     1379                      .NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 )  )             &
     1380                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1381                                        SQRT( 0.25_wp * dy**2 +                &
     1382                                       ( zu(k) - zw(k-1) )**2 ) )
     1383!
     1384!--              yz-distance (vertical edges, up)
     1385                 IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 )  .OR.          &
     1386                      .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 )  )             &
     1387                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1388                                        SQRT( 0.25_wp * dy**2 +                &
     1389                                       ( zw(k) - zu(k) )**2 ) )
     1390!   
     1391!--              xz-distance (vertical edges, down)
     1392                 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 )  .OR.          &
     1393                      .NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 )  )             &
     1394                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1395                                        SQRT( 0.25_wp * dx**2 +                &
     1396                                       ( zu(k) - zw(k-1) )**2 ) )
     1397!
     1398!--              xz-distance (vertical edges, up)
     1399                 IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 )  .OR.          &
     1400                      .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 )  )             &
     1401                  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1402                                        SQRT( 0.25_wp * dx**2 +                &
     1403                                       ( zw(k) - zu(k) )**2 ) )
     1404!
     1405!--             xy-distance (horizontal edges)
     1406                IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 )  .OR.           &
     1407                     .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.           &
     1408                     .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.           &
     1409                     .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) )               &
     1410                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1411                                        SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
     1412!
     1413!--             xyz distance (vertical and horizontal edges, down)
     1414                IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 )  .OR.         &
     1415                     .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.         &
     1416                     .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.         &
     1417                     .NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) )             &
     1418                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1419                                        SQRT( 0.25_wp * ( dx**2 + dy**2 )      &
     1420                                              +  ( zu(k) - zw(k-1) )**2  ) )
     1421!
     1422!--             xyz distance (vertical and horizontal edges, up)
     1423                IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 )  .OR.         &
     1424                     .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.         &
     1425                     .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.         &
     1426                     .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) )             &
     1427                   l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
     1428                                        SQRT( 0.25_wp * ( dx**2 + dy**2 )      &
     1429                                              +  ( zw(k) - zu(k) )**2  ) )
     1430                 
     1431             ENDIF
     1432          ENDDO
     1433       ENDDO
     1434    ENDDO
     1435!
     1436!-- Set lateral boundary conditions for l_wall
     1437    CALL exchange_horiz( l_wall, nbgp )     
     1438
     1439 END SUBROUTINE init_mixing_length
     1440
     1441! Description:
     1442! -----------------------------------------------------------------------------!
     1443!> Reads topography information from file or sets generic topography. Moreover,
     1444!> all topography-relevant topography arrays are initialized, and grid flags
     1445!> are set. 
     1446!------------------------------------------------------------------------------!
     1447 SUBROUTINE init_topo( topo )
     1448
     1449    USE arrays_3d,                                                             &
     1450        ONLY:  zw
     1451       
     1452    USE control_parameters,                                                    &
     1453        ONLY:  bc_lr_cyc, bc_ns_cyc, building_height, building_length_x,       &
     1454               building_length_y, building_wall_left, building_wall_south,     &
     1455               canyon_height, canyon_wall_left, canyon_wall_south,             &
     1456               canyon_width_x, canyon_width_y, dp_level_ind_b, dz,             &
     1457               message_string, ocean, topography, topography_grid_convention,  &
     1458               tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,   &
     1459               tunnel_wall_depth
     1460         
     1461    USE grid_variables,                                                        &
     1462        ONLY:  dx, dy
     1463       
     1464    USE indices,                                                               &
     1465        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
     1466               nzb, nzt
     1467   
     1468    USE kinds
     1469
     1470    USE pegrid
     1471
     1472    USE surface_mod,                                                           &
     1473        ONLY:  get_topography_top_index
    3001474
    3011475    IMPLICIT NONE
     
    3161490    INTEGER(iwp) ::  cys           !< index for south canyon wall
    3171491    INTEGER(iwp) ::  i             !< index variable along x
    318     INTEGER(iwp) ::  id_topo       !< NetCDF id of topograhy input file
    319     INTEGER(iwp) ::  ii            !< loop variable for reading topography file
    320     INTEGER(iwp) ::  inc           !< incremental parameter for coarsening grid level
     1492    INTEGER(iwp) ::  it            !< Number of iterations required to fill all problematic holes
    3211493    INTEGER(iwp) ::  j             !< index variable along y
    3221494    INTEGER(iwp) ::  k             !< index variable along z
    323     INTEGER(iwp) ::  k_top         !< topography top index on local PE
    324     INTEGER(iwp) ::  l             !< loop variable
    325     INTEGER(iwp) ::  nxl_l         !< index of left PE boundary for multigrid level
    326     INTEGER(iwp) ::  nxr_l         !< index of right PE boundary for multigrid level
    327     INTEGER(iwp) ::  nyn_l         !< index of north PE boundary for multigrid level
    328     INTEGER(iwp) ::  nys_l         !< index of south PE boundary for multigrid level
    329     INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
    330     INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
    331     INTEGER(iwp) ::  nzt_l         !< index of top PE boundary for multigrid level
    332     INTEGER(iwp) ::  num_hole      !< number of holes (in topography) resolved by only one grid point
    333     INTEGER(iwp) ::  num_hole_l    !< number of holes (in topography) resolved by only one grid point on local PE     
    334     INTEGER(iwp) ::  num_wall      !< number of surrounding vertical walls for a single grid point
    335     INTEGER(iwp) ::  skip_n_rows   !< counting variable to skip rows while reading topography file   
    3361495    INTEGER(iwp) ::  hv_in         !< heavyside function to model inner tunnel surface
    3371496    INTEGER(iwp) ::  hv_out        !< heavyside function to model outer tunnel surface
     
    3461505    INTEGER(iwp) ::  td            !< tunnel wall depth
    3471506    INTEGER(iwp) ::  th            !< height of outer tunnel wall
    348                                      
    349     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
    350     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp        !< dummy to calculate topography indices on u- and v-grid
    351 
    352     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags
    353 
    354     LOGICAL  ::  netcdf_extend = .FALSE. !< Flag indicating wether netcdf topography input file or not
    355 
    356     REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file   
    357     REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
    358 
    359     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  oro_height    !< input variable for terrain height
    360     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  topo_height   !< input variable for topography height
    361 
    362     INTEGER( KIND=1 ), DIMENSION(:,:,:), ALLOCATABLE ::  topo_3d_read !< input variable for 3D topography
    363 
    364 !
    365 !-- Calculation of horizontal array bounds including ghost layers
    366     nxlg = nxl - nbgp
    367     nxrg = nxr + nbgp
    368     nysg = nys - nbgp
    369     nyng = nyn + nbgp
    370 
    371 !
    372 !-- Allocate grid arrays
    373     ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1),        &
    374               dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
    375 
    376 !
    377 !-- Compute height of u-levels from constant grid length and dz stretch factors
    378     IF ( dz == -1.0_wp )  THEN
    379        message_string = 'missing dz'
    380        CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
    381     ELSEIF ( dz <= 0.0_wp )  THEN
    382        WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
    383        CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
    384     ENDIF
    385 
    386 !
    387 !-- Define the vertical grid levels
    388     IF ( .NOT. ocean )  THEN
    389 !
    390 !--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
    391 !--    The second u-level (k=1) corresponds to the top of the
    392 !--    Prandtl-layer.
    393 
    394        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
    395           zu(0) = 0.0_wp
    396       !    zu(0) = - dz * 0.5_wp
    397        ELSE
    398           zu(0) = - dz * 0.5_wp
    399        ENDIF
    400        zu(1) =   dz * 0.5_wp
    401 
    402        dz_stretch_level_index = nzt+1
    403        dz_stretched = dz
    404        DO  k = 2, nzt+1
    405           IF ( dz_stretch_level <= zu(k-1)  .AND.  dz_stretched < dz_max )  THEN
    406              dz_stretched = dz_stretched * dz_stretch_factor
    407              dz_stretched = MIN( dz_stretched, dz_max )
    408              IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
    409           ENDIF
    410           zu(k) = zu(k-1) + dz_stretched
    411        ENDDO
    412 
    413 !
    414 !--    Compute the w-levels. They are always staggered half-way between the
    415 !--    corresponding u-levels. In case of dirichlet bc for u and v at the
    416 !--    ground the first u- and w-level (k=0) are defined at same height (z=0).
    417 !--    The top w-level is extrapolated linearly.
    418        zw(0) = 0.0_wp
    419        DO  k = 1, nzt
    420           zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
    421        ENDDO
    422        zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
    423 
    424     ELSE
    425 !
    426 !--    Grid for ocean with free water surface is at k=nzt (w-grid).
    427 !--    In case of neumann bc at the ground the first first u-level (k=0) lies
    428 !--    below the first w-level (k=0). In case of dirichlet bc the first u- and
    429 !--    w-level are defined at same height, but staggered from the second level.
    430 !--    The second u-level (k=1) corresponds to the top of the Prandtl-layer.
    431        zu(nzt+1) =   dz * 0.5_wp
    432        zu(nzt)   = - dz * 0.5_wp
    433 
    434        dz_stretch_level_index = 0
    435        dz_stretched = dz
    436        DO  k = nzt-1, 0, -1
    437 !
    438 !--       The default value of dz_stretch_level is positive, thus the first
    439 !--       condition is always true. Hence, the second condition is necessary.
    440           IF ( dz_stretch_level >= zu(k+1)  .AND.  dz_stretch_level <= 0.0  &
    441                .AND.  dz_stretched < dz_max )  THEN
    442              dz_stretched = dz_stretched * dz_stretch_factor
    443              dz_stretched = MIN( dz_stretched, dz_max )
    444              IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
    445           ENDIF
    446           zu(k) = zu(k+1) - dz_stretched
    447        ENDDO
    448 
    449 !
    450 !--    Compute the w-levels. They are always staggered half-way between the
    451 !--    corresponding u-levels, except in case of dirichlet bc for u and v
    452 !--    at the ground. In this case the first u- and w-level are defined at
    453 !--    same height. The top w-level (nzt+1) is not used but set for
    454 !--    consistency, since w and all scalar variables are defined up tp nzt+1.
    455        zw(nzt+1) = dz
    456        zw(nzt)   = 0.0_wp
    457        DO  k = 0, nzt
    458           zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
    459        ENDDO
    460 
    461 !
    462 !--    In case of dirichlet bc for u and v the first u- and w-level are defined
    463 !--    at same height.
    464        IF ( ibc_uv_b == 0 ) THEN
    465           zu(0) = zw(0)
    466        ENDIF
    467 
    468     ENDIF
    469 
    470 !
    471 !-- Compute grid lengths.
    472     DO  k = 1, nzt+1
    473        dzu(k)  = zu(k) - zu(k-1)
    474        ddzu(k) = 1.0_wp / dzu(k)
    475        dzw(k)  = zw(k) - zw(k-1)
    476        ddzw(k) = 1.0_wp / dzw(k)
    477     ENDDO
    478 
    479     DO  k = 1, nzt
    480        dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
    481     ENDDO
    482    
    483 !   
    484 !-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
    485 !-- everywhere. For the actual grid, the grid spacing at the lowest level
    486 !-- is only dz/2, but should be dz. Therefore, an additional array
    487 !-- containing with appropriate grid information is created for these
    488 !-- solvers.
    489     IF ( psolver(1:9) /= 'multigrid' )  THEN
    490        ALLOCATE( ddzu_pres(1:nzt+1) )
    491        ddzu_pres = ddzu
    492        ddzu_pres(1) = ddzu_pres(2)  ! change for lowest level
    493     ENDIF
    494 
    495 !
    496 !-- Compute the reciprocal values of the horizontal grid lengths.
    497     ddx = 1.0_wp / dx
    498     ddy = 1.0_wp / dy
    499     dx2 = dx * dx
    500     dy2 = dy * dy
    501     ddx2 = 1.0_wp / dx2
    502     ddy2 = 1.0_wp / dy2
    503 
    504 !
    505 !-- Compute the grid-dependent mixing length.
    506     DO  k = 1, nzt
    507        l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
    508     ENDDO
    509 
    510 !
    511 !-- Allocate outer and inner index arrays for topography and set
    512 !-- defaults.                   
    513     ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg),                             &
    514               nzb_s_outer(nysg:nyng,nxlg:nxrg),                             &
    515               nzb_u_inner(nysg:nyng,nxlg:nxrg),                             &
    516               nzb_u_outer(nysg:nyng,nxlg:nxrg),                             &
    517               nzb_v_inner(nysg:nyng,nxlg:nxrg),                             &
    518               nzb_v_outer(nysg:nyng,nxlg:nxrg),                             &
    519               nzb_w_inner(nysg:nyng,nxlg:nxrg),                             &
    520               nzb_w_outer(nysg:nyng,nxlg:nxrg),                             &
    521               nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                        &
    522               nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                        &
    523               nzb_local(nysg:nyng,nxlg:nxrg),                               &
    524               nzb_tmp(nysg:nyng,nxlg:nxrg),                                 &
    525               wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    526 
    527     ALLOCATE( topo_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    528     topo_3d    = 0
    529 
    530     ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    531 
    532     nzb_s_inner = nzb;  nzb_s_outer = nzb
    533     nzb_u_inner = nzb;  nzb_u_outer = nzb
    534     nzb_v_inner = nzb;  nzb_v_outer = nzb
    535     nzb_w_inner = nzb;  nzb_w_outer = nzb
    536 
    537 !
    538 !-- Define vertical gridpoint from (or to) which on the usual finite difference
    539 !-- form (which does not use surface fluxes) is applied
    540     IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
    541        nzb_diff = nzb + 2
    542     ELSE
    543        nzb_diff = nzb + 1
    544     ENDIF
    545 
    546     nzb_diff_s_inner = nzb_diff;  nzb_diff_s_outer = nzb_diff
    547 
    548 !
    549 !-- Initialize near-wall mixing length l_wall only in the vertical direction
    550 !-- for the moment,
    551 !-- multiplication with wall_adjustment_factor near the end of this routine
    552     l_wall(nzb,:,:)   = l_grid(1)
    553     DO  k = nzb+1, nzt
    554        l_wall(k,:,:)  = l_grid(k)
    555     ENDDO
    556     l_wall(nzt+1,:,:) = l_grid(nzt)
    557 
    558     DO  k = 1, nzt
    559        IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.  &
    560             l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
    561           WRITE( message_string, * ) 'grid anisotropy exceeds ', &
    562                                      'threshold given by only local', &
    563                                      ' &horizontal reduction of near_wall ', &
    564                                      'mixing length l_wall', &
    565                                      ' &starting from height level k = ', k, '.'
    566           CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
    567           EXIT
    568        ENDIF
    569     ENDDO
     1507
     1508    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local         !< index for topography top at cell-center
     1509    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
     1510
     1511
    5701512!
    5711513!-- Set outer and inner index arrays for non-flat topography.
     
    5781520
    5791521       CASE ( 'flat' )
    580 !
    581 !--       nzb_local is required for the multigrid solver
    582           nzb_local = 0
    583 !
     1522!   
    5841523!--       Initialilize 3D topography array, used later for initializing flags
    585           topo_3d(nzb+1:nzt+1,:,:) = IBSET( topo_3d(nzb+1:nzt+1,:,:), 0 )
    586 !
    587 !--       level of detail is required for output routines
    588           lod       = 1
     1524          topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
    5891525
    5901526       CASE ( 'single_building' )
     
    5971533          IF ( ABS( zw(bh)   - building_height ) == &
    5981534               ABS( zw(bh+1) - building_height )    )  bh = bh + 1
    599 
    6001535          IF ( building_wall_left == 9999999.9_wp )  THEN
    6011536             building_wall_left = ( nx + 1 - blx ) / 2 * dx
     
    6051540
    6061541          IF ( building_wall_south == 9999999.9_wp )  THEN
    607              building_wall_south = ( ny + 1 - bly ) / 2 * dy
     1542              building_wall_south = ( ny + 1 - bly ) / 2 * dy
    6081543          ENDIF
    6091544          bys = NINT( building_wall_south / dy )
     
    6121547!
    6131548!--       Building size has to meet some requirements
    614           IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR.  &
     1549          IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR.       &
    6151550               ( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) )  THEN
    6161551             WRITE( message_string, * ) 'inconsistent building parameters:',   &
     
    6201555          ENDIF
    6211556
     1557          ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
    6221558!
    6231559!--       Define the building.
    624           nzb_local = 0
    6251560          IF ( bxl <= nxr  .AND.  bxr >= nxl  .AND.                            &
    626                bys <= nyn  .AND.  byn >= nys )                                 &        
     1561               bys <= nyn  .AND.  byn >= nys )                                 &
    6271562             nzb_local(MAX(nys,bys):MIN(nyn,byn),MAX(nxl,bxl):MIN(nxr,bxr)) = bh
    628 
    629           CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    630 !
    631 !--       Set bit array to mask topography
    632           DO  i = nxlg, nxrg
    633              DO  j = nysg, nyng
    634 
    635                 topo_3d(nzb_local(j,i)+1:nzt+1,j,i) =                          &
    636                                  IBSET( topo_3d(nzb_local(j,i)+1:nzt+1,j,i), 0 )
     1563!
     1564!--       Set bit array on basis of nzb_local
     1565          DO  i = nxl, nxr
     1566             DO  j = nys, nyn
     1567                topo(nzb_local(j,i)+1:nzt+1,j,i) =                             &
     1568                                 IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 )
    6371569             ENDDO
    6381570          ENDDO
    639 !
    640 !--       level of detail is required for output routines. Here, 2D topography.
    641           lod = 1
    642 
    643           CALL exchange_horiz_int( topo_3d, nbgp )
     1571       
     1572          DEALLOCATE( nzb_local )
     1573
     1574          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
    6441575
    6451576       CASE ( 'single_street_canyon' )
     
    6561587             cxl = NINT( canyon_wall_left / dx )
    6571588             cxr = cxl + cwx
    658 
    6591589          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
    6601590!
     
    6661596             cys = NINT( canyon_wall_south / dy )
    6671597             cyn = cys + cwy
    668 
     1598     
    6691599          ELSE
    6701600             
     
    6771607          IF ( ABS( zw(ch)   - canyon_height ) == &
    6781608               ABS( zw(ch+1) - canyon_height )    )  ch = ch + 1
    679 
    6801609          dp_level_ind_b = ch
    6811610!
     
    6831612          IF ( canyon_width_x /= 9999999.9_wp )  THEN
    6841613             IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.        &
    685                ( ch < 3 ) )  THEN
     1614                  ( ch < 3 ) )  THEN
    6861615                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
    6871616                                           '&cxl=', cxl, 'cxr=', cxr,          &
     
    6921621          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
    6931622             IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.        &
    694                ( ch < 3 ) )  THEN
     1623                  ( ch < 3 ) )  THEN
    6951624                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
    6961625                                           '&cys=', cys, 'cyn=', cyn,          &
     
    7081637          ENDIF
    7091638
     1639          ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )
    7101640          nzb_local = ch
    7111641          IF ( canyon_width_x /= 9999999.9_wp )  THEN
     
    7161646                nzb_local(MAX(nys,cys+1):MIN(nyn,cyn-1),:) = 0
    7171647          ENDIF
    718 
    719           CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    720 !
    721 !--       Set bit array to mask topography
    722           DO  i = nxlg, nxrg
    723              DO  j = nysg, nyng
    724                 topo_3d(nzb_local(j,i)+1:nzt+1,j,i) =                          &
    725                                  IBSET( topo_3d(nzb_local(j,i)+1:nzt+1,j,i), 0 )
     1648!
     1649!--       Set bit array on basis of nzb_local
     1650          DO  i = nxl, nxr
     1651             DO  j = nys, nyn
     1652                topo(nzb_local(j,i)+1:nzt+1,j,i) =                             &
     1653                                 IBSET( topo(nzb_local(j,i)+1:nzt+1,j,i), 0 )
    7261654             ENDDO
    7271655          ENDDO
    728 !
    729 !--       level of detail is required for output routines. Here, 2D topography.
    730           lod = 1
    731 
    732           CALL exchange_horiz_int( topo_3d, nbgp )
     1656          DEALLOCATE( nzb_local )
     1657
     1658          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
    7331659
    7341660       CASE ( 'tunnel' )
     
    7431669!
    7441670!--       Tunnel-wall depth
    745           IF ( tunnel_wall_depth == 9999999.9_wp )  THEN
     1671          IF ( tunnel_wall_depth == 9999999.9_wp )  THEN 
    7461672             td = MAX ( dx, dy, dz )
    7471673          ELSE
     
    7751701                                      ( tunnel_width_x * 0.5_wp - td ) )
    7761702             txe_in  = INT( ( nx + 1 ) * 0.5_wp * dx +                         &
    777                                       ( tunnel_width_x * 0.5_wp - td ) )
     1703                                   ( tunnel_width_x * 0.5_wp - td ) )
    7781704
    7791705             tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_length * 0.5_wp )
     
    7821708             tye_in  = tye_out
    7831709          ENDIF
    784           IF ( tunnel_width_x /= 9999999.9_wp  .AND.                           &
    785                tunnel_width_x - 2.0_wp * tunnel_wall_depth <= 2.0_wp * dx )  THEN
     1710          IF ( tunnel_width_x /= 9999999.9_wp  .AND.                           &   
     1711               tunnel_width_x - 2.0_wp * td <= 2.0_wp * dx )                   &
     1712          THEN
    7861713             message_string = 'Tunnel width too small'
    7871714             CALL message( 'init_grid', 'PA0175', 1, 2, 0, 6, 0 )
    7881715          ENDIF
    7891716          IF ( tunnel_width_y /= 9999999.9_wp  .AND.                           &
    790                tunnel_width_y - 2.0_wp * tunnel_wall_depth <= 2.0_wp * dy )  THEN
     1717               tunnel_width_y - 2.0_wp * td <= 2.0_wp * dy )                   &
     1718          THEN
    7911719             message_string = 'Tunnel width too small'
    7921720             CALL message( 'init_grid', 'PA0455', 1, 2, 0, 6, 0 )
     
    8081736             tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_width_y * 0.5_wp )
    8091737             tys_in  = INT( ( ny + 1 ) * 0.5_wp * dy -                         &
    810                                      ( tunnel_width_y * 0.5_wp - td ) )
     1738                                        ( tunnel_width_y * 0.5_wp - td ) )
    8111739             tye_in  = INT( ( ny + 1 ) * 0.5_wp * dy +                         &
    8121740                                     ( tunnel_width_y * 0.5_wp - td ) )
    8131741          ENDIF
    8141742
    815           topo_3d = 0
     1743          topo = 0
    8161744          DO  i = nxl, nxr
    8171745             DO  j = nys, nyn
     
    8251753                            ( ( SIGN( 1.0_wp, j * dy - tys_out ) + 1.0_wp )    &
    8261754                            - ( SIGN( 1.0_wp, j * dy - tye_out ) + 1.0_wp ) )
    827 !
     1755!    
    8281756!--             Use heaviside function to model inner tunnel surface
    8291757                hv_in  = ( th - td ) * 0.5_wp *                                &
     
    8371765!--             Set flags at x-y-positions without any tunnel surface
    8381766                IF ( hv_out - hv_in == 0.0_wp )  THEN
    839                    topo_3d(nzb+1:nzt+1,j,i) = IBSET( topo_3d(nzb+1:nzt+1,j,i), 0 )
     1767                   topo(nzb+1:nzt+1,j,i) = IBSET( topo(nzb+1:nzt+1,j,i), 0 )
    8401768!
    8411769!--             Set flags at x-y-positions with tunnel surfaces
     
    8461774                      IF ( hv_out - hv_in == th )  THEN
    8471775                         IF ( zw(k) <= hv_out )  THEN
    848                             topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1776                            topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
    8491777                         ELSE
    850                             topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 )
     1778                            topo(k,j,i) = IBSET( topo(k,j,i), 0 )
    8511779                         ENDIF
    8521780                      ENDIF
     
    8551783                      IF ( hv_out - hv_in == td )  THEN
    8561784                         IF ( zw(k) <= hv_in )  THEN
    857                             topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 )
     1785                            topo(k,j,i) = IBSET( topo(k,j,i), 0 )
    8581786                         ELSEIF ( zw(k) > hv_in  .AND.  zw(k) <= hv_out )  THEN
    859                             topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1787                            topo(k,j,i) = IBCLR( topo(k,j,i), 0 )
    8601788                         ELSEIF ( zw(k) > hv_out )  THEN
    861                             topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 )
     1789                            topo(k,j,i) = IBSET( topo(k,j,i), 0 )
    8621790                         ENDIF
    8631791                      ENDIF
     
    8671795          ENDDO
    8681796
    869           nzb_local = 0
    870 !
    871 !--       level of detail is required for output routines. Here, 3D topography.
    872           lod = 2
    873 
    874           CALL exchange_horiz_int( topo_3d, nbgp )
     1797          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
    8751798
    8761799       CASE ( 'read_from_file' )
    877 
    878           ALLOCATE ( oro_height(nys:nyn,nxl:nxr)  )
    879           ALLOCATE ( topo_height(nys:nyn,nxl:nxr) )
    880           oro_height  = 0.0_wp
    881           topo_height = 0.0_wp
    882 
    883           DO  ii = 0, io_blocks-1
    884              IF ( ii == io_group )  THEN
    885 
    886 !
    887 !--             Arbitrary irregular topography data in PALM format (exactly
    888 !--             matching the grid size and total domain size).
    889 !--             First, check if NetCDF file for topography exist or not.
    890 !--             This case, read topography from NetCDF, else read it from
    891 !--             ASCII file.
    892 #if defined ( __netcdf )
    893                 INQUIRE( FILE='TOPOGRAPHY_DATA_NC'//TRIM( coupling_char ),     &
    894                          EXIST=netcdf_extend )
    895 !
    896 !--             NetCDF branch   
    897                 IF ( netcdf_extend )  THEN
    898 !
    899 !--                Open file in read-only mode
    900                    CALL netcdf_open_read_file( 'TOPOGRAPHY_DATA_NC',           &
    901                                                id_topo, 20 )  !Error number still need to be set properly
    902 
    903 !
    904 !--                Read terrain height. Reading is done PE-wise, i.e. each
    905 !--                processor reads its own domain. Reading is realized
    906 !--                via looping over x-dimension, i.e. calling
    907 !--                netcdf_get_variable reads topography along y for given x.
    908 !--                Orography is 2D.
    909                    DO  i = nxl, nxr
    910                       CALL netcdf_get_variable( id_topo, 'orography_0',        &
    911                                                 i, oro_height(:,i), 20 )  !Error number still need to be set properly
    912                    ENDDO
    913 !
    914 !--                Read attribute lod (level of detail), required for variable
    915 !--                buildings_0
    916                    CALL netcdf_get_attribute( id_topo, "lod", lod, .FALSE.,    &
    917                                               20, 'buildings_0' ) !Error number still need to be set properly
    918 !
    919 !--                Read building height
    920 !--                2D for lod = 1, 3D for lod = 2
    921                    IF ( lod == 1 )  THEN
    922                       DO  i = nxl, nxr
    923                          CALL netcdf_get_variable( id_topo, 'buildings_0',     &
    924                                                    i, topo_height(:,i), 20 )  !Error number still need to be set properly
    925                       ENDDO
    926 
    927                    ELSEIF ( lod == 2 )  THEN
    928 !
    929 !--                   Allocate 1-byte integer dummy array to read 3D topography
    930                       ALLOCATE( topo_3d_read(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    931 !
    932 !--                   Read data PE-wise. Read yz-slices.
    933                       DO  i = nxl, nxr
    934                          DO  j = nys, nyn
    935                             CALL netcdf_get_variable( id_topo, 'buildings_0',  &
    936                                                       i, j, topo_3d_read(:,j,i), 20 )  !Error number still need to be set properly
    937                          ENDDO
    938                       ENDDO
    939                    ELSE
    940                       message_string = 'NetCDF attribute lod ' //              &
    941                                        '(level of detail) is not set properly.'
    942                       CALL message( 'init_grid', 'PA0457', 1, 2, 0, 6, 0 )
    943                    ENDIF
    944 !
    945 !--                On file, 3D topography grid points is classified with 1,
    946 !--                atmosphere is classified with 0, contrary to the internal
    947 !--                treatment. Hence, conversion is required. Moreover, set
    948 !--                topography array to zero at lowest grid level.
    949                    topo_3d = MERGE( 0, 1, topo_3d_read == 1 )
    950                    topo_3d(nzb,:,:) = 0
    951 !
    952 !--                Deallocate dummy array
    953                    DEALLOCATE( topo_3d_read )
    954 !
    955 !--                Close topography input file
    956                    CALL netcdf_close_file( id_topo, 20 )
    957 #endif
    958 
    959 !
    960 !--             ASCII branch. Please note, reading of 3D topography is not
    961 !--             supported in ASCII format. Further, no distinction is made
    962 !--             between orography and buildings
    963                 ELSE
    964 
    965                    OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ),    &
    966                          STATUS='OLD', FORM='FORMATTED', ERR=10 )
    967 !
    968 !--                Read topography PE-wise. Rows are read from nyn to nys, columns
    969 !--                are read from nxl to nxr. At first, ny-nyn rows need to be skipped.
    970                    skip_n_rows = 0
    971                    DO WHILE ( skip_n_rows < ny - nyn )
    972                       READ( 90, * ) 
    973                       skip_n_rows = skip_n_rows + 1
    974                    ENDDO
    975 !
    976 !--                Read data from nyn to nys and nxl to nxr. Therefore, skip
    977 !--                column until nxl-1 is reached
    978                    DO  j = nyn, nys, -1
    979                       READ( 90, *, ERR=11, END=11 )                            &
    980                                               ( dum, i = 0, nxl-1 ),           &
    981                                               ( topo_height(j,i), i = nxl, nxr )
    982                    ENDDO
    983 
    984                    GOTO 12
    985          
    986  10                message_string = 'file TOPOGRAPHY'//TRIM( coupling_char )// &
    987                                     ' does not exist'
    988                    CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 )
    989 
    990  11                message_string = 'errors in file TOPOGRAPHY_DATA'//         &
    991                                     TRIM( coupling_char )
    992                    CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 )
    993 
    994  12                CLOSE( 90 )
    995 
    996                 ENDIF
    997 
    998              ENDIF
    999 #if defined( __parallel )
    1000              CALL MPI_BARRIER( comm2d, ierr )
    1001 #endif
    1002           ENDDO
    1003 
    1004 !
    1005 !--       Calculate the index height of the topography
    1006           nzb_local = 0
    1007           DO  i = nxl, nxr
    1008              DO  j = nys, nyn
    1009                 IF ( .NOT. ocean )  THEN
    1010                    nzb_local(j,i) = MINLOC( ABS( zw - topo_height(j,i)         &
    1011                                                     - oro_height(j,i) ), 1 ) - 1
    1012                    IF ( ABS( zw(nzb_local(j,i)  ) - topo_height(j,i)           &
    1013                                                   - oro_height(j,i)  ) ==      &
    1014                         ABS( zw(nzb_local(j,i)+1) - topo_height(j,i)           &
    1015                                                   - oro_height(j,i)  )    )    &
    1016                       nzb_local(j,i) = nzb_local(j,i) + 1
    1017                 ELSE
    1018                    nzb_local(j,i) = MINLOC( ABS( zw - zw(0)                    &
    1019                                                     - topo_height(j,i)         &
    1020                                                     - oro_height(j,i) ), 1 ) - 1
    1021                    IF ( ABS( zw(nzb_local(j,i)  ) - zw(0)                      &
    1022                                                   - topo_height(j,i)           &
    1023                                                   - oro_height(j,i)  )  ==     &
    1024                         ABS( zw(nzb_local(j,i)+1) - zw(0)                      &
    1025                                                   - topo_height(j,i)           &
    1026                                                   - oro_height(j,i)  )    )    &
    1027                       nzb_local(j,i) = nzb_local(j,i) + 1
    1028                 ENDIF
    1029 
    1030              ENDDO
    1031           ENDDO
    1032 
    1033           DEALLOCATE ( oro_height  )
    1034           DEALLOCATE ( topo_height )
    1035 !
    1036 !--       Filter topography, i.e. fill holes resolved by only one grid point. 
    1037 !--       Such holes are suspected to lead to velocity blow-ups as continuity
    1038 !--       equation on discrete grid cannot be fulfilled in such case.
    1039 !--       For now, check only for holes and fill them to the lowest height level
    1040 !--       of the directly adjoining grid points along x- and y- direction.
    1041 !--       Before checking for holes, set lateral boundary conditions for
    1042 !--       topography. After hole-filling, boundary conditions must be set again!
    1043           CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    1044          
     1800!
     1801!--       Note, topography information have been already read. 
     1802!--       If required, further process topography, i.e. reference buildings on
     1803!--       top of orography and set temporary 3D topography array, which is
     1804!--       used later to set grid flags. Calling of this rouinte is also
     1805!--       required in case of ASCII input, even though no distinction between
     1806!--       terrain- and building height is made in this case. 
     1807          CALL process_topography( topo )
     1808!
     1809!--       Filter holes resolved by only one grid-point
     1810          CALL filter_topography( topo )
     1811!
     1812!--       Exchange ghost-points, as well as add cyclic or Neumann boundary
     1813!--       conditions.
     1814          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
     1815!
     1816!--       Set lateral boundary conditions for topography on all ghost layers         
    10451817          IF ( .NOT. bc_ns_cyc )  THEN
    10461818             IF ( nys == 0  )  THEN
    1047                 nzb_local(-1,:)   = nzb_local(0,:)
    1048                 nzb_local(-2,:)   = nzb_local(0,:)
    1049                 nzb_local(-3,:)   = nzb_local(0,:)
    1050              ELSEIF ( nyn == ny )  THEN
    1051                 nzb_local(ny+1,:) = nzb_local(ny,:)
    1052                 nzb_local(ny+2,:) = nzb_local(ny,:)
    1053                 nzb_local(ny+3,:) = nzb_local(ny,:)
     1819                DO  i = 1, nbgp         
     1820                   topo(:,nys-i,:) = topo(:,nys,:)
     1821                ENDDO
     1822             ENDIF
     1823             IF ( nyn == ny )  THEN
     1824                DO  i = 1, nbgp         
     1825                   topo(:,nyn+i,:) = topo(:,nyn,:)
     1826                ENDDO
    10541827             ENDIF
    10551828          ENDIF
     
    10571830          IF ( .NOT. bc_lr_cyc )  THEN
    10581831             IF ( nxl == 0  )  THEN
    1059                 nzb_local(:,-1)   = nzb_local(:,0)
    1060                 nzb_local(:,-2)   = nzb_local(:,0)
    1061                 nzb_local(:,-3)   = nzb_local(:,0)
    1062              ELSEIF ( nxr == nx )  THEN
    1063                 nzb_local(:,nx+1) = nzb_local(:,nx)
    1064                 nzb_local(:,nx+2) = nzb_local(:,nx)         
    1065                 nzb_local(:,nx+3) = nzb_local(:,nx)         
    1066              ENDIF         
     1832                DO  i = 1, nbgp
     1833                   topo(:,:,nxl-i) = topo(:,:,nxl)
     1834                ENDDO
     1835             ENDIF
     1836             IF ( nxr == nx )  THEN
     1837                DO  i = 1, nbgp
     1838                   topo(:,:,nxr+i) = topo(:,:,nxr)
     1839                ENDDO
     1840             ENDIF
    10671841          ENDIF
    10681842
    1069           num_hole_l = 0
    1070           DO i = nxl, nxr
    1071              DO j = nys, nyn
    1072 
    1073                 num_wall = 0
    1074 
    1075                 IF ( nzb_local(j-1,i) > nzb_local(j,i) )                       &
    1076                    num_wall = num_wall + 1
    1077                 IF ( nzb_local(j+1,i) > nzb_local(j,i) )                       &
    1078                    num_wall = num_wall + 1
    1079                 IF ( nzb_local(j,i-1) > nzb_local(j,i) )                       &
    1080                    num_wall = num_wall + 1
    1081                 IF ( nzb_local(j,i+1) > nzb_local(j,i) )                       &
    1082                    num_wall = num_wall + 1
    1083 
    1084                 IF ( num_wall == 4 )  THEN
    1085                    nzb_local(j,i) = MIN( nzb_local(j-1,i), nzb_local(j+1,i),   &
    1086                                          nzb_local(j,i-1), nzb_local(j,i+1) )
    1087                    num_hole_l     = num_hole_l + 1
    1088                 ENDIF
    1089              ENDDO
    1090           ENDDO
    1091 !
    1092 !--       Count the total number of holes, required for informative message.
    1093 #if defined( __parallel )
    1094           CALL MPI_ALLREDUCE( num_hole_l, num_hole, 1, MPI_INTEGER, MPI_SUM,   &
    1095                               comm2d, ierr )
    1096 #else
    1097           num_hole = num_hole_l
    1098 #endif   
    1099 !
    1100 !--       Create an informative message if any hole was removed.
    1101           IF ( num_hole > 0 )  THEN
    1102              WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//&
    1103                                                   'one grid point were filled'
    1104              CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 )
    1105           ENDIF
    1106 
    1107 !
    1108 !--       Set bit array to mask topography. Only required for lod = 1
    1109           IF ( lod == 1 )  THEN
    1110              DO  i = nxlg, nxrg
    1111                 DO  j = nysg, nyng
    1112                    topo_3d(nzb_local(j,i)+1:nzt+1,j,i) =                          &
    1113                                  IBSET( topo_3d(nzb_local(j,i)+1:nzt+1,j,i), 0 )
    1114                 ENDDO
    1115              ENDDO
    1116           ENDIF
    1117 !
    1118 !--       Exchange ghost-points, as well as add cyclic or Neumann boundary
    1119 !--       conditions.
    1120           CALL exchange_horiz_int( topo_3d, nbgp )
    1121           CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    1122          
    1123           IF ( .NOT. bc_ns_cyc )  THEN
    1124              IF ( nys == 0  )  topo_3d(:,-1,:)   = topo_3d(:,0,:)
    1125              IF ( nyn == ny )  topo_3d(:,ny+1,:) = topo_3d(:,ny,:)
    1126 
    1127              IF ( nys == 0  )  nzb_local(-1,:)   = nzb_local(0,:)
    1128              IF ( nyn == ny )  nzb_local(ny+1,:) = nzb_local(ny,:)
    1129           ENDIF
    1130 
    1131           IF ( .NOT. bc_lr_cyc )  THEN
    1132              IF ( nxl == 0  )  topo_3d(:,:,-1)   = topo_3d(:,:,0)
    1133              IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)     
    1134 
    1135              IF ( nxl == 0  )  nzb_local(:,-1)   = nzb_local(:,0)
    1136              IF ( nxr == nx )  nzb_local(:,nx+1) = nzb_local(:,nx)         
    1137           ENDIF
    1138 
    11391843
    11401844       CASE DEFAULT
    1141 !
     1845!   
    11421846!--       The DEFAULT case is reached either if the parameter topography
    11431847!--       contains a wrong character string or if the user has defined a special
    11441848!--       case in the user interface. There, the subroutine user_init_grid
    11451849!--       checks which of these two conditions applies.
    1146           CALL user_init_grid( topo_3d )
     1850          CALL user_init_grid( topo )
     1851          CALL filter_topography( topo )
    11471852
    11481853    END SELECT
    11491854!
    11501855!-- Consistency checks and index array initialization are only required for
    1151 !-- non-flat topography, also the initialization of topography height arrays
    1152 !-- zu_s_inner and zw_w_inner
     1856!-- non-flat topography.
    11531857    IF ( TRIM( topography ) /= 'flat' )  THEN
    1154 #if defined( __parallel )
    1155        CALL MPI_ALLREDUCE( MAXVAL( nzb_local ), nzb_local_max, 1, MPI_INTEGER, &
    1156                            MPI_MAX, comm2d, ierr )
    1157        CALL MPI_ALLREDUCE( MINVAL( nzb_local ), nzb_local_min, 1, MPI_INTEGER, &
    1158                            MPI_MIN, comm2d, ierr )                           
    1159 #else
    1160        nzb_local_max = MAXVAL( nzb_local )
    1161        nzb_local_min = MINVAL( nzb_local )
    1162 #endif
    1163 
    1164 !
    1165 !--    Consistency checks
    1166        IF ( nzb_local_min < 0  .OR.  nzb_local_max  > nz + 1 )  THEN
    1167           WRITE( message_string, * ) 'nzb_local values are outside the',       &
    1168                                 'model domain',                                &
    1169                                 '&MINVAL( nzb_local ) = ', nzb_local_min,      &
    1170                                 '&MAXVAL( nzb_local ) = ', nzb_local_max
    1171           CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
    1172        ENDIF
    11731858!
    11741859!--    In case of non-flat topography, check whether the convention how to
     
    11841869!--          defined in init_grid.
    11851870             WRITE( message_string, * )                                        &
    1186                   'The value for "topography_grid_convention" ',               &
    1187                   'is not set. Its default value is & only valid for ',        &
    1188                   '"topography" = ''single_building'', ',                      &
    1189                   '''single_street_canyon'' & or ''read_from_file''.',         &
    1190                   ' & Choose ''cell_edge'' or ''cell_center''.'
     1871               'The value for "topography_grid_convention" ',                  &
     1872               'is not set. Its default value is & only valid for ',           &
     1873               '"topography" = ''single_building'', ',                         &
     1874               '''single_street_canyon'' & or ''read_from_file''.',            &
     1875               ' & Choose ''cell_edge'' or ''cell_center''.'
    11911876             CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 )
    11921877          ELSE
     
    12041889                TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
    12051890          WRITE( message_string, * )                                           &
    1206                'The value for "topography_grid_convention" is ',               &
    1207                'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
     1891            'The value for "topography_grid_convention" is ',                  &
     1892            'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
    12081893          CALL message( 'init_grid', 'PA0240', 1, 2, 0, 6, 0 )
    12091894       ENDIF
    12101895
    1211 !
    1212 !--    In case of non-flat topography, check whether the convention how to
    1213 !--    define the topography grid has been set correctly, or whether the default
    1214 !--    is applicable. If this is not possible, abort.
    1215        IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
    1216           IF ( TRIM( topography ) /= 'single_building' .AND.                   &
    1217                TRIM( topography ) /= 'single_street_canyon' .AND.              &
    1218                TRIM( topography ) /= 'read_from_file' )  THEN
    1219 !--          The default value is not applicable here, because it is only valid
    1220 !--          for the two standard cases 'single_building' and 'read_from_file'
    1221 !--          defined in init_grid.
    1222              WRITE( message_string, * )                                        &
    1223                   'The value for "topography_grid_convention" ',               &
    1224                   'is not set. Its default value is & only valid for ',        &
    1225                   '"topography" = ''single_building'', ',                      &
    1226                   '''single_street_canyon'' & or ''read_from_file''.',         &
    1227                   ' & Choose ''cell_edge'' or ''cell_center''.'
    1228              CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 )
    1229           ELSE
    1230 !--          The default value is applicable here.
    1231 !--          Set convention according to topography.
    1232              IF ( TRIM( topography ) == 'single_building' .OR.                 &
    1233                   TRIM( topography ) == 'single_street_canyon' )  THEN
    1234                 topography_grid_convention = 'cell_edge'
    1235              ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
    1236                 topography_grid_convention = 'cell_center'
    1237              ENDIF
    1238           ENDIF
    1239        ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.        &
    1240                 TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
    1241           WRITE( message_string, * )                                           &
    1242                'The value for "topography_grid_convention" is ',               &
    1243                'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
    1244           CALL message( 'init_grid', 'PA0240', 1, 2, 0, 6, 0 )
    1245        ENDIF
    12461896
    12471897       IF ( topography_grid_convention == 'cell_edge' )  THEN
     
    12591909          DO  j = nys+1, nyn+1
    12601910             DO  i = nxl-1, nxr
    1261                 nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) )
    1262              ENDDO
    1263           ENDDO
    1264 !
    1265 !--       Exchange ghost points
    1266           CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    1267 
    1268           DO  i = nxl, nxr+1
    1269              DO  j = nys-1, nyn
    1270                 nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) )
    1271              ENDDO
    1272           ENDDO 
    1273 !
    1274 !--       Exchange ghost points         
    1275           CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    1276 !
    1277 !--       Apply cell-edge convention also for 3D topo array. The former setting
    1278 !--       of nzb_local will be removed later.
    1279           DO  j = nys+1, nyn+1
    1280              DO  i = nxl-1, nxr
    12811911                DO  k = nzb, nzt+1
    1282                    IF ( BTEST( topo_3d(k,j,i), 0 )  .OR.                       &
    1283                         BTEST( topo_3d(k,j,i+1), 0 ) )                         &
    1284                       topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 )
     1912                   IF ( BTEST( topo(k,j,i), 0 )  .OR.                          &
     1913                        BTEST( topo(k,j,i+1), 0 ) )                            &
     1914                       topo(k,j,i) = IBSET( topo(k,j,i), 0 )
    12851915                ENDDO
    12861916             ENDDO
    12871917          ENDDO     
    1288           CALL exchange_horiz_int( topo_3d, nbgp )   
     1918          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
    12891919
    12901920          DO  i = nxl, nxr+1
    12911921             DO  j = nys-1, nyn
    12921922                DO  k = nzb, nzt+1
    1293                    IF ( BTEST( topo_3d(k,j,i), 0 )  .OR.                       &
    1294                         BTEST( topo_3d(k,j+1,i), 0 ) )                         &
    1295                       topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 )
     1923                   IF ( BTEST( topo(k,j,i), 0 )  .OR.                          &
     1924                        BTEST( topo(k,j+1,i), 0 ) )                            &
     1925                      topo(k,j,i) = IBSET( topo(k,j,i), 0 )
    12961926                ENDDO
    12971927             ENDDO
    12981928          ENDDO 
    1299           CALL exchange_horiz_int( topo_3d, nbgp )
     1929          CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
    13001930   
    13011931       ENDIF
    1302      
    1303 !
    1304 !--    Initialize index arrays nzb_s_inner and nzb_w_inner
    1305 
    1306        nzb_s_inner = nzb_local
    1307        nzb_w_inner = nzb_local
    1308 
    1309 !
    1310 !--    Initialize remaining index arrays:
    1311 !--    first pre-initialize them with nzb_s_inner...
    1312        nzb_u_inner = nzb_s_inner
    1313        nzb_u_outer = nzb_s_inner
    1314        nzb_v_inner = nzb_s_inner
    1315        nzb_v_outer = nzb_s_inner
    1316        nzb_w_outer = nzb_s_inner
    1317        nzb_s_outer = nzb_s_inner
    1318 
    1319 !
    1320 !--    ...then extend pre-initialized arrays in their according directions
    1321 !--    based on nzb_local using nzb_tmp as a temporary global index array
    1322 
    1323 !
    1324 !--    nzb_s_outer:
    1325 !--    extend nzb_local east-/westwards first, then north-/southwards
    1326        nzb_tmp = nzb_local
    1327        DO  j = nys, nyn
    1328           DO  i = nxl, nxr
    1329              nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i),             &
    1330                                  nzb_local(j,i+1) )
    1331           ENDDO
    1332        ENDDO
    1333        
    1334        CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
    1335        
    1336        DO  i = nxl, nxr
    1337           DO  j = nys, nyn
    1338              nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
    1339                                      nzb_tmp(j+1,i) )
    1340           ENDDO
    1341 !
    1342 !--       non-cyclic boundary conditions (overwritten by call of
    1343 !--       exchange_horiz_2d_int below in case of cyclic boundary conditions)
    1344           IF ( nys == 0 )  THEN
    1345              j = -1
    1346              nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
    1347           ENDIF
    1348           IF ( nyn == ny )  THEN
    1349              j = ny + 1
    1350              nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
    1351           ENDIF
    1352        ENDDO
    1353 !
    1354 !--    nzb_w_outer:
    1355 !--    identical to nzb_s_outer
    1356        nzb_w_outer = nzb_s_outer
    1357 
    1358 !
    1359 !--    nzb_u_inner:
    1360 !--    extend nzb_local rightwards only
    1361        nzb_tmp = nzb_local
    1362        DO  j = nys, nyn
    1363           DO  i = nxl, nxr
    1364              nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
    1365           ENDDO
    1366        ENDDO
    1367        
    1368        CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
    1369        
    1370        nzb_u_inner = nzb_tmp
    1371 !
    1372 !--    nzb_u_outer:
    1373 !--    extend current nzb_tmp (nzb_u_inner) north-/southwards
    1374        DO  i = nxl, nxr
    1375           DO  j = nys, nyn
    1376              nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
    1377                                      nzb_tmp(j+1,i) )
    1378           ENDDO
    1379 !
    1380 !--       non-cyclic boundary conditions (overwritten by call of
    1381 !--       exchange_horiz_2d_int below in case of cyclic boundary conditions)
    1382           IF ( nys == 0 )  THEN
    1383              j = -1
    1384              nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
    1385           ENDIF
    1386           IF ( nyn == ny )  THEN
    1387              j = ny + 1
    1388              nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
    1389           ENDIF
    1390        ENDDO
    1391 
    1392 !
    1393 !--    nzb_v_inner:
    1394 !--    extend nzb_local northwards only
    1395        nzb_tmp = nzb_local
    1396        DO  i = nxl, nxr
    1397           DO  j = nys, nyn
    1398              nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
    1399           ENDDO
    1400        ENDDO
    1401        
    1402        CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )     
    1403        nzb_v_inner = nzb_tmp
    1404 
    1405 !
    1406 !--    nzb_v_outer:
    1407 !--    extend current nzb_tmp (nzb_v_inner) right-/leftwards
    1408        DO  j = nys, nyn
    1409           DO  i = nxl, nxr
    1410              nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i),             &
    1411                                      nzb_tmp(j,i+1) )
    1412           ENDDO
    1413 !
    1414 !--       non-cyclic boundary conditions (overwritten by call of
    1415 !--       exchange_horiz_2d_int below in case of cyclic boundary conditions)
    1416           IF ( nxl == 0 )  THEN
    1417              i = -1
    1418              nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
    1419           ENDIF
    1420           IF ( nxr == nx )  THEN
    1421              i = nx + 1
    1422              nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
    1423           ENDIF
    1424        ENDDO
    1425 
    1426 !
    1427 !--    Exchange of lateral boundary values (parallel computers) and cyclic
    1428 !--    boundary conditions, if applicable.
    1429 !--    Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
    1430 !--    they do not require exchange and are not included here.
    1431        CALL exchange_horiz_2d_int( nzb_u_inner, nys, nyn, nxl, nxr, nbgp )
    1432        CALL exchange_horiz_2d_int( nzb_u_outer, nys, nyn, nxl, nxr, nbgp )
    1433        CALL exchange_horiz_2d_int( nzb_v_inner, nys, nyn, nxl, nxr, nbgp )
    1434        CALL exchange_horiz_2d_int( nzb_v_outer, nys, nyn, nxl, nxr, nbgp )
    1435        CALL exchange_horiz_2d_int( nzb_w_outer, nys, nyn, nxl, nxr, nbgp )
    1436        CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp )
    1437 
    14381932    ENDIF
    1439 !
    1440 !-- Deallocate temporary array, as it might be reused for different
    1441 !-- grid-levels further below.
    1442     DEALLOCATE( nzb_tmp )
    1443 
    1444 !
    1445 !-- Determine the maximum level of topography. It is used for
    1446 !-- steering the degradation of order of the applied advection scheme.
    1447 !-- In case of non-cyclic lateral boundaries, the order of the advection
    1448 !-- scheme has to be reduced up to nzt (required at the lateral boundaries).
    1449     k_top = 0
     1933
     1934
     1935 END SUBROUTINE init_topo
     1936
     1937 SUBROUTINE set_topo_flags(topo)
     1938
     1939    USE control_parameters,                                                    &
     1940        ONLY:  bc_lr_cyc, bc_ns_cyc, constant_flux_layer, land_surface,        &
     1941               use_surface_fluxes, use_top_fluxes, urban_surface
     1942
     1943    USE indices,                                                               &
     1944        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
     1945               nzb, nzt, wall_flags_0
     1946
     1947    USE kinds
     1948
     1949    IMPLICIT NONE
     1950
     1951    INTEGER(iwp) ::  i             !< index variable along x
     1952    INTEGER(iwp) ::  j             !< index variable along y
     1953    INTEGER(iwp) ::  k             !< index variable along z
     1954
     1955    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
     1956
     1957    ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1958    wall_flags_0 = 0
     1959!
     1960!-- Set-up topography flags. First, set flags only for s, u, v and w-grid.
     1961!-- Further special flags will be set in following loops.
    14501962    DO  i = nxl, nxr
    14511963       DO  j = nys, nyn
    1452           DO  k = nzb, nzt + 1
    1453              k_top = MAX( k_top, MERGE( k, 0,                                  &
    1454                                         .NOT. BTEST( topo_3d(k,j,i), 0 ) ) )
     1964          DO  k = nzb, nzt+1
     1965!
     1966!--          scalar grid
     1967             IF ( BTEST( topo(k,j,i), 0 ) )                                 &
     1968                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
     1969!
     1970!--          u grid
     1971             IF ( BTEST( topo(k,j,i),   0 )  .AND.                          &
     1972                  BTEST( topo(k,j,i-1), 0 ) )                               &
     1973                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
     1974!
     1975!--          v grid
     1976             IF ( BTEST( topo(k,j,i),   0 )  .AND.                          &
     1977                  BTEST( topo(k,j-1,i), 0 ) )                               &
     1978                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
     1979
    14551980          ENDDO
    1456        ENDDO
    1457     ENDDO
    1458 #if defined( __parallel )
    1459     CALL MPI_ALLREDUCE( k_top + 1, nzb_max, 1, MPI_INTEGER,                    & !is +1 really necessary here?
    1460                         MPI_MAX, comm2d, ierr )
    1461 #else
    1462     nzb_max = k_top + 1
    1463 #endif
    1464     IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR.             &
    1465          inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s .OR.             &
    1466          nest_domain )                                                         &
    1467     THEN
    1468        nzb_max = nzt
    1469     ENDIF
    1470 !
    1471 !-- Finally, if topography extents up to the model top, limit nzb_max to nzt.
    1472     nzb_max = MIN( nzb_max, nzt )
    1473 
    1474 !
    1475 !-- Set the individual index arrays which define the k index from which on
    1476 !-- the usual finite difference form (which does not use surface fluxes) is
    1477 !-- applied
    1478     IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
    1479        nzb_diff_s_inner   = nzb_s_inner + 2
    1480        nzb_diff_s_outer   = nzb_s_outer + 2
    1481     ELSE
    1482        nzb_diff_s_inner   = nzb_s_inner + 1
    1483        nzb_diff_s_outer   = nzb_s_outer + 1
    1484     ENDIF
    1485 
    1486 !
    1487 !-- Set-up topography flags. First, set flags only for s, u, v and w-grid.
    1488 !-- Further special flags will be set in following loops.
    1489     wall_flags_0 = 0
    1490     DO  j = nys, nyn
    1491        DO  i = nxl, nxr
    1492           DO  k = nzb, nzt+1
    1493 !
    1494 !--          scalar grid
    1495              IF ( BTEST( topo_3d(k,j,i), 0 ) )                                 &
    1496                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
    1497 !
    1498 !--          v grid
    1499              IF ( BTEST( topo_3d(k,j,i),   0 )  .AND.                          &
    1500                   BTEST( topo_3d(k,j-1,i), 0 ) )                               &
    1501                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
    1502 !
    1503 !--     To do: set outer arrays on basis of topo_3d array, adjust for downward-facing walls
    1504 !--          s grid outer array
    1505              IF ( k >= nzb_s_outer(j,i) )                                      &
    1506                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
    1507 !
    1508 !--          s grid outer array
    1509              IF ( k >= nzb_u_outer(j,i) )                                      &
    1510                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 )
    1511 !
    1512 !--          s grid outer array
    1513              IF ( k >= nzb_v_outer(j,i) )                                      &
    1514                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
    1515 !
    1516 !--          w grid outer array
    1517              IF ( k >= nzb_w_outer(j,i) )                                      &
    1518                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
    1519           ENDDO
    15201981
    15211982          DO k = nzb, nzt
    15221983!
    15231984!--          w grid
    1524              IF ( BTEST( topo_3d(k,j,i),   0 )  .AND.                          &
    1525                   BTEST( topo_3d(k+1,j,i), 0 ) )                               &
     1985             IF ( BTEST( topo(k,j,i),   0 )  .AND.                          &
     1986                  BTEST( topo(k+1,j,i), 0 ) )                               &
    15261987                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
    15271988          ENDDO
     
    15301991       ENDDO
    15311992    ENDDO
    1532 !
    1533 !-- u grid. Note, reverse
    1534 !-- memory access is required for setting flag on u-grid
    1535     DO  j = nys, nyn
    1536        DO  i = nxl, nxr
     1993
     1994    CALL exchange_horiz_int ( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
     1995!
     1996!-- Set outer array for scalars to mask near-surface grid points in
     1997!-- production_e
     1998    DO i = nxl, nxr
     1999       DO j = nys, nyn
    15372000          DO k = nzb, nzt+1
    1538              IF ( BTEST( topo_3d(k,j,i),   0 )  .AND.                          &
    1539                   BTEST( topo_3d(k,j,i-1), 0 ) )                               &
    1540                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
     2001             IF ( BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.                       &
     2002                  BTEST( wall_flags_0(k,j+1,i), 0 )  .AND.                       &
     2003                  BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.                       &
     2004                  BTEST( wall_flags_0(k,j-1,i-1), 0 )  .AND.                       &
     2005                  BTEST( wall_flags_0(k,j+1,i-1), 0 )  .AND.                       &
     2006                  BTEST( wall_flags_0(k,j-1,i+1), 0 )  .AND.                       &
     2007                  BTEST( wall_flags_0(k,j+1,i+1), 0 ) )                            &
     2008                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
    15412009          ENDDO
    15422010       ENDDO
     
    15762044             wall_flags_0(nzt+1,j,i) = IBCLR( wall_flags_0(nzt+1,j,i), 9 )
    15772045
     2046
    15782047          DO k = nzb+1, nzt
    15792048!
     
    16342103             IF ( BTEST( wall_flags_0(k-1,j,i), 0 )  .AND.                     &
    16352104            .NOT. BTEST( wall_flags_0(k,j,i), 0   ) )                          &
    1636                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
     2105                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
    16372106!
    16382107!--          Downward facing wall on u grid
     
    16652134                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
    16662135
    1667 !
     2136!   
    16682137!--          Upward facing wall on v grid
    16692138             IF ( .NOT. BTEST( wall_flags_0(k,j,i),   2 )  .AND.               &
    16702139                        BTEST( wall_flags_0(k+1,j,i), 2 ) )                    &
    16712140                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
    1672 
     2141   
    16732142!
    16742143!--          Upward facing wall on w grid
     
    16812150                  BTEST( wall_flags_0(k,j,i), 12 ) .OR.                        &
    16822151                  BTEST( wall_flags_0(k,j,i), 13 ) )                           &
    1683                 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
     2152                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
    16842153!
    16852154!--          Special flag on scalar grid, nzb_diff_s_inner - 1, required for
     
    16882157                IF ( BTEST( wall_flags_0(k,j,i),   0 )  .AND.                  &
    16892158                     BTEST( wall_flags_0(k+1,j,i), 0 ) )                       &
    1690                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
     2159                  wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
    16912160             ELSE
    16922161                IF ( BTEST( wall_flags_0(k,j,i), 22 ) )                        &
    16932162                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
    16942163             ENDIF
    1695 
     2164   
    16962165
    16972166          ENDDO
     
    17012170    ENDDO
    17022171!
     2172!-- Finally, set identification flags indicating natural terrain or buildings.
     2173!-- Natural terrain grid points.
     2174    IF ( land_surface )  THEN
     2175       DO i = nxl, nxr
     2176          DO j = nys, nyn
     2177             DO k = nzb, nzt+1
     2178!
     2179!--             Natural terrain grid point
     2180                IF ( BTEST( topo(k,j,i), 1 ) )                                 &
     2181                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
     2182             ENDDO
     2183          ENDDO
     2184       ENDDO
     2185    ENDIF
     2186!
     2187!-- Building grid points.
     2188    IF ( urban_surface )  THEN
     2189       DO i = nxl, nxr
     2190          DO j = nys, nyn
     2191             DO k = nzb, nzt+1
     2192                IF ( BTEST( topo(k,j,i), 2 ) )                                 &
     2193                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
     2194             ENDDO
     2195          ENDDO
     2196       ENDDO
     2197    ENDIF
     2198
     2199!
    17032200!-- Exchange ghost points for wall flags
    1704     CALL exchange_horiz_int( wall_flags_0, nbgp )
     2201    CALL exchange_horiz_int( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
    17052202!
    17062203!-- Set boundary conditions also for flags. Can be interpreted as Neumann
    17072204!-- boundary conditions for topography.
    17082205    IF ( .NOT. bc_ns_cyc )  THEN
    1709        IF ( nys == 0  )  wall_flags_0(:,-1,:)   = wall_flags_0(:,0,:)
    1710        IF ( nyn == ny )  wall_flags_0(:,ny+1,:) = wall_flags_0(:,ny,:)
     2206       IF ( nys == 0  )  THEN
     2207          DO  i = 1, nbgp     
     2208             wall_flags_0(:,nys-i,:)   = wall_flags_0(:,nys,:)
     2209          ENDDO
     2210       ENDIF
     2211       IF ( nyn == ny )  THEN
     2212          DO  i = 1, nbgp 
     2213             wall_flags_0(:,nyn+i,:) = wall_flags_0(:,nyn,:)
     2214          ENDDO
     2215       ENDIF
    17112216    ENDIF
    17122217    IF ( .NOT. bc_lr_cyc )  THEN
    1713        IF ( nxl == 0  )  wall_flags_0(:,:,-1)   = wall_flags_0(:,:,0)
    1714        IF ( nxr == nx )  wall_flags_0(:,:,nx+1) = wall_flags_0(:,:,nx)           
     2218       IF ( nxl == 0  )  THEN
     2219          DO  i = 1, nbgp   
     2220             wall_flags_0(:,:,nxl-i)   = wall_flags_0(:,:,nxl)
     2221          ENDDO
     2222       ENDIF
     2223       IF ( nxr == nx )  THEN
     2224          DO  i = 1, nbgp   
     2225             wall_flags_0(:,:,nxr+i) = wall_flags_0(:,:,nxr)     
     2226          ENDDO
     2227       ENDIF     
    17152228    ENDIF
    17162229
    1717 !
    1718 !-- Initialize boundary conditions via surface type
    1719     CALL init_bc
    1720 !
    1721 !-- Allocate and set topography height arrays required for data output
    1722     IF ( TRIM( topography ) /= 'flat' )  THEN
    1723 !
    1724 !--    Allocate and set the arrays containing the topography height
    1725 
    1726        IF ( nxr == nx  .AND.  nyn /= ny )  THEN
    1727           ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
    1728                     zw_w_inner(nxl:nxr+1,nys:nyn) )
    1729        ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
    1730           ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
    1731                     zw_w_inner(nxl:nxr,nys:nyn+1) )
    1732        ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
    1733           ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
    1734                     zw_w_inner(nxl:nxr+1,nys:nyn+1) )
    1735        ELSE
    1736           ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
    1737                     zw_w_inner(nxl:nxr,nys:nyn) )
    1738        ENDIF
    1739 
    1740        zu_s_inner   = 0.0_wp
    1741        zw_w_inner   = 0.0_wp
    1742 !
    1743 !--    Determine local topography height on scalar and w-grid. Note, setting
    1744 !--    lateral boundary values is not necessary, realized via wall_flags_0
    1745 !--    array. Further, please note that loop bounds are different from
    1746 !--    nxl to nxr and nys to nyn on south and right model boundary, hence,
    1747 !--    use intrinsic lbound and ubound functions to infer array bounds.
    1748        DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
    1749           DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
    1750 !
    1751 !--          Topography height on scalar grid. Therefore, determine index of
    1752 !--          upward-facing surface element on scalar grid.
    1753              zu_s_inner(i,j) = zu( get_topography_top_index( j, i, 's' ) )
    1754 !
    1755 !--          Topography height on w grid. Therefore, determine index of
    1756 !--          upward-facing surface element on w grid.
    1757              zw_w_inner(i,j) = zw( get_topography_top_index( j, i, 's' ) )
    1758           ENDDO
    1759        ENDDO
    1760 
    1761     ENDIF
    1762 
    1763 !
    1764 !-- Calculate wall flag arrays for the multigrid method.
    1765 !-- Please note, wall flags are only applied in the not cache-optimized
    1766 !-- version.
    1767     IF ( psolver == 'multigrid_noopt' )  THEN
    1768 
    1769 !
    1770 !--    Gridpoint increment of the current level.
    1771        inc = 1
    1772        DO  l = maximum_grid_level, 1 , -1
    1773 !
    1774 !--       Set grid_level as it is required for exchange_horiz_2d_int
    1775           grid_level = l
    1776 
    1777           nxl_l = nxl_mg(l)
    1778           nxr_l = nxr_mg(l)
    1779           nys_l = nys_mg(l)
    1780           nyn_l = nyn_mg(l)
    1781           nzt_l = nzt_mg(l)
    1782 !
    1783 !--       Assign the flag level to be calculated
    1784           SELECT CASE ( l )
    1785              CASE ( 1 )
    1786                 flags => wall_flags_1
    1787              CASE ( 2 )
    1788                 flags => wall_flags_2
    1789              CASE ( 3 )
    1790                 flags => wall_flags_3
    1791              CASE ( 4 )
    1792                 flags => wall_flags_4
    1793              CASE ( 5 )
    1794                 flags => wall_flags_5
    1795              CASE ( 6 )
    1796                 flags => wall_flags_6
    1797              CASE ( 7 )
    1798                 flags => wall_flags_7
    1799              CASE ( 8 )
    1800                 flags => wall_flags_8
    1801              CASE ( 9 )
    1802                 flags => wall_flags_9
    1803              CASE ( 10 )
    1804                 flags => wall_flags_10
    1805           END SELECT
    1806 
    1807 !
    1808 !--       Depending on the grid level, set the respective bits in case of
    1809 !--       neighbouring walls
    1810 !--       Bit 0:  wall to the bottom
    1811 !--       Bit 1:  wall to the top (not realized in remaining PALM code so far)
    1812 !--       Bit 2:  wall to the south
    1813 !--       Bit 3:  wall to the north
    1814 !--       Bit 4:  wall to the left
    1815 !--       Bit 5:  wall to the right
    1816 !--       Bit 6:  inside building
    1817 
    1818           flags = 0
    1819 
    1820 !
    1821 !--       In case of masking method, flags are not set and multigrid method
    1822 !--       works like FFT-solver
    1823           IF ( .NOT. masking_method )  THEN
    1824 
    1825 !
    1826 !--          Allocate temporary array for topography heights on coarser grid
    1827 !--          level. Please note, 2 ghoist points are required, in order to
    1828 !--          calculate flags() on the interior ghost point.
    1829              ALLOCATE( nzb_tmp(nys_l-2:nyn_l+2,nxl_l-2:nxr_l+2) )
    1830              nzb_tmp = 0
    1831              
    1832              DO  i = nxl_l, nxr_l
    1833                 DO  j = nys_l, nyn_l
    1834                    nzb_tmp(j,i) = nzb_local(j*inc,i*inc)
    1835                 ENDDO
    1836              ENDDO
    1837 !
    1838 !--          Exchange ghost points on respective multigrid level. 2 ghost points
    1839 !--          are required, in order to calculate flags on
    1840 !--          nys_l-1 / nyn_l+1 / nxl_l-1 / nxr_l+1. The alternative would be to
    1841 !--          exchange 3D-INTEGER array flags on the respective multigrid level.
    1842              CALL exchange_horiz_2d_int( nzb_tmp, nys_l, nyn_l, nxl_l, nxr_l, 2 )
    1843 !
    1844 !--          Set non-cyclic boundary conditions on respective multigrid level
    1845              IF ( .NOT. bc_ns_cyc )  THEN
    1846                 IF ( inflow_s .OR. outflow_s .OR. nest_bound_s  )  THEN
    1847                    nzb_tmp(-2,:) = nzb_tmp(0,:)
    1848                    nzb_tmp(-1,:) = nzb_tmp(0,:)
    1849                 ENDIF
    1850                 IF ( inflow_n .OR. outflow_n .OR. nest_bound_n )  THEN
    1851                    nzb_tmp(nyn_l+2,:) = nzb_tmp(nyn_l,:)
    1852                    nzb_tmp(nyn_l+1,:) = nzb_tmp(nyn_l,:)
    1853                 ENDIF
    1854              ENDIF
    1855              IF ( .NOT. bc_lr_cyc )  THEN
    1856                 IF ( inflow_l .OR. outflow_l .OR. nest_bound_l  )  THEN
    1857                    nzb_tmp(:,-2) = nzb_tmp(:,0)
    1858                    nzb_tmp(:,-1) = nzb_tmp(:,0)
    1859                 ENDIF
    1860                 IF ( inflow_r .OR. outflow_r .OR. nest_bound_r )  THEN
    1861                    nzb_tmp(:,nxr_l+1) = nzb_tmp(:,nxr_l)   
    1862                    nzb_tmp(:,nxr_l+2) = nzb_tmp(:,nxr_l)     
    1863                 ENDIF       
    1864              ENDIF
    1865                        
    1866              DO  i = nxl_l-1, nxr_l+1
    1867                 DO  j = nys_l-1, nyn_l+1
    1868                    DO  k = nzb, nzt_l+1     
    1869 !
    1870 !--                   Inside/outside building (inside building does not need
    1871 !--                   further tests for walls)
    1872                       IF ( k*inc <= nzb_tmp(j,i) )  THEN
    1873 
    1874                          flags(k,j,i) = IBSET( flags(k,j,i), 6 )
    1875 
    1876                       ELSE
    1877 !
    1878 !--                      Bottom wall
    1879                          IF ( (k-1)*inc <= nzb_tmp(j,i) )  THEN
    1880                             flags(k,j,i) = IBSET( flags(k,j,i), 0 )
    1881                          ENDIF
    1882 !
    1883 !--                      South wall
    1884                          IF ( k*inc <= nzb_tmp(j-1,i) )  THEN
    1885                             flags(k,j,i) = IBSET( flags(k,j,i), 2 )
    1886                          ENDIF
    1887 !
    1888 !--                      North wall
    1889                          IF ( k*inc <= nzb_tmp(j+1,i) )  THEN
    1890                             flags(k,j,i) = IBSET( flags(k,j,i), 3 )
    1891                          ENDIF
    1892 !
    1893 !--                      Left wall
    1894                          IF ( k*inc <= nzb_tmp(j,i-1) )  THEN
    1895                             flags(k,j,i) = IBSET( flags(k,j,i), 4 )
    1896                          ENDIF
    1897 !
    1898 !--                      Right wall
    1899                          IF ( k*inc <= nzb_tmp(j,i+1) )  THEN
    1900                             flags(k,j,i) = IBSET( flags(k,j,i), 5 )
    1901                          ENDIF
    1902 
    1903                       ENDIF
    1904                            
    1905                    ENDDO
    1906                 ENDDO
    1907              ENDDO
    1908 
    1909              DEALLOCATE( nzb_tmp )
    1910 
    1911           ENDIF
    1912 
    1913           inc = inc * 2
    1914 
    1915        ENDDO
    1916 !
    1917 !--    Reset grid_level to "normal" grid
    1918        grid_level = 0
    1919        
    1920     ENDIF
    1921 !
    1922 !-- Allocate flags needed for masking walls. Even though these flags are only
    1923 !-- required in the ws-scheme, the arrays need to be allocated here as they are
    1924 !-- used in OpenACC directives.
    1925     ALLOCATE( advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
    1926               advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    1927     advc_flags_1 = 0
    1928     advc_flags_2 = 0
    1929 !
    1930 !-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e.
    1931 !-- to decrease the numerical stencil appropriately.
    1932     IF ( momentum_advec == 'ws-scheme'  .OR.  scalar_advec == 'ws-scheme'      &
    1933        )  THEN
    1934        CALL ws_init_flags
    1935     ENDIF
    1936 
    1937 !
    1938 !-- In case of topography: limit near-wall mixing length l_wall further:
    1939 !-- Go through all points of the subdomain one by one and look for the closest
    1940 !-- surface
    1941     DO  i = nxl, nxr
    1942        DO  j = nys, nyn
    1943           DO  k = nzb+1, nzt
    1944 !
    1945 !--          Check if current gridpoint belongs to the atmosphere
    1946              IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
    1947 !
    1948 !--             Check for neighbouring grid-points.
    1949 !--             Vertical distance, down
    1950                 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )                 &
    1951                    l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
    1952 !
    1953 !--             Vertical distance, up
    1954                 IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )                 &
    1955                    l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
    1956 !
    1957 !--             y-distance
    1958                 IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .OR.             &
    1959                      .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )                 &
    1960                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
    1961 !
    1962 !--             x-distance
    1963                 IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .OR.             &
    1964                      .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )                 &
    1965                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
    1966 !
    1967 !--              yz-distance (vertical edges, down)
    1968                  IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 )  .OR.          &
    1969                       .NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 )  )             &
    1970                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    1971                                         SQRT( 0.25_wp * dy**2 +                &
    1972                                        ( zu(k) - zw(k-1) )**2 ) )
    1973 !
    1974 !--              yz-distance (vertical edges, up)
    1975                  IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 )  .OR.          &
    1976                       .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 )  )             &
    1977                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    1978                                         SQRT( 0.25_wp * dy**2 +                &
    1979                                        ( zw(k) - zu(k) )**2 ) )
    1980 !
    1981 !--              xz-distance (vertical edges, down)
    1982                  IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 )  .OR.          &
    1983                       .NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 )  )             &
    1984                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    1985                                         SQRT( 0.25_wp * dx**2 +                &
    1986                                        ( zu(k) - zw(k-1) )**2 ) )
    1987 !
    1988 !--              xz-distance (vertical edges, up)
    1989                  IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 )  .OR.          &
    1990                       .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 )  )             &
    1991                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    1992                                         SQRT( 0.25_wp * dx**2 +                &
    1993                                        ( zw(k) - zu(k) )**2 ) )
    1994 !
    1995 !--             xy-distance (horizontal edges)
    1996                 IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 )  .OR.           &
    1997                      .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.           &
    1998                      .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.           &
    1999                      .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) )               &
    2000                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    2001                                         SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
    2002 !
    2003 !--             xyz distance (vertical and horizontal edges, down)
    2004                 IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 )  .OR.         &
    2005                      .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.         &
    2006                      .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.         &
    2007                      .NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) )             &
    2008                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    2009                                         SQRT( 0.25_wp * ( dx**2 + dy**2 )      &
    2010                                               +  ( zu(k) - zw(k-1) )**2  ) )
    2011 !
    2012 !--             xyz distance (vertical and horizontal edges, up)
    2013                 IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 )  .OR.         &
    2014                      .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.         &
    2015                      .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.         &
    2016                      .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) )             &
    2017                    l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),              &
    2018                                         SQRT( 0.25_wp * ( dx**2 + dy**2 )      &
    2019                                               +  ( zw(k) - zu(k) )**2  ) )
    2020                  
    2021              ENDIF
    2022           ENDDO
    2023        ENDDO
    2024     ENDDO
    2025 !
    2026 !-- Set lateral boundary conditions for l_wall
    2027     CALL exchange_horiz( l_wall, nbgp )     
    2028 
    2029 !
    2030 !-- Vertical nesting: communicate vertical grid level arrays between fine and
    2031 !-- coarse grid
    2032     IF ( vnested )  CALL vnest_init_grid
    2033 
    2034  END SUBROUTINE init_grid
     2230
     2231 END SUBROUTINE set_topo_flags
     2232
     2233
     2234
Note: See TracChangeset for help on using the changeset viewer.