Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 4413)
+++ /palm/trunk/SOURCE/Makefile	(revision 4414)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# Move dependencies for init grid from advection scheme and multigrid solver
+# to module_interface
+# 
+# 4411 2020-02-18 14:28:02Z maronga
 # Added output routines for WTM
 # 
@@ -664,11 +668,8 @@
 	vertical_nesting_mod.o
 init_grid.o: \
-	advec_ws.o \
-	mod_kinds.o \
-	modules.o \
-	netcdf_data_input_mod.o \
-	netcdf_interface_mod.o \
-	poismg_noopt_mod.o \
-	surface_mod.o \
+	mod_kinds.o \
+	modules.o \
+	netcdf_data_input_mod.o \
+	netcdf_interface_mod.o \
 	vertical_nesting_mod.o
 init_masks.o: \
@@ -769,8 +770,10 @@
 	ocean_mod.o \
 	plant_canopy_model_mod.o \
+	poismg_noopt_mod.o \
 	radiation_model_mod.o \
 	salsa_mod.o \
 	spectra_mod.o \
 	surface_data_output_mod.o \
+	surface_mod.o \
 	synthetic_turbulence_generator_mod.o \
 	turbulence_closure_mod.o \
Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 4413)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 4414)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Move call for initialization of control flags to ws_init
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
@@ -379,4 +382,20 @@
 
        ENDIF
+!
+!--    Initialize the flag arrays controlling degradation near walls, i.e. 
+!--    to decrease the numerical stencil appropriately. The order of the scheme
+!--    is degraded near solid walls as well as near non-cyclic inflow and outflow
+!--    boundaries. Do this separately for momentum and scalars.
+       IF ( ws_scheme_mom )  CALL ws_init_flags_momentum
+
+       IF ( ws_scheme_sca )  THEN
+          ALLOCATE( advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          advc_flags_s = 0
+          CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l,     &
+                                     bc_dirichlet_n  .OR.  bc_radiation_n,     &
+                                     bc_dirichlet_r  .OR.  bc_radiation_r,     &
+                                     bc_dirichlet_s  .OR.  bc_radiation_s,     &
+                                     advc_flags_s )
+       ENDIF
 
     END SUBROUTINE ws_init
@@ -398,9 +417,9 @@
        INTEGER(iwp) ::  k_pp  !< dummy index along z
        INTEGER(iwp) ::  k_ppp !< dummy index along z
-       
+
        LOGICAL      ::  flag_set !< steering variable for advection flags
-   
+
+       ALLOCATE( advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        advc_flags_m = 0
-
 !
 !--    Set advc_flags_m to steer the degradation of the advection scheme in advec_ws
@@ -850,8 +869,8 @@
        INTEGER(iwp) ::  k_pp  !< dummy index along z
        INTEGER(iwp) ::  k_ppp !< dummy index along z
-       
+
        INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::&
                                                   advc_flag !< flag array to control order of scalar advection
-       
+
        LOGICAL ::  flag_set     !< steering variable for advection flags
        LOGICAL ::  non_cyclic_l !< flag that indicates non-cyclic boundary on the left
@@ -859,5 +878,5 @@
        LOGICAL ::  non_cyclic_r !< flag that indicates non-cyclic boundary on the right
        LOGICAL ::  non_cyclic_s !< flag that indicates non-cyclic boundary on the south
-       
+
        LOGICAL, OPTIONAL ::  extensive_degrad !< flag indicating that extensive degradation is required, e.g. for 
                                               !< passive scalars nearby topography along the horizontal directions,
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 4413)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 4414)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! - Remove deprecated topography arrays nzb_s_inner, nzb_u_inner, etc.
+! - Move initialization of boundary conditions and multigrid into an extra
+!   module interface.
+! 
+! 4386 2020-01-27 15:07:30Z Giersch
 ! Allocation statements, comments, naming of variables revised and _wp added to 
 ! real type values
@@ -123,49 +128,47 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE init_grid
- 
-    USE advec_ws,                                                              &
-        ONLY:  ws_init_flags_momentum, ws_init_flags_scalar
 
     USE arrays_3d,                                                             &
         ONLY:  dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, x, xu, y, yv, zu, zw
-        
+
     USE control_parameters,                                                    &
-        ONLY:  bc_lr_cyc, bc_ns_cyc,                                           &
-               bc_dirichlet_l,                                                 &
-               bc_dirichlet_n,                                                 &
-               bc_dirichlet_r,                                                 &
-               bc_dirichlet_s,                                                 &
-               bc_radiation_l,                                                 &
-               bc_radiation_n,                                                 &
-               bc_radiation_r,                                                 &
-               bc_radiation_s,                                                 &
-               constant_flux_layer, dz, dz_max, dz_stretch_factor,             &
+        ONLY:  constant_flux_layer, dz, dz_max, dz_stretch_factor,             &
                dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,&
                dz_stretch_level_end_index, dz_stretch_level_start_index,       &
                dz_stretch_level_start, ibc_uv_b, message_string,               &
-               momentum_advec, number_stretch_level_end,                       &
-               number_stretch_level_start, ocean_mode, psolver, scalar_advec,  &
-               symmetry_flag, topography, use_surface_fluxes
-         
+               number_stretch_level_end,                                       &
+               number_stretch_level_start,                                     &
+               ocean_mode,                                                     &
+               psolver,                                                        &
+               symmetry_flag,                                                  &
+               topography,                                                     &
+               use_surface_fluxes
+
     USE grid_variables,                                                        &
         ONLY:  ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner
-        
+
     USE indices,                                                               &
-        ONLY:  advc_flags_m,                                                   &
-               advc_flags_s,                                                   &
-               nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
-               nzb, nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer,              &
-               nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,                 &
-               nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner,             &
-               nzb_w_outer, nzt, topo_top_ind, topo_min_level
-    
+        ONLY:  nbgp,                                                           &
+               nx,                                                             &
+               nxl,                                                            &
+               nxlg,                                                           &
+               nxr,                                                            &
+               nxrg,                                                           &
+               ny,                                                             &
+               nyn,                                                            &
+               nyng,                                                           &
+               nys,                                                            &
+               nysg,                                                           &
+               nz,                                                             &
+               nzb,                                                            &
+               nzb_diff,                                                       &
+               nzb_max,                                                        &
+               nzt,                                                            &
+               topo_top_ind,                                                   &
+               topo_min_level
+
     USE kinds
 
     USE pegrid
-
-    USE poismg_noopt_mod
-
-    USE surface_mod,                                                           &
-        ONLY:  init_bc
 
     USE vertical_nesting_mod,                                                  &
@@ -182,7 +185,4 @@
     INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
     INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
-                                      
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local  !< index for topography top at cell-center
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp    !< dummy to calculate topography indices on u- and v-grid
 
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
@@ -644,29 +644,5 @@
 !
 !-- Set flags to mask topography on the grid. 
-    CALL set_topo_flags( topo )    
-!
-!-- Calculate wall flag arrays for the multigrid method.
-!-- Please note, wall flags are only applied in the non-optimized version.
-    IF ( psolver == 'multigrid_noopt' )  CALL poismg_noopt_init
-
-!
-!-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e. 
-!-- to decrease the numerical stencil appropriately. The order of the scheme
-!-- is degraded near solid walls as well as near non-cyclic inflow and outflow
-!-- boundaries. Do this separately for momentum and scalars.
-    IF ( momentum_advec == 'ws-scheme' )  THEN
-       ALLOCATE( advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       CALL ws_init_flags_momentum
-    ENDIF
-    IF ( scalar_advec == 'ws-scheme'   )  THEN
-       ALLOCATE( advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       advc_flags_s = 0
-       
-       CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l,        &
-                                  bc_dirichlet_n  .OR.  bc_radiation_n,        &
-                                  bc_dirichlet_r  .OR.  bc_radiation_r,        &
-                                  bc_dirichlet_s  .OR.  bc_radiation_s,        &
-                                  advc_flags_s )
-    ENDIF
+    CALL set_topo_flags( topo )
 
 !
@@ -704,70 +680,5 @@
     topo_min_level = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
 #endif
-!
-!-- Initialize boundary conditions via surface type
-    CALL init_bc
-
-!
-!-- Allocate and set topography height arrays required for data output
-    IF ( TRIM( topography ) /= 'flat' )  THEN
-!
-!--    Allocate and set the arrays containing the topography height
-       IF ( nxr == nx  .AND.  nyn /= ny )  THEN
-          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn) )                             
-          ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn) )
-       ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
-          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1) )                             
-          ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn+1) )
-       ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
-          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1) ) 
-          ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn+1) )
-       ELSE
-          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn) )
-          ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn) )
-       ENDIF
-
-       zu_s_inner   = 0.0_wp
-       zw_w_inner   = 0.0_wp
-!
-!--    Determine local topography height on scalar and w-grid. Note, setting
-!--    lateral boundary values is not necessary, realized via wall_flags_static_0 
-!--    array. Further, please note that loop bounds are different from 
-!--    nxl to nxr and nys to nyn on south and right model boundary, hence, 
-!--    use intrinsic lbound and ubound functions to infer array bounds.
-       DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
-          DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
-!
-!--          Topography height on scalar grid. Therefore, determine index of 
-!--          upward-facing surface element on scalar grid. 
-             zu_s_inner(i,j) = zu(topo_top_ind(j,i,0))
-!
-!--          Topography height on w grid. Therefore, determine index of 
-!--          upward-facing surface element on w grid.
-             zw_w_inner(i,j) = zw(topo_top_ind(j,i,3))
-          ENDDO
-       ENDDO
-    ENDIF
-
-!
-!-- In the following, calculate 2D index arrays. Note, these will be removed
-!-- soon. 
-!-- Allocate outer and inner index arrays for topography and set
-!-- defaults.                    
-    ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_s_outer(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_u_inner(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_u_outer(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_v_inner(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_v_outer(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_w_inner(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_w_outer(nysg:nyng,nxlg:nxrg) )                               
-    ALLOCATE( nzb_diff_s_inner(nysg:nyng,nxlg:nxrg) )                          
-    ALLOCATE( nzb_diff_s_outer(nysg:nyng,nxlg:nxrg) )                          
-    ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )                                 
-    ALLOCATE( nzb_tmp(nysg:nyng,nxlg:nxrg) )
-!
-!-- Initialize 2D-index arrays. Note, these will be removed soon! 
-    nzb_local(nys:nyn,nxl:nxr) = topo_top_ind(nys:nyn,nxl:nxr,0)
-    CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
+
 !
 !-- Check topography for consistency with model domain. Therefore, use
@@ -792,10 +703,4 @@
        ENDIF
     ENDIF
-
-    nzb_s_inner = nzb;  nzb_s_outer = nzb
-    nzb_u_inner = nzb;  nzb_u_outer = nzb
-    nzb_v_inner = nzb;  nzb_v_outer = nzb
-    nzb_w_inner = nzb;  nzb_w_outer = nzb
-
 !
 !-- Define vertical gridpoint from (or to) which on the usual finite difference
@@ -807,170 +712,42 @@
     ENDIF
 
-    nzb_diff_s_inner = nzb_diff;  nzb_diff_s_outer = nzb_diff
-!
-!-- Set Neumann conditions for topography. Will be removed soon.
-    IF ( .NOT. bc_ns_cyc )  THEN
-       IF ( nys == 0  )  THEN
-          DO  i = 1, nbgp  
-             nzb_local(nys-i,:)   = nzb_local(nys,:)
-          ENDDO
-       ELSEIF ( nyn == ny )  THEN
-          DO  i = 1, nbgp  
-             nzb_local(ny+i,:) = nzb_local(ny,:)
-          ENDDO
-       ENDIF
-    ENDIF
-
-    IF ( .NOT. bc_lr_cyc )  THEN
-       IF ( nxl == 0  )  THEN
-          DO  i = 1, nbgp  
-             nzb_local(:,nxl-i)   = nzb_local(:,nxl)
-          ENDDO
-       ELSEIF ( nxr == nx )  THEN
-          DO  i = 1, nbgp  
-             nzb_local(:,nx+i) = nzb_local(:,nx)
-          ENDDO 
-       ENDIF          
-    ENDIF
-!
-!-- Initialization of 2D index arrays, will be removed soon!
-!-- Initialize nzb_s_inner and nzb_w_inner
-    nzb_s_inner = nzb_local
-    nzb_w_inner = nzb_local
-
-!
-!-- Initialize remaining index arrays:
-!-- first pre-initialize them with nzb_s_inner...
-    nzb_u_inner = nzb_s_inner
-    nzb_u_outer = nzb_s_inner
-    nzb_v_inner = nzb_s_inner
-    nzb_v_outer = nzb_s_inner
-    nzb_w_outer = nzb_s_inner
-    nzb_s_outer = nzb_s_inner
-
-!
-!-- nzb_s_outer: 
-!-- extend nzb_local east-/westwards first, then north-/southwards
-    nzb_tmp = nzb_local
-    DO  j = nys, nyn
-       DO  i = nxl, nxr
-          nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i),             &
-                              nzb_local(j,i+1) )
+    IF ( TRIM( topography ) /= 'flat' )  THEN
+!
+!--    Allocate and set the arrays containing the topography height (for output
+!--    reasons only).
+       IF ( nxr == nx  .AND.  nyn /= ny )  THEN
+          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
+                    zw_w_inner(nxl:nxr+1,nys:nyn) )
+       ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
+          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
+                    zw_w_inner(nxl:nxr,nys:nyn+1) )
+       ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
+          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
+                    zw_w_inner(nxl:nxr+1,nys:nyn+1) )
+       ELSE
+          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
+                    zw_w_inner(nxl:nxr,nys:nyn) )
+       ENDIF
+
+       zu_s_inner   = 0.0_wp
+       zw_w_inner   = 0.0_wp
+!
+!--    Determine local topography height on scalar and w-grid. Note, setting
+!--    lateral boundary values is not necessary, realized via wall_flags_static_0
+!--    array. Further, please note that loop bounds are different from
+!--    nxl to nxr and nys to nyn on south and right model boundary, hence,
+!--    use intrinsic lbound and ubound functions to infer array bounds.
+       DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
+          DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
+!
+!--          Topography height on scalar grid. Therefore, determine index of
+!--          upward-facing surface element on scalar grid.
+             zu_s_inner(i,j) = zu(topo_top_ind(j,i,0))
+!
+!--          Topography height on w grid. Therefore, determine index of
+!--          upward-facing surface element on w grid.
+             zw_w_inner(i,j) = zw(topo_top_ind(j,i,3))
+          ENDDO
        ENDDO
-    ENDDO
-       
-    CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
-     
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-          nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
-                                  nzb_tmp(j+1,i) )
-       ENDDO
-!
-!--    non-cyclic boundary conditions (overwritten by call of
-!--    exchange_horiz_2d_int below in case of cyclic boundary conditions)
-       IF ( nys == 0 )  THEN
-          j = -1
-          nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
-       ENDIF
-       IF ( nyn == ny )  THEN
-          j = ny + 1
-          nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
-       ENDIF
-    ENDDO
-!
-!-- nzb_w_outer: 
-!-- identical to nzb_s_outer
-    nzb_w_outer = nzb_s_outer
-!
-!-- nzb_u_inner: 
-!-- extend nzb_local rightwards only
-    nzb_tmp = nzb_local
-    DO  j = nys, nyn
-       DO  i = nxl, nxr
-          nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
-       ENDDO
-    ENDDO
-       
-    CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
-       
-    nzb_u_inner = nzb_tmp
-!
-!-- nzb_u_outer: 
-!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-          nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i),             &
-                                  nzb_tmp(j+1,i) )
-       ENDDO
-!
-!--    non-cyclic boundary conditions (overwritten by call of
-!--    exchange_horiz_2d_int below in case of cyclic boundary conditions)
-       IF ( nys == 0 )  THEN
-          j = -1
-          nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
-       ENDIF
-       IF ( nyn == ny )  THEN
-          j = ny + 1
-          nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
-       ENDIF
-    ENDDO
-
-!
-!-- nzb_v_inner:
-!-- extend nzb_local northwards only
-    nzb_tmp = nzb_local
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-          nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
-       ENDDO
-    ENDDO
-       
-    CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )      
-    nzb_v_inner = nzb_tmp
-
-!
-!-- nzb_v_outer: 
-!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
-    DO  j = nys, nyn
-       DO  i = nxl, nxr
-          nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i),                &
-                                  nzb_tmp(j,i+1) )
-       ENDDO
-!
-!--    non-cyclic boundary conditions (overwritten by call of
-!--    exchange_horiz_2d_int below in case of cyclic boundary conditions)
-       IF ( nxl == 0 )  THEN
-          i = -1
-          nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
-       ENDIF
-       IF ( nxr == nx )  THEN
-          i = nx + 1
-          nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
-       ENDIF
-    ENDDO
-
-!
-!-- Exchange of lateral boundary values (parallel computers) and cyclic
-!-- boundary conditions, if applicable.
-!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
-!-- they do not require exchange and are not included here.
-    CALL exchange_horiz_2d_int( nzb_u_inner, nys, nyn, nxl, nxr, nbgp )
-    CALL exchange_horiz_2d_int( nzb_u_outer, nys, nyn, nxl, nxr, nbgp )
-    CALL exchange_horiz_2d_int( nzb_v_inner, nys, nyn, nxl, nxr, nbgp )
-    CALL exchange_horiz_2d_int( nzb_v_outer, nys, nyn, nxl, nxr, nbgp )
-    CALL exchange_horiz_2d_int( nzb_w_outer, nys, nyn, nxl, nxr, nbgp )
-    CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp )
-
-!
-!-- Set the individual index arrays which define the k index from which on
-!-- the usual finite difference form (which does not use surface fluxes) is
-!-- applied 
-    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
-       nzb_diff_s_inner   = nzb_s_inner + 2
-       nzb_diff_s_outer   = nzb_s_outer + 2
-    ELSE
-       nzb_diff_s_inner   = nzb_s_inner + 1
-       nzb_diff_s_outer   = nzb_s_outer + 1
     ENDIF
 !
@@ -2798,8 +2575,9 @@
                                     BTEST( wall_flags_total_0(:,:,:), ibit )   &
                                        ), DIM = 1                              &
-                                ) - 1                            
-        
+                                ) - 1
+
  END SUBROUTINE set_topo_flags
 
 
 
+
Index: /palm/trunk/SOURCE/module_interface.f90
===================================================================
--- /palm/trunk/SOURCE/module_interface.f90	(revision 4413)
+++ /palm/trunk/SOURCE/module_interface.f90	(revision 4414)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add module interface for basic initialization of numerics.
+! 
+! 4411 2020-02-18 14:28:02Z maronga
 ! Added output routines for WTM
 ! 
@@ -425,4 +428,7 @@
         ONLY:  particle_advection
 
+    USE poismg_noopt_mod,                                                      &
+        ONLY:  poismg_noopt_init
+
     USE plant_canopy_model_mod,                                                &
          ONLY: pcm_parin,                                                      &
@@ -487,4 +493,7 @@
                surface_data_output_wrd_global
 
+    USE surface_mod,                                                           &
+        ONLY:  init_bc
+
     USE synthetic_turbulence_generator_mod,                                    &
         ONLY:  stg_parin,                                                      &
@@ -562,4 +571,5 @@
        module_interface_init,                                                  &
        module_interface_init_checks,                                           &
+       module_interface_init_numerics,                                         &
        module_interface_init_output,                                           &
        module_interface_header,                                                &
@@ -620,4 +630,8 @@
        MODULE PROCEDURE module_interface_init_checks
     END INTERFACE module_interface_init_checks
+
+    INTERFACE module_interface_init_numerics
+       MODULE PROCEDURE module_interface_init_numerics
+    END INTERFACE module_interface_init_numerics
 
     INTERFACE module_interface_init_output
@@ -1110,4 +1124,22 @@
 ! Description:
 ! ------------
+!> Initialize boundary conditions and numerical schemes.
+!------------------------------------------------------------------------------!
+ SUBROUTINE module_interface_init_numerics
+
+!
+!-- Initialize boundary conditions via surface type
+    CALL init_bc
+!
+!-- Calculate wall flag arrays for the multigrid method.
+!-- Please note, wall flags are only applied in the non-optimized version.
+    CALL poismg_noopt_init
+ 
+ END SUBROUTINE module_interface_init_numerics
+
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Initialize data output
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 4413)
+++ /palm/trunk/SOURCE/modules.f90	(revision 4414)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! - nzb_diff_s_inner, nzb_diff_s_outer, nzb_inner,nzb_outer, nzb_s_inner, 
+!   nzb_s_outer, nzb_u_inner, nzb_u_outer, nzb_v_inner, nzb_v_outer, 
+!   nzb_w_inner, nzb_w_outer
+!
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
@@ -1033,5 +1039,5 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zu_s_inner  !< height of topography top on scalar grid 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_w_inner  !< height of topography top on w grid 
-                                              
+
     SAVE
 
@@ -1093,16 +1099,4 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner   !< number of horizontal grid points which are non-topography
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mg_loc_ind        !< internal array to store index bounds of all PEs of that multigrid level where data is collected to PE0
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_diff_s_inner  !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_diff_s_outer  !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_inner         !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_outer         !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_s_inner       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_s_outer       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_u_inner       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_u_outer       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_v_inner       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_v_outer       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_w_inner       !< will be removed
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_w_outer       !< will be removed
 
     INTEGER(iwp), DIMENSION(:,:,:), POINTER ::  flags  !< pointer to wall_flags_1-10
Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 4413)
+++ /palm/trunk/SOURCE/palm.f90	(revision 4414)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Call to module_interface_init_numerics
+! 
+! 4400 2020-02-10 20:32:41Z suehring
 ! Add interface to initialize data output with dom
 ! 
@@ -93,5 +96,6 @@
 
     USE module_interface,                                                      &
-        ONLY:  module_interface_init_output,                                   &
+        ONLY:  module_interface_init_numerics,                                 &
+               module_interface_init_output,                                   &
                module_interface_last_actions
 
@@ -112,5 +116,5 @@
         ONLY:  nested_run, pmci_child_initialize, pmci_init,                   &
                pmci_modelconfiguration, pmci_parent_initialize
-               
+
     USE surface_data_output_mod,                                               &
         ONLY:  surface_data_output_last_action
@@ -261,4 +265,8 @@
 !-- topography information if required
     CALL init_grid
+!
+!-- Initialize boundary conditions and numerics such as the multigrid solver or 
+!-- the advection routine
+    CALL module_interface_init_numerics
 !
 !-- Read global attributes if available.  
Index: /palm/trunk/SOURCE/poismg_noopt_mod.f90
===================================================================
--- /palm/trunk/SOURCE/poismg_noopt_mod.f90	(revision 4413)
+++ /palm/trunk/SOURCE/poismg_noopt_mod.f90	(revision 4414)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Remove double-declared use only construct.
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
@@ -84,14 +87,7 @@
 
     SUBROUTINE poismg_noopt( r )
-    
 
        USE arrays_3d,                                                          &
            ONLY:  d, p_loc
-
-       USE control_parameters,                                                 &
-           ONLY:  bc_lr_cyc, bc_ns_cyc, gathered_size, grid_level_count,       &
-                  ibc_p_t, maximum_grid_level, message_string, mgcycles,       &
-                  mg_cycles, mg_switch_to_pe0_level, residual_limit,           &
-                  subdomain_size
 
        USE control_parameters,                                                 &
@@ -1655,5 +1651,6 @@
 
        USE control_parameters,                                                 &
-           ONLY:  bc_lr_cyc, bc_ns_cyc, masking_method, maximum_grid_level
+           ONLY:  bc_lr_cyc, bc_ns_cyc, masking_method, maximum_grid_level,    &
+                  psolver
 
        USE indices,                                                            &
@@ -1678,4 +1675,6 @@
 
        INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo_tmp
+
+       IF ( psolver /= 'multigrid_noopt' )  RETURN
 !
 !--    Gridpoint increment of the current level. 
