Index: /palm/trunk/SCRIPTS/subjob
===================================================================
--- /palm/trunk/SCRIPTS/subjob	(revision 1574)
+++ /palm/trunk/SCRIPTS/subjob	(revision 1575)
@@ -23,5 +23,5 @@
 # Current revisions:
 # ------------------
-# 
+# mpp2-queues added to lccrayh
 # 
 # Former revisions:
@@ -417,5 +417,5 @@
                  esac;;
         (lccrayh) case  $ndq  in
-                     (dataq|mpp1q|mpp1testq|smp1q|smp1testq|specialm1q)   error=false;;
+                     (dataq|mpp1q|mpp1testq|mpp2q|mpp2testq|smp1q|smp1testq|specialm1q)   error=false;;
                      (*)                                     error=true;;
                  esac;;
Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 1574)
+++ /palm/trunk/SOURCE/Makefile	(revision 1575)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# +poismg_fast
 # 
 # Former revisions:
@@ -220,5 +220,5 @@
 	init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
 	interaction_droplets_ptq.f90 land_surface_model.f90 local_flush.f90 \
-        local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
+   local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
 	local_tremain_ini.f90 lpm.f90 lpm_advec.f90 lpm_boundary_conds.f90 \
 	lpm_calc_liquid_water_content.f90 lpm_collision_kernels.f90 \
@@ -231,7 +231,7 @@
 	lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \
 	ls_forcing.f90 message.f90 microphysics.f90 modules.f90 mod_kinds.f90 \
-	mod_particle_attributes.f90 netcdf.f90 nudging.f90 \
-	package_parin.f90 palm.f90 parin.f90 plant_canopy_model.f90 poisfft.f90 \
-	poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 production_e.f90 \
+	mod_particle_attributes.f90 netcdf.f90 nudging.f90 package_parin.f90 \
+	palm.f90 parin.f90 plant_canopy_model.f90 poisfft.f90 poismg.f90 \
+	poismg_fast.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 production_e.f90 \
 	prognostic_equations.f90 progress_bar.f90 radiation_model.f90 \
 	random_function.f90 random_gauss.f90 random_generator_parallel.f90 \
@@ -249,6 +249,5 @@
 	user_data_output_mask.f90 user_define_netcdf_grid.f90 \
 	user_dvrp_coltab.f90 user_header.f90 user_init.f90 \
-	user_init_3d_model.f90 user_init_grid.f90 \
-        user_init_land_surface.f90 \
+	user_init_3d_model.f90 user_init_grid.f90 user_init_land_surface.f90 \
 	user_init_plant_canopy.f90 user_last_actions.f90 user_lpm_advec.f90 \
 	user_lpm_init.f90 user_lpm_set_attributes.f90 user_module.f90 \
@@ -413,6 +412,7 @@
 poisfft.o: modules.o cpulog.o fft_xy.o mod_kinds.o tridia_solver.o
 poismg.o: modules.o cpulog.o mod_kinds.o
+poismg_fast.o: modules.o cpulog.o mod_kinds.o
 prandtl_fluxes.o: modules.o mod_kinds.o
-pres.o: modules.o cpulog.o mod_kinds.o poisfft.o
+pres.o: modules.o cpulog.o mod_kinds.o poisfft.o poismg_fast.o
 print_1d.o: modules.o cpulog.o mod_kinds.o
 production_e.o: modules.o mod_kinds.o wall_fluxes.o
Index: /palm/trunk/SOURCE/calc_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/calc_spectra.f90	(revision 1574)
+++ /palm/trunk/SOURCE/calc_spectra.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! adjustments for psolver-queries
 ! 
 ! Former revisions:
@@ -113,5 +113,5 @@
 !
 !-- Reallocate array d in required size
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
        DEALLOCATE( d )
        ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1574)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! multigrid_fast added as allowed pressure solver
 ! 
 ! Former revisions:
@@ -716,6 +716,6 @@
 !-- Check whether there are any illegal values
 !-- Pressure solver:
-    IF ( psolver /= 'poisfft'  .AND. &
-         psolver /= 'sor'  .AND.  psolver /= 'multigrid' )  THEN
+    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
+         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_fast' )  THEN
        message_string = 'unknown solver for perturbation pressure: psolver' // &
                         ' = "' // TRIM( psolver ) // '"'
@@ -723,5 +723,5 @@
     ENDIF
 
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
        IF ( cycle_mg == 'w' )  THEN
           gamma_mg = 2
@@ -1715,5 +1715,5 @@
 !-- and tools do not work with non-cyclic boundary conditions.
     IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
-       IF ( psolver /= 'multigrid' )  THEN
+       IF ( psolver(1:9) /= 'multigrid' )  THEN
           message_string = 'non-cyclic lateral boundaries do not allow ' // &
                            'psolver = "' // TRIM( psolver ) // '"'
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 1574)
+++ /palm/trunk/SOURCE/header.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! adjustments for psolver-queries, output of seed_follows_topography
 ! 
 ! Former revisions:
@@ -227,4 +227,5 @@
                particles_per_point, pdx, pdy, pdz,  psb, psl, psn, psr, pss,   &
                pst, radius, radius_classes, random_start_position,             &
+               seed_follows_topography,                                        &
                total_number_of_particles, use_particle_tails,                  &
                use_sgs_for_particles, total_number_of_tails,                   &
@@ -424,6 +425,6 @@
     ELSEIF ( psolver == 'sor' )  THEN
        WRITE ( io, 112 )  nsor_ini, nsor, omega_sor
-    ELSEIF ( psolver == 'multigrid' )  THEN
-       WRITE ( io, 135 )  cycle_mg, maximum_grid_level, ngsrb
+    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
+       WRITE ( io, 135 )  TRIM(psolver), cycle_mg, maximum_grid_level, ngsrb
        IF ( mg_cycles == -1 )  THEN
           WRITE ( io, 140 )  residual_limit
@@ -1788,4 +1789,5 @@
        IF ( use_sgs_for_particles )  WRITE ( io, 488 )  dt_min_part
        IF ( random_start_position )  WRITE ( io, 481 )
+       IF ( seed_follows_topography )  WRITE ( io, 496 )
        IF ( particles_per_point > 1 )  WRITE ( io, 489 )  particles_per_point
        WRITE ( io, 495 )  total_number_of_particles
@@ -1907,5 +1909,5 @@
 133 FORMAT ('     Precipitation parameterization via Kessler-Scheme')
 134 FORMAT (' --> Additional prognostic equation for a passive scalar')
-135 FORMAT (' --> Solve perturbation pressure via multigrid method (', &
+135 FORMAT (' --> Solve perturbation pressure via ',A,' method (', &
                   A,'-cycle)'/ &
             '     number of grid levels:                   ',I2/ &
@@ -2326,4 +2328,6 @@
                     F8.2,' s'/)
 495 FORMAT ('       Number of particles in total domain: ',I10/)
+496 FORMAT ('       Initial vertical particle positions are interpreted ', &
+                    'as relative to the given topography')
 500 FORMAT (//' 1D-Model parameters:'/                           &
               ' -------------------'//                           &
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1574)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! adjustments for psolver-queries
 ! 
 ! Former revisions:
@@ -356,5 +356,5 @@
     IF ( psolver == 'sor' )  THEN
        ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-    ELSEIF ( psolver == 'multigrid' )  THEN
+    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
 !
 !--    For performance reasons, multigrid is using one ghost layer only
@@ -1135,5 +1135,5 @@
 !
 !--    In case of iterative solvers, p must get an initial value
-       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
+       IF ( psolver(1:9) == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
 
 !
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 1574)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! adjustments for psolver-queries
 ! 
 ! Former revisions:
@@ -332,5 +332,5 @@
 !-- containing with appropriate grid information is created for these
 !-- solvers.
-    IF ( psolver /= 'multigrid' )  THEN
+    IF ( psolver(1:9) /= 'multigrid' )  THEN
        ALLOCATE( ddzu_pres(1:nzt+1) )
        ddzu_pres = ddzu
@@ -341,5 +341,5 @@
 !-- In case of multigrid method, compute grid lengths and grid factors for the
 !-- grid levels
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
 
        ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
@@ -1047,5 +1047,5 @@
 !
 !-- Calculate wall flag arrays for the multigrid method
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
 !
 !--    Gridpoint increment of the current level
Index: /palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- /palm/trunk/SOURCE/init_pegrid.f90	(revision 1574)
+++ /palm/trunk/SOURCE/init_pegrid.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! adjustments for psolver-queries, calculation of ngp_xz added
 ! 
 ! Former revisions:
@@ -761,5 +761,5 @@
 !-- Calculate number of grid levels necessary for the multigrid poisson solver
 !-- as well as the gridpoint indices on each level
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
 
 !
@@ -977,6 +977,6 @@
 !-- Do these calculations for the model grid and (if necessary) also
 !-- for the coarser grid levels used in the multigrid method
-    ALLOCATE ( ngp_yz(0:maximum_grid_level), type_xz(0:maximum_grid_level),&
-               type_yz(0:maximum_grid_level) )
+    ALLOCATE ( ngp_xz(0:maximum_grid_level), ngp_yz(0:maximum_grid_level),     &
+               type_xz(0:maximum_grid_level), type_yz(0:maximum_grid_level) )
 
     nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
@@ -1007,9 +1007,10 @@
 !
 !-- Definition of MPI-datatypes for multigrid method (coarser level grids)
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
 !    
 !--    Definition of MPI-datatyoe as above, but only 1 ghost level is used
        DO  i = maximum_grid_level, 1 , -1
 
+          ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
           ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
 
@@ -1129,5 +1130,5 @@
 !
 !-- Allocate wall flag arrays used in the multigrid solver
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
 
        DO  i = maximum_grid_level, 1, -1
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 1574)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! initial vertical particle position is allowed to follow the topography
 ! 
 ! Former revisions:
@@ -101,5 +101,6 @@
 
     USE indices,                                                               &
-        ONLY:  nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb, nzt
+        ONLY:  nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb,    &
+               nzb_w_inner, nzt
 
     USE kinds
@@ -126,5 +127,5 @@
                 prt_count, psb, psl, psn, psr, pss, pst,                       &
                 radius, random_start_position, read_particles_from_restartfile,&
-                skip_particles_for_tail, sort_count,                           &
+                seed_follows_topography, skip_particles_for_tail, sort_count,  &
                 tail_mask, total_number_of_particles, total_number_of_tails,   &
                 use_particle_tails, use_sgs_for_particles,                     &
@@ -649,5 +650,5 @@
                    pos_x = psl(i)
 
-                   DO WHILE ( pos_x <= psr(i) )
+            xloop: DO WHILE ( pos_x <= psr(i) )
 
                       IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.  &
@@ -738,7 +739,20 @@
                                tmp_particle%tail_id    = 0
                             ENDIF
+!
+!--                         Determine the grid indices of the particle position
                             ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
                             jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
                             kp = tmp_particle%z / dz + 1 + offset_ocean_nzt_m1
+
+                            IF ( seed_follows_topography )  THEN
+!
+!--                            Particle height is given relative to topography
+                               kp = kp + nzb_w_inner(jp,ip)
+                               tmp_particle%z = tmp_particle%z + zw(kp)
+                               IF ( kp > nzt )  THEN
+                                  pos_x = pos_x + pdx(i)
+                                  CYCLE xloop
+                               ENDIF
+                            ENDIF
 
                             local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
@@ -758,5 +772,5 @@
                       pos_x = pos_x + pdx(i)
 
-                   ENDDO
+                   ENDDO xloop
 
                 ENDIF
Index: /palm/trunk/SOURCE/mod_particle_attributes.f90
===================================================================
--- /palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 1574)
+++ /palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +seed_follows_topography
 ! 
 ! Former revisions:
@@ -69,4 +69,5 @@
                 particle_advection = .FALSE., random_start_position = .FALSE., &
                 read_particles_from_restartfile = .TRUE.,                      &
+                seed_follows_topography = .FALSE.,                             &
                 uniform_particles = .TRUE., use_kernel_tables = .FALSE.,       &
                 use_particle_tails = .FALSE., use_sgs_for_particles = .FALSE., &
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1574)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +ngp_xz
 ! 
 ! Former revisions:
@@ -1304,5 +1304,5 @@
 
 
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp_yz, type_xz, type_yz
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp_xz, ngp_yz, type_xz, type_yz
 
     LOGICAL ::  left_border_pe  = .FALSE., north_border_pe = .FALSE., &
Index: /palm/trunk/SOURCE/package_parin.f90
===================================================================
--- /palm/trunk/SOURCE/package_parin.f90	(revision 1574)
+++ /palm/trunk/SOURCE/package_parin.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +seed_follows_topography in particles_par
 ! 
 ! Former revisions:
@@ -129,7 +129,7 @@
                psb, psl, psn, psr, pss, pst, radius, radius_classes,           &
                random_start_position, read_particles_from_restartfile,         &
-               skip_particles_for_tail, use_particle_tails,                    &
-               use_sgs_for_particles, vertical_particle_advection,             &
-               write_particle_statistics
+               seed_follows_topography, skip_particles_for_tail,               &
+               use_particle_tails, use_sgs_for_particles,                      &
+               vertical_particle_advection, write_particle_statistics
 
     USE plant_canopy_model_mod,                                                &
@@ -212,4 +212,5 @@
                                   radius_classes, random_start_position,       &
                                   read_particles_from_restartfile,             &
+                                  seed_follows_topography,                     &
                                   skip_particles_for_tail,                     &
                                   use_particle_tails, use_sgs_for_particles,   &
Index: /palm/trunk/SOURCE/poismg_fast.f90
===================================================================
--- /palm/trunk/SOURCE/poismg_fast.f90	(revision 1575)
+++ /palm/trunk/SOURCE/poismg_fast.f90	(revision 1575)
@@ -0,0 +1,2880 @@
+ MODULE poismg_mod
+
+!--------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the terms
+! of the GNU General Public License as published by the Free Software Foundation,
+! either version 3 of the License, or (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2014 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
+!            still does not give the expected speedup! Further work required.
+!
+! Current revisions:
+! -----------------
+! Initial revision.
+! Routine re-written and optimised based on poismg.
+!
+! Following optimisations have been made:
+! - vectorisation (for Intel-CPUs) of the red-black algorithm by resorting
+!   array elements with even and odd indices
+! - explicit boundary conditions for building walls removed (solver is
+!   running through the buildings
+! - reduced data transfer in case of ghost point exchange, because only
+!   "red" or "black" data points need to be exchanged. This is not applied
+!   for coarser grid levels, since for then the transfer time is latency bound
+! 
+! Former revisions:
+! -----------------
+! $Id$
+!
+!
+!
+! Description:
+! ------------
+! Solves the Poisson equation for the perturbation pressure with a multigrid
+! V- or W-Cycle scheme.
+!
+! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
+! September 2000 - July 2001. It has been optimised for speed by Klaus
+! Ketelsen in November 2014.
+!------------------------------------------------------------------------------!
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point_s
+    USE kinds
+    USE MPI
+    USE pegrid
+
+    PRIVATE
+
+    INTEGER, SAVE                             ::  ind_even_odd    !: border index between even and odd k index
+    INTEGER, DIMENSION(:), SAVE, ALLOCATABLE  ::  even_odd_level  !: stores ind_even_odd for all MG levels
+
+    REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE ::  f1_mg_b, f2_mg_b, f3_mg_b  !: blocked version of f1_mg ...
+
+    INTERFACE poismg_fast
+       MODULE PROCEDURE poismg_fast
+    END INTERFACE poismg_fast
+
+    INTERFACE sort_k_to_even_odd_blocks
+       MODULE PROCEDURE sort_k_to_even_odd_blocks
+       MODULE PROCEDURE sort_k_to_even_odd_blocks_int
+       MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
+    END INTERFACE sort_k_to_even_odd_blocks
+
+    PUBLIC poismg_fast
+
+ CONTAINS
+
+    SUBROUTINE poismg_fast( r )
+
+       USE arrays_3d,                                                          &
+           ONLY:  d, p_loc
+
+       USE control_parameters,                                                 &
+           ONLY:  gathered_size, grid_level, grid_level_count,                 &
+                  maximum_grid_level, message_string, mgcycles, mg_cycles,     &
+                  mg_switch_to_pe0_level, residual_limit, subdomain_size
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,&
+                  nyng, nyn_mg, nzb, nzt, nzt_mg
+
+       IMPLICIT NONE
+
+       REAL(wp) ::  maxerror          !:
+       REAL(wp) ::  maximum_mgcycles  !:
+       REAL(wp) ::  residual_norm     !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r  !:
+
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p3  !:
+
+
+       CALL cpu_log( log_point_s(29), 'poismg_fast', 'start' )
+!
+!--    Initialize arrays and variables used in this subroutine
+
+!--    If the number of grid points of the gathered grid, which is collected
+!--    on PE0, is larger than the number of grid points of an PE, than array
+!--    p3 will be enlarged.
+       IF ( gathered_size > subdomain_size )  THEN
+          ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(            &
+                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
+                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
+                    mg_switch_to_pe0_level)+1) )
+       ELSE
+          ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+
+       p3 = 0.0_wp
+
+ 
+!
+!--    Ghost boundaries have to be added to divergence array.
+!--    Exchange routine needs to know the grid level!
+       grid_level = maximum_grid_level
+       CALL exchange_horiz( d, 1)
+       d(nzb,:,:) = d(nzb+1,:,:)
+
+!
+!--    Initiation of the multigrid scheme. Does n cycles until the
+!--    residual is smaller than the given limit. The accuracy of the solution
+!--    of the poisson equation will increase with the number of cycles.
+!--    If the number of cycles is preset by the user, this number will be
+!--    carried out regardless of the accuracy.
+       grid_level_count =  0
+       mgcycles         =  0
+       IF ( mg_cycles == -1 )  THEN
+          maximum_mgcycles = 0
+          residual_norm    = 1.0_wp
+       ELSE
+          maximum_mgcycles = mg_cycles
+          residual_norm    = 0.0_wp
+       ENDIF
+
+!
+!--    Initial settings for sorting k-dimension from sequential order (alternate
+!--    even/odd) into blocks of even and odd or vice versa
+       CALL init_even_odd_blocks
+
+!
+!--    Sort input arrays in even/odd blocks along k-dimension
+       CALL sort_k_to_even_odd_blocks( d, grid_level )
+       CALL sort_k_to_even_odd_blocks( p_loc, grid_level )
+
+!
+!--    The complete multigrid cycles are running in block mode, i.e. over
+!--    seperate data blocks of even and odd indices
+       DO WHILE ( residual_norm > residual_limit  .OR. &
+                  mgcycles < maximum_mgcycles )
+ 
+          CALL next_mg_level_fast( d, p_loc, p3, r)
+
+!
+!--       Calculate the residual if the user has not preset the number of
+!--       cycles to be performed
+          IF ( maximum_mgcycles == 0 )  THEN
+             CALL resid_fast( d, p_loc, r )
+             maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
+
+#if defined( __parallel )
+             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+                CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL,      &
+                                    MPI_SUM, comm2d, ierr)
+#else
+                residual_norm = maxerror
+#endif
+             residual_norm = SQRT( residual_norm )
+          ENDIF
+
+          mgcycles = mgcycles + 1
+
+!
+!--       If the user has not limited the number of cycles, stop the run in case
+!--       of insufficient convergence
+          IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
+             message_string = 'no sufficient convergence within 1000 cycles'
+             CALL message( 'poismg_fast', 'PA0283', 1, 2, 0, 6, 0 )
+          ENDIF
+
+       ENDDO
+
+       DEALLOCATE( p3 )
+!
+!--    Result has to be sorted back from even/odd blocks to sequential order
+       CALL sort_k_to_sequential( p_loc )
+!
+!--    Unset the grid level. Variable is used to determine the MPI datatypes for
+!--    ghost point exchange
+       grid_level = 0
+
+       CALL cpu_log( log_point_s(29), 'poismg_fast', 'stop' )
+
+    END SUBROUTINE poismg_fast
+
+
+
+    SUBROUTINE resid_fast( f_mg, p_mg, r )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Computes the residual of the perturbation pressure.
+!------------------------------------------------------------------------------!
+
+       USE arrays_3d,                                                          &
+           ONLY:  f1_mg, f2_mg, f3_mg
+
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
+                  inflow_n, inflow_r, inflow_s, masking_method, outflow_l,     &
+                  outflow_n, outflow_r, outflow_s, topography
+
+       USE grid_variables,                                                     &
+           ONLY:  ddx2_mg, ddy2_mg
+
+       USE indices,                                                            &
+           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
+                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
+                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
+                  nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i
+       INTEGER(iwp) ::  j
+       INTEGER(iwp) ::  k
+       INTEGER(iwp) ::  l
+       INTEGER(iwp) :: km1      !:
+       INTEGER(iwp) :: kp1      !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !:
+
+!
+!--    Calculate the residual
+       l = grid_level
+
+       CALL cpu_log( log_point_s(53), 'resid', 'start' )
+       IF ( topography == 'flat'  .OR.  masking_method )  THEN
+          !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
+          !$OMP DO
+          DO  i = nxl_mg(l), nxr_mg(l)
+             DO  j = nys_mg(l), nyn_mg(l)
+                !DIR$ IVDEP
+                DO k = ind_even_odd+1, nzt_mg(l)
+                   km1 = k-ind_even_odd-1
+                   kp1 = k-ind_even_odd
+                   r(k,j,i) = f_mg(k,j,i)                                 &
+                      - ddx2_mg(l) *                                      &
+                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                 &
+                      - ddy2_mg(l) *                                      &
+                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                  &
+                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                      &
+                      - f3_mg_b(k,l) * p_mg(km1,j,i)                      &
+                      + f1_mg_b(k,l) * p_mg(k,j,i)
+                ENDDO
+                !DIR$ IVDEP
+                DO k = nzb+1, ind_even_odd
+                   km1 = k+ind_even_odd
+                   kp1 = k+ind_even_odd+1
+                   r(k,j,i) = f_mg(k,j,i)                                 &
+                      - ddx2_mg(l) *                                      &
+                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                 &
+                      - ddy2_mg(l) *                                      &
+                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                  &
+                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                      &
+                      - f3_mg_b(k,l) * p_mg(km1,j,i)                      &
+                      + f1_mg_b(k,l) * p_mg(k,j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+          !$OMP END PARALLEL
+       ELSE
+!
+!--       Choose flag array of this level
+          SELECT CASE ( l )
+             CASE ( 1 )
+                flags => wall_flags_1
+             CASE ( 2 )
+                flags => wall_flags_2
+             CASE ( 3 )
+                flags => wall_flags_3
+             CASE ( 4 )
+                flags => wall_flags_4
+             CASE ( 5 )
+                flags => wall_flags_5
+             CASE ( 6 )
+                flags => wall_flags_6
+             CASE ( 7 )
+                flags => wall_flags_7
+             CASE ( 8 )
+                flags => wall_flags_8
+             CASE ( 9 )
+                flags => wall_flags_9
+             CASE ( 10 )
+                flags => wall_flags_10
+          END SELECT
+
+          !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
+          !$OMP DO
+          DO  i = nxl_mg(l), nxr_mg(l)
+             DO  j = nys_mg(l), nyn_mg(l)
+
+                !DIR$ IVDEP
+                DO  k = ind_even_odd+1, nzt_mg(l)
+                   km1 = k-ind_even_odd-1
+                   kp1 = k-ind_even_odd
+
+                   r(k,j,i) = f_mg(k,j,i)                                         &
+                    - ddx2_mg(l) *                                                &
+                    ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 ))   &
+                    + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
+                    - ddy2_mg(l)                                                  &
+                    * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
+                    + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
+                    - f2_mg(k,l) * p_mg(kp1,j,i)                                  &
+                    - f3_mg(k,l) *                                                &
+                      MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 ))   &
+                    + f1_mg(k,l) * p_mg(k,j,i)
+                ENDDO
+
+                !DIR$ IVDEP
+                DO  k = nzb+1, ind_even_odd
+                   km1 = k+ind_even_odd
+                   kp1 = k+ind_even_odd+1
+
+                   r(k,j,i) = f_mg(k,j,i)                                         &
+                    - ddx2_mg(l) *                                                &
+                    ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 ))   &
+                    + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
+                    - ddy2_mg(l)                                                  &
+                    * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
+                    + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
+                    - f2_mg(k,l) * p_mg(kp1,j,i)                                  &
+                    - f3_mg(k,l) *                                                &
+                      MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 ))   &
+                    + f1_mg(k,l) * p_mg(k,j,i)
+                ENDDO
+!
+!--             The residual within topography should be zero
+                WHERE( BTEST(flags(nzb+1:nzt_mg(l),j,i), 6) )
+                   r(nzb+1:nzt_mg(l),j,i) = 0.0
+                END WHERE
+
+             ENDDO
+          ENDDO
+          !$OMP END PARALLEL
+
+       ENDIF
+
+!
+!--    Horizontal boundary conditions
+       CALL exchange_horiz( r, 1)
+
+       IF ( .NOT. bc_lr_cyc )  THEN
+          IF ( inflow_l .OR. outflow_l )  r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
+          IF ( inflow_r .OR. outflow_r )  r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
+       ENDIF
+
+       IF ( .NOT. bc_ns_cyc )  THEN
+          IF ( inflow_n .OR. outflow_n )  r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
+          IF ( inflow_s .OR. outflow_s )  r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
+       ENDIF
+
+!
+!--    Boundary conditions at bottom and top of the domain.outflow_l, outflow_n,
+!--    These points are not handled by the above loop. Points may be within
+!--    buildings, but that doesn't matter.
+       IF ( ibc_p_b == 1 )  THEN
+          r(nzb,:,: ) = r(nzb+1,:,:)
+       ELSE
+          r(nzb,:,: ) = 0.0_wp
+       ENDIF
+
+       IF ( ibc_p_t == 1 )  THEN
+          r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
+       ELSE
+          r(nzt_mg(l)+1,:,: ) = 0.0_wp
+       ENDIF
+
+       CALL cpu_log( log_point_s(53), 'resid', 'stop' )
+
+    END SUBROUTINE resid_fast
+
+
+
+    SUBROUTINE restrict_fast( f_mg, r )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Interpolates the residual on the next coarser grid with "full weighting"
+! scheme
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
+                  inflow_n, inflow_r, inflow_s, masking_method, outflow_l,     &
+                  outflow_n, outflow_r, outflow_s, topography
+
+       USE indices,                                                            &
+           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
+                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
+                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
+                  nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i    !:
+       INTEGER(iwp) ::  ic   !:
+       INTEGER(iwp) ::  j    !:
+       INTEGER(iwp) ::  jc   !:
+       INTEGER(iwp) ::  k    !:
+       INTEGER(iwp) ::  kc   !:
+       INTEGER(iwp) ::  l    !:
+       INTEGER(iwp) ::  km1  !:
+       INTEGER(iwp) ::  kp1  !:
+
+       REAL(wp) ::  rkjim    !:
+       REAL(wp) ::  rkjip    !:
+       REAL(wp) ::  rkjmi    !:
+       REAL(wp) ::  rkjmim   !:
+       REAL(wp) ::  rkjmip   !:
+       REAL(wp) ::  rkjpi    !:
+       REAL(wp) ::  rkjpim   !:
+       REAL(wp) ::  rkjpip   !:
+       REAL(wp) ::  rkmji    !:
+       REAL(wp) ::  rkmjim   !:
+       REAL(wp) ::  rkmjip   !:
+       REAL(wp) ::  rkmjmi   !:
+       REAL(wp) ::  rkmjmim  !:
+       REAL(wp) ::  rkmjmip  !:
+       REAL(wp) ::  rkmjpi   !:
+       REAL(wp) ::  rkmjpim  !:
+       REAL(wp) ::  rkmjpip  !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,                         &
+                           nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,      &
+                           nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r !:
+
+!
+!--    Interpolate the residual
+       l = grid_level
+
+       CALL cpu_log( log_point_s(54), 'restrict', 'start' )
+
+       IF ( topography == 'flat'  .OR.  masking_method )  THEN
+!
+!--       No wall treatment
+          !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
+          DO  ic = nxl_mg(l), nxr_mg(l)
+             i = 2*ic
+             !$OMP DO SCHEDULE( STATIC )
+             DO  jc = nys_mg(l), nyn_mg(l)
+!
+!--             Calculation for the first point along k
+                j  = 2*jc
+                k   = ind_even_odd+1
+                kp1 = k-ind_even_odd
+                kc  = k-ind_even_odd
+                f_mg(kc,jc,ic) =  1.0_wp / 64.0_wp * (                         &
+                                  8.0_wp * r(k,j,i)                            &
+                               +  4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
+                                             r(k,j+1,i)     + r(k,j-1,i)     ) &
+                               +  2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
+                                             r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
+                               + 16.0_wp * r(k,j,i)                            &
+                               +  4.0_wp * r(kp1,j,i)                          &
+                               +  2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
+                                             r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
+                               +           ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
+                                             r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
+                                                      )
+!
+!--             Calculation for the other points along k
+                !DIR$ IVDEP
+                DO k = ind_even_odd+2, nzt_mg(l+1)    ! Fine grid at this point
+                   km1 = k-ind_even_odd-1
+                   kp1 = k-ind_even_odd
+                   kc  = k-ind_even_odd               ! Coarse grid index
+
+                   f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
+                                  8.0_wp * r(k,j,i)                            &
+                                + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
+                                             r(k,j+1,i)     + r(k,j-1,i)     ) &
+                                + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
+                                             r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
+                                + 4.0_wp * r(km1,j,i)                          &
+                                + 2.0_wp * ( r(km1,j,i-1)   + r(km1,j,i+1)   + &
+                                             r(km1,j+1,i)   + r(km1,j-1,i)   ) &
+                                +          ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
+                                             r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
+                                + 4.0_wp * r(kp1,j,i)                          &
+                                + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
+                                             r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
+                                +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
+                                             r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
+                                           )
+                ENDDO
+             ENDDO
+             !$OMP ENDDO nowait
+          ENDDO
+          !$OMP END PARALLEL
+
+       ELSE
+!
+!--       Choose flag array of the upper level
+          SELECT CASE ( l+1 )
+             CASE ( 1 )
+                flags => wall_flags_1
+             CASE ( 2 )
+                flags => wall_flags_2
+             CASE ( 3 )
+                flags => wall_flags_3
+             CASE ( 4 )
+                flags => wall_flags_4
+             CASE ( 5 )
+                flags => wall_flags_5
+             CASE ( 6 )
+                flags => wall_flags_6
+             CASE ( 7 )
+                flags => wall_flags_7
+             CASE ( 8 )
+                flags => wall_flags_8
+             CASE ( 9 )
+                flags => wall_flags_9
+             CASE ( 10 )
+                flags => wall_flags_10
+          END SELECT
+    
+          !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1,rkjim,rkjip,rkjpi,    &
+          !$OMP                   rkjmim,rkjpim,rkjmip,rkjpip,rkmji,rkmjim,    &
+          !$OMP                   rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim,        &
+          !$OMP                   rkmjmip,rkmjpip)
+          !$OMP DO
+          DO  ic = nxl_mg(l), nxr_mg(l)
+             i = 2*ic
+             DO  jc = nys_mg(l), nyn_mg(l)
+                j = 2*jc
+                !DIR$ IVDEP
+                DO k = ind_even_odd+1, nzt_mg(l+1)      !fine grid at this point
+
+                    km1 = k-ind_even_odd-1
+                    kp1 = k-ind_even_odd
+                    kc  = k-ind_even_odd                !Coarse grid index
+!
+!--                 Use implicit Neumann BCs if the respective gridpoint is
+!--                 inside the building
+                    rkjim   = MERGE(r(k,j,i),r(k,j,i-1),BTEST( flags(k,j,i-1), 6))
+                    rkjip   = MERGE(r(k,j,i),r(k,j,i+1),BTEST( flags(k,j,i+1), 6))
+                    rkjpi   = MERGE(r(k,j,i),r(k,j+1,i),BTEST( flags(k,j+1,i), 6))
+                    rkjmi   = MERGE(r(k,j,i),r(k,j-1,i),BTEST( flags(k,j-1,i), 6))
+
+                    rkjmim  = MERGE(r(k,j,i),r(k,j-1,i-1),BTEST( flags(k,j-1,i-1), 6))
+                    rkjpim  = MERGE(r(k,j,i),r(k,j+1,i-1),BTEST( flags(k,j+1,i-1), 6))
+                    rkjmip  = MERGE(r(k,j,i),r(k,j-1,i+1),BTEST( flags(k,j-1,i+1), 6))
+                    rkjpip  = MERGE(r(k,j,i),r(k,j+1,i+1),BTEST( flags(k,j+1,i+1), 6))
+
+                    rkmji   = MERGE(r(k,j,i),r(km1,j,i)  ,BTEST( flags(km1,j,i)  , 6))
+                    rkmjim  = MERGE(r(k,j,i),r(km1,j,i-1),BTEST( flags(km1,j,i-1), 6))
+                    rkmjip  = MERGE(r(k,j,i),r(km1,j,i+1),BTEST( flags(km1,j,i+1), 6))
+                    rkmjpi  = MERGE(r(k,j,i),r(km1,j+1,i),BTEST( flags(km1,j+1,i), 6))
+                    rkmjmi  = MERGE(r(k,j,i),r(km1,j-1,i),BTEST( flags(km1,j-1,i), 6))
+
+                    rkmjmim = MERGE(r(k,j,i),r(km1,j-1,i-1),BTEST( flags(km1,j-1,i-1), 6))
+                    rkmjpim = MERGE(r(k,j,i),r(km1,j+1,i-1),BTEST( flags(km1,j+1,i-1), 6))
+                    rkmjmip = MERGE(r(k,j,i),r(km1,j-1,i+1),BTEST( flags(km1,j-1,i+1), 6))
+                    rkmjpip = MERGE(r(k,j,i),r(km1,j+1,i+1),BTEST( flags(km1,j+1,i+1), 6))
+
+                    f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
+                           8.0_wp * r(k,j,i)                                   &
+                         + 4.0_wp * ( rkjim   + rkjip   +  rkjpi   + rkjmi   ) &
+                         + 2.0_wp * ( rkjmim  + rkjpim  +  rkjmip  + rkjpip  ) &
+                         + 4.0_wp * rkmji                                      &
+                         + 2.0_wp * ( rkmjim  + rkmjip  +  rkmjpi  + rkmjmi  ) &
+                         +          ( rkmjmim + rkmjpim +  rkmjmip + rkmjpip ) &
+                         + 4.0_wp * r(kp1,j,i)                                 &
+                         + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   +        &
+                                      r(kp1,j+1,i)   + r(kp1,j-1,i)   )        &
+                         +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) +        &
+                                      r(kp1,j-1,i+1) + r(kp1,j+1,i+1) )        &
+                                                        )
+                ENDDO
+             ENDDO
+          ENDDO
+          !$OMP END PARALLEL
+
+       ENDIF
+
+!
+!--    Horizontal boundary conditions
+       CALL exchange_horiz( f_mg, 1)
+
+       IF ( .NOT. bc_lr_cyc )  THEN
+          IF (inflow_l .OR. outflow_l)  f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
+          IF (inflow_r .OR. outflow_r)  f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
+       ENDIF
+
+       IF ( .NOT. bc_ns_cyc )  THEN
+          IF (inflow_n .OR. outflow_n)  f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
+          IF (inflow_s .OR. outflow_s)  f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
+       ENDIF
+
+!
+!--    Boundary conditions at bottom and top of the domain.
+!--    These points are not handled by the above loop. Points may be within
+!--    buildings, but that doesn't matter.
+       IF ( ibc_p_b == 1 )  THEN
+          f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
+       ELSE
+          f_mg(nzb,:,: ) = 0.0_wp
+       ENDIF
+
+       IF ( ibc_p_t == 1 )  THEN
+          f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
+       ELSE
+          f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
+       ENDIF
+
+       CALL cpu_log( log_point_s(54), 'restrict', 'stop' )
+!
+!--    Why do we need a sorting here?
+       CALL sort_k_to_even_odd_blocks( f_mg , l)
+
+    END SUBROUTINE restrict_fast
+
+
+
+    SUBROUTINE prolong_fast( p, temp )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Interpolates the correction of the perturbation pressure 
+! to the next finer grid.
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
+                  inflow_n, inflow_r, inflow_s, outflow_l, outflow_n,          &
+                  outflow_r, outflow_s
+
+       USE indices,                                                            &
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i   !:
+       INTEGER(iwp) ::  j   !:
+       INTEGER(iwp) ::  k   !:
+       INTEGER(iwp) ::  l   !:
+       INTEGER(iwp) ::  kp1 !:
+       INTEGER(iwp) ::  ke  !: Index for prolog even
+       INTEGER(iwp) ::  ko  !: Index for prolog odd
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
+                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
+                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p  !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp  !:
+
+
+       CALL cpu_log( log_point_s(55), 'prolong', 'start' )
+
+!
+!--    First, store elements of the coarser grid on the next finer grid
+       l = grid_level
+       ind_even_odd = even_odd_level(grid_level-1)
+
+       !$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko)
+       !$OMP DO
+       DO  i = nxl_mg(l-1), nxr_mg(l-1)
+          DO  j = nys_mg(l-1), nyn_mg(l-1)
+
+             !DIR$ IVDEP
+             DO k = ind_even_odd+1, nzt_mg(l-1)
+                kp1 = k - ind_even_odd
+                ke  = 2 * ( k-ind_even_odd - 1 ) + 1
+                ko  = 2 * k - 1
+!
+!--             Points of the coarse grid are directly stored on the next finer
+!--             grid
+                temp(ko,2*j,2*i)   = p(k,j,i)
+!
+!--             Points between two coarse-grid points
+                temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
+                temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
+                temp(ke,2*j,2*i)   = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
+!
+!--             Points in the center of the planes stretched by four points
+!--             of the coarse grid cube
+                temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
+                                                   p(k,j+1,i) + p(k,j+1,i+1) )
+                temp(ke,2*j,2*i+1)   = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
+                                                   p(kp1,j,i) + p(kp1,j,i+1) )
+                temp(ke,2*j+1,2*i)   = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +   &
+                                                   p(kp1,j,i) + p(kp1,j+1,i) )
+!
+!--             Points in the middle of coarse grid cube
+                temp(ke,2*j+1,2*i+1) = 0.125_wp *                              &
+                                               ( p(k,j,i)     + p(k,j,i+1)   + &
+                                                 p(k,j+1,i)   + p(k,j+1,i+1) + &
+                                                 p(kp1,j,i)   + p(kp1,j,i+1) + &
+                                                 p(kp1,j+1,i) + p(kp1,j+1,i+1) )
+             ENDDO
+
+             !DIR$ IVDEP
+             DO k = nzb+1, ind_even_odd
+                kp1 = k + ind_even_odd + 1
+                ke  = 2 * k
+                ko  = 2 * ( k + ind_even_odd )
+!
+!--             Points of the coarse grid are directly stored on the next finer
+!--             grid
+                temp(ko,2*j,2*i)   = p(k,j,i)
+!
+!--             Points between two coarse-grid points
+                temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
+                temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
+                temp(ke,2*j,2*i)   = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
+!
+!--             Points in the center of the planes stretched by four points
+!--             of the coarse grid cube
+                temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
+                                                   p(k,j+1,i) + p(k,j+1,i+1) )
+                temp(ke,2*j,2*i+1)   = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
+                                                   p(kp1,j,i) + p(kp1,j,i+1) )
+                temp(ke,2*j+1,2*i)   = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +   &
+                                                   p(kp1,j,i) + p(kp1,j+1,i) )
+!
+!--             Points in the middle of coarse grid cube
+                temp(ke,2*j+1,2*i+1) = 0.125_wp *                              &
+                                               ( p(k,j,i)     + p(k,j,i+1)   + &
+                                                 p(k,j+1,i)   + p(k,j+1,i+1) + &
+                                                 p(kp1,j,i)   + p(kp1,j,i+1) + &
+                                                 p(kp1,j+1,i) + p(kp1,j+1,i+1) )
+             ENDDO
+
+          ENDDO
+       ENDDO
+       !$OMP END PARALLEL
+
+       ind_even_odd = even_odd_level(grid_level)
+!
+!--    Horizontal boundary conditions
+       CALL exchange_horiz( temp, 1)
+
+       IF ( .NOT. bc_lr_cyc )  THEN
+          IF (inflow_l .OR. outflow_l)  temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
+          IF (inflow_r .OR. outflow_r)  temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
+       ENDIF
+
+       IF ( .NOT. bc_ns_cyc )  THEN
+          IF (inflow_n .OR. outflow_n)  temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
+          IF (inflow_s .OR. outflow_s)  temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
+       ENDIF
+
+!
+!--    Bottom and top boundary conditions
+       IF ( ibc_p_b == 1 )  THEN
+          temp(nzb,:,: ) = temp(ind_even_odd+1,:,:)
+       ELSE
+          temp(nzb,:,: ) = 0.0_wp
+       ENDIF
+
+       IF ( ibc_p_t == 1 )  THEN
+          temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:)
+       ELSE
+          temp(nzt_mg(l)+1,:,: ) = 0.0_wp
+       ENDIF
+
+       CALL cpu_log( log_point_s(55), 'prolong', 'stop' )
+
+    END SUBROUTINE prolong_fast
+
+
+
+    SUBROUTINE redblack_fast( f_mg, p_mg )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
+! 3D-Red-Black decomposition (GS-RB) is used.
+!------------------------------------------------------------------------------!
+
+       USE arrays_3d,                                                          &
+           ONLY:  f1_mg, f2_mg, f3_mg
+
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
+                  inflow_n, inflow_r, inflow_s, masking_method, ngsrb,         &
+                  outflow_l, outflow_n, outflow_r, outflow_s, topography
+
+       USE grid_variables,                                                     &
+           ONLY:  ddx2_mg, ddy2_mg
+
+       USE indices,                                                            &
+           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
+                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
+                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
+                  nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) :: color    !:
+       INTEGER(iwp) :: i        !:
+       INTEGER(iwp) :: ic       !:
+       INTEGER(iwp) :: j        !:
+       INTEGER(iwp) :: jc       !:
+       INTEGER(iwp) :: jj       !:
+       INTEGER(iwp) :: k        !:
+       INTEGER(iwp) :: l        !:
+       INTEGER(iwp) :: n        !:
+       INTEGER(iwp) :: km1      !:
+       INTEGER(iwp) :: kp1      !:
+
+       LOGICAL :: unroll        !:
+
+       REAL(wp) ::  wall_left   !:
+       REAL(wp) ::  wall_north  !:
+       REAL(wp) ::  wall_right  !:
+       REAL(wp) ::  wall_south  !:
+       REAL(wp) ::  wall_total  !:
+       REAL(wp) ::  wall_top    !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
+
+       l = grid_level
+
+!
+!--    Choose flag array of this level
+       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
+
+          SELECT CASE ( l )
+             CASE ( 1 )
+                flags => wall_flags_1
+             CASE ( 2 )
+                flags => wall_flags_2
+             CASE ( 3 )
+                flags => wall_flags_3
+             CASE ( 4 )
+                flags => wall_flags_4
+             CASE ( 5 )
+                flags => wall_flags_5
+             CASE ( 6 )
+                flags => wall_flags_6
+             CASE ( 7 )
+                flags => wall_flags_7
+             CASE ( 8 )
+                flags => wall_flags_8
+             CASE ( 9 )
+                flags => wall_flags_9
+             CASE ( 10 )
+                flags => wall_flags_10
+          END SELECT
+
+       ENDIF
+
+       unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
+                  MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
+
+       DO  n = 1, ngsrb
+       
+          DO  color = 1, 2
+!
+!--          No wall treatment in case of flat topography
+             IF ( topography == 'flat'  .OR.  masking_method )  THEN
+
+                IF ( .NOT. unroll )  THEN
+
+                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
+!
+!--                Without unrolling of loops, no cache optimization
+                   !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
+                   !$OMP DO
+                   DO  i = nxl_mg(l), nxr_mg(l), 2
+                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+                         !DIR$ IVDEP
+                         DO  k = ind_even_odd+1, nzt_mg(l)
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+   
+                   !$OMP DO
+                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+                          !DIR$ IVDEP
+                          DO  k = ind_even_odd+1, nzt_mg(l)
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+  
+                   !$OMP DO
+                   DO  i = nxl_mg(l), nxr_mg(l), 2
+                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                   !$OMP DO
+                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   !$OMP END PARALLEL
+
+                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
+
+                ELSE
+!
+!--                Loop unrolling along y, only one i loop for better cache use
+                   CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
+
+                   !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
+                   !$OMP DO
+                   DO  ic = nxl_mg(l), nxr_mg(l), 2
+                      DO  jc = nys_mg(l), nyn_mg(l), 4
+                         i  = ic
+                         jj = jc+2-color
+                         !DIR$ IVDEP
+                         DO  k = ind_even_odd+1, nzt_mg(l)
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+                            j   = jj
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+
+                         i  = ic+1
+                         jj = jc+color-1
+                         !DIR$ IVDEP
+                         DO  k = ind_even_odd+1, nzt_mg(l)
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+                            j   = jj
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+
+                         i  = ic
+                         jj = jc+color-1
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+                            j   = jj
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+
+                         i  = ic+1
+                         jj = jc+2-color
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+                            j   = jj
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+
+                      ENDDO
+                   ENDDO
+                   !$OMP END PARALLEL
+
+                   CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
+
+                ENDIF
+
+             ELSE
+
+                IF ( .NOT. unroll )  THEN
+
+                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
+!
+!--                Without unrolling of loops, no cache optimization /
+!--                vectorization. Therefore, the non-vectorized IBITS function
+!--                is still used.
+                   !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
+                   !$OMP DO
+                   DO  i = nxl_mg(l), nxr_mg(l), 2
+                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+                         !DIR$ IVDEP
+                         DO  k = ind_even_odd+1, nzt_mg(l)
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+
+                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
+                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
+                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
+                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
+                               + f3_mg(k,l) *                                  &
+                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+   
+                   !$OMP DO
+                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+                          !DIR$ IVDEP
+                          DO  k = ind_even_odd+1, nzt_mg(l)
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+
+                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
+                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
+                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
+                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
+                               + f3_mg(k,l) *                                  &
+                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                   !$OMP DO
+                   DO  i = nxl_mg(l), nxr_mg(l), 2
+                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
+                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
+                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
+                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
+                               + f3_mg(k,l) *                                  &
+                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                   !$OMP DO
+                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+
+                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
+                               ddx2_mg(l) *                                    &
+                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
+                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
+                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
+                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
+                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
+                               + f3_mg(k,l) *                                  &
+                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
+                               - f_mg(k,j,i)               )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   !$OMP END PARALLEL
+
+                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
+
+                ELSE
+!
+!--                Loop unrolling along y, only one i loop for better cache use
+                   CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
+
+                   !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
+                   !$OMP DO
+                   DO  ic = nxl_mg(l), nxr_mg(l), 2
+                      DO  jc = nys_mg(l), nyn_mg(l), 4
+                         i  = ic
+                         jj = jc+2-color
+                         !DIR$ IVDEP
+                         DO  k = ind_even_odd+1, nzt_mg(l)
+
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+
+                            j = jj
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                         ENDDO
+
+                         i  = ic+1
+                         jj = jc+color-1
+                         !DIR$ IVDEP
+                         DO  k = ind_even_odd+1, nzt_mg(l)
+
+                            km1 = k-ind_even_odd-1
+                            kp1 = k-ind_even_odd
+
+                            j =jj
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                         ENDDO
+
+                         i  = ic
+                         jj = jc+color-1
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+
+                            j =jj
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                         ENDDO
+
+                         i  = ic+1
+                         jj = jc+2-color
+                         !DIR$ IVDEP
+                         DO  k = nzb+1, ind_even_odd
+
+                            km1 = k+ind_even_odd
+                            kp1 = k+ind_even_odd+1
+
+                            j =jj
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                            j = jj+2
+                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
+                                ddx2_mg(l) *                                                 &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
+                                + ddy2_mg(l) *                                               &
+                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
+                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
+                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
+                                + f3_mg(k,l) *                                               &
+                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
+                                - f_mg(k,j,i)               )
+
+                         ENDDO
+
+                      ENDDO
+                   ENDDO
+                   !$OMP END PARALLEL
+
+                   CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
+
+                ENDIF
+
+             ENDIF
+
+!
+!--          Horizontal boundary conditions
+             CALL special_exchange_horiz( p_mg, color )
+
+              IF ( .NOT. bc_lr_cyc )  THEN
+                IF ( inflow_l .OR. outflow_l )  THEN
+                   p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
+                ENDIF
+                IF ( inflow_r .OR. outflow_r )  THEN
+                   p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
+                ENDIF
+             ENDIF
+
+             IF ( .NOT. bc_ns_cyc )  THEN
+                IF ( inflow_n .OR. outflow_n )  THEN
+                   p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
+                ENDIF
+                IF ( inflow_s .OR. outflow_s )  THEN
+                   p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
+                ENDIF
+             ENDIF
+
+!
+!--          Bottom and top boundary conditions
+             IF ( ibc_p_b == 1 )  THEN
+                p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:)
+             ELSE
+                p_mg(nzb,:,: ) = 0.0_wp
+             ENDIF
+
+             IF ( ibc_p_t == 1 )  THEN
+                p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:)
+             ELSE
+                p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
+             ENDIF
+
+          ENDDO
+       ENDDO
+!
+!--    Set pressure within topography and at the topography surfaces
+       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
+
+          !$OMP PARALLEL PRIVATE (i,j,k,kp1,wall_left,wall_north,wall_right,   &
+          !$OMP                   wall_south,wall_top,wall_total)
+          !$OMP DO
+          DO  i = nxl_mg(l), nxr_mg(l)
+             DO  j = nys_mg(l), nyn_mg(l)
+                !DIR$ IVDEP
+                DO  k = ind_even_odd+1, nzt_mg(l)
+                   km1 = k-ind_even_odd-1
+                   kp1 = k-ind_even_odd
+!
+!--                First, set pressure inside topography to zero
+                   p_mg(k,j,i) = p_mg(k,j,i) *                                 &
+                                 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
+!
+!--                Second, determine if the gridpoint inside topography is
+!--                adjacent to a wall and set its value to a value given by the
+!--                average of those values obtained from Neumann boundary
+!--                condition
+                   wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
+                   wall_right = IBITS( flags(k,j,i+1), 4, 1 )
+                   wall_south = IBITS( flags(k,j-1,i), 3, 1 )
+                   wall_north = IBITS( flags(k,j+1,i), 2, 1 )
+                   wall_top   = IBITS( flags(kp1,j,i), 0, 1 )
+                   wall_total = wall_left + wall_right + wall_south +          &
+                                wall_north + wall_top
+
+                   IF ( wall_total > 0.0_wp )  THEN
+                      p_mg(k,j,i) = 1.0_wp / wall_total *                      &
+                                    ( wall_left  * p_mg(k,j,i-1) +             &
+                                      wall_right * p_mg(k,j,i+1) +             &
+                                      wall_south * p_mg(k,j-1,i) +             &
+                                      wall_north * p_mg(k,j+1,i) +             &
+                                      wall_top   * p_mg(kp1,j,i) )
+                   ENDIF
+                ENDDO
+
+                !DIR$ IVDEP
+                DO  k = nzb+1, ind_even_odd
+                   km1 = k+ind_even_odd
+                   kp1 = k+ind_even_odd+1
+!
+!--                First, set pressure inside topography to zero
+                   p_mg(k,j,i) = p_mg(k,j,i) *                                 &
+                                 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
+!
+!--                Second, determine if the gridpoint inside topography is
+!--                adjacent to a wall and set its value to a value given by the
+!--                average of those values obtained from Neumann boundary
+!--                condition
+                   wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
+                   wall_right = IBITS( flags(k,j,i+1), 4, 1 )
+                   wall_south = IBITS( flags(k,j-1,i), 3, 1 )
+                   wall_north = IBITS( flags(k,j+1,i), 2, 1 )
+                   wall_top   = IBITS( flags(kp1,j,i), 0, 1 )
+                   wall_total = wall_left + wall_right + wall_south +          &
+                                wall_north + wall_top
+
+                   IF ( wall_total > 0.0_wp )  THEN
+                      p_mg(k,j,i) = 1.0_wp / wall_total *                      &
+                                    ( wall_left  * p_mg(k,j,i-1) +             &
+                                      wall_right * p_mg(k,j,i+1) +             &
+                                      wall_south * p_mg(k,j-1,i) +             &
+                                      wall_north * p_mg(k,j+1,i) +             &
+                                      wall_top   * p_mg(kp1,j,i) )
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+          !$OMP END PARALLEL
+
+!
+!--       One more time horizontal boundary conditions
+          CALL exchange_horiz( p_mg, 1)
+
+       ENDIF
+
+    END SUBROUTINE redblack_fast
+
+
+
+    SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Sort k-Dimension from sequential into blocks of even and odd.
+! This is required to vectorize the red-black subroutine.
+! Version for 3D-REAL arrays
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  grid_level
+
+       USE indices,                                                            &
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN) ::  glevel
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1,                               &
+                           nys_mg(glevel)-1:nyn_mg(glevel)+1,                  &
+                           nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::  p_mg  !:
+!
+!--    Local variables
+       INTEGER(iwp) :: i        !:
+       INTEGER(iwp) :: j        !:
+       INTEGER(iwp) :: k        !:
+       INTEGER(iwp) :: l        !:
+       INTEGER(iwp) :: ind      !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !:
+
+
+       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
+
+       l = glevel
+       ind_even_odd = even_odd_level(l)
+
+       !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
+       !$OMP DO
+       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
+          DO  j = nys_mg(l)-1, nyn_mg(l)+1
+
+!
+!--          Sort the data with even k index
+             ind = nzb-1
+             DO  k = nzb, nzt_mg(l), 2
+                ind = ind + 1
+                tmp(ind) = p_mg(k,j,i)
+             ENDDO
+!
+!--          Sort the data with odd k index
+             DO  k = nzb+1, nzt_mg(l)+1, 2
+                ind = ind + 1
+                tmp(ind) = p_mg(k,j,i)
+             ENDDO
+
+             p_mg(:,j,i) = tmp
+
+          ENDDO
+       ENDDO
+       !$OMP END PARALLEL
+
+       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
+
+    END SUBROUTINE sort_k_to_even_odd_blocks
+
+
+
+    SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Sort k-Dimension from sequential into blocks of even and odd.
+! This is required to vectorize the red-black subroutine.
+! Version for 1D-REAL arrays
+!------------------------------------------------------------------------------!
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN) ::  glevel
+
+       REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) ::  f_mg
+       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  f_mg_b
+
+!
+!--    Local variables
+       INTEGER(iwp) :: ind   !:
+       INTEGER(iwp) :: k     !:
+
+
+       ind = nzb - 1
+!
+!--    Sort the data with even k index
+       DO  k = nzb, nzt_mg(glevel), 2
+          ind = ind + 1
+          IF ( k >= nzb+1  .AND.  k <= nzt_mg(glevel) )  THEN
+             f_mg_b(ind) = f_mg(k)
+          ENDIF
+       ENDDO
+!
+!--    Sort the data with odd k index
+       DO  k = nzb+1, nzt_mg(glevel)+1, 2
+          ind = ind + 1
+          IF( k >= nzb+1  .AND.  k <= nzt_mg(glevel) )  THEN
+             f_mg_b(ind) = f_mg(k)
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE sort_k_to_even_odd_blocks_1d
+
+
+
+    SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Sort k-Dimension from sequential into blocks of even and odd.
+! This is required to vectorize the red-black subroutine.
+! Version for 2D-INTEGER arrays
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  grid_level
+
+       USE indices,                                                            &
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN) ::  glevel
+
+       INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
+                               nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
+                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::  i_mg  !:
+!
+!--    Local variables
+       INTEGER(iwp) :: i        !:
+       INTEGER(iwp) :: j        !:
+       INTEGER(iwp) :: k        !:
+       INTEGER(iwp) :: l        !:
+       INTEGER(iwp) :: ind      !:
+       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp
+
+
+       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
+
+       l = glevel
+       ind_even_odd = even_odd_level(l)
+
+       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
+          DO  j = nys_mg(l)-1, nyn_mg(l)+1
+
+!
+!--          Sort the data with even k index
+             ind = nzb-1
+             DO  k = nzb, nzt_mg(l), 2
+                ind = ind + 1
+                tmp(ind) = i_mg(k,j,i)
+             ENDDO
+!
+!++          ATTENTION: Check reason for this error. Remove it or replace WRITE
+!++                     by PALM message
+             IF ( ind /= ind_even_odd )  THEN
+                WRITE (0,*) 'ERROR ==> illegal ind_even_odd ',ind,ind_even_odd,l
+                CALL MPI_ABORT(MPI_COMM_WORLD,i,j)
+             ENDIF
+!
+!--          Sort the data with odd k index
+             DO  k = nzb+1, nzt_mg(l)+1, 2
+                ind = ind + 1
+                tmp(ind) = i_mg(k,j,i)
+             ENDDO
+
+             i_mg(:,j,i) = tmp
+
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
+
+    END SUBROUTINE sort_k_to_even_odd_blocks_int
+
+
+
+    SUBROUTINE sort_k_to_sequential( p_mg )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Sort k-dimension from blocks of even and odd into sequential
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  grid_level
+
+       USE indices,                                                            &
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
+!
+!--    Local variables
+       INTEGER(iwp) :: i        !:
+       INTEGER(iwp) :: j        !:
+       INTEGER(iwp) :: k        !:
+       INTEGER(iwp) :: l        !:
+       INTEGER(iwp) :: ind      !:
+
+       REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) ::  tmp
+
+
+       l = grid_level
+
+       !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
+       !$OMP DO
+       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
+          DO  j = nys_mg(l)-1, nyn_mg(l)+1
+
+             ind = nzb - 1
+             tmp = p_mg(:,j,i)
+             DO  k = nzb, nzt_mg(l), 2
+                ind = ind + 1
+                p_mg(k,j,i) = tmp(ind)
+             ENDDO
+
+             DO  k = nzb+1, nzt_mg(l)+1, 2
+                ind = ind + 1
+                p_mg(k,j,i) = tmp(ind)
+             ENDDO
+          ENDDO
+       ENDDO
+       !$OMP END PARALLEL
+
+    END SUBROUTINE sort_k_to_sequential
+
+
+
+    SUBROUTINE mg_gather_fast( f2, f2_sub )
+
+       USE control_parameters,                                                 &
+           ONLY:  grid_level
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i       !:
+       INTEGER(iwp) ::  il      !:
+       INTEGER(iwp) ::  ir      !:
+       INTEGER(iwp) ::  j       !:
+       INTEGER(iwp) ::  jn      !:
+       INTEGER(iwp) ::  js      !:
+       INTEGER(iwp) ::  k       !:
+       INTEGER(iwp) ::  nwords  !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2    !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2_l  !:
+
+       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
+                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
+                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub  !:
+
+
+#if defined( __parallel )
+       CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
+
+       f2_l = 0.0_wp
+
+!
+!--    Store the local subdomain array on the total array
+       js = mg_loc_ind(3,myid)
+       IF ( south_border_pe )  js = js - 1
+       jn = mg_loc_ind(4,myid)
+       IF ( north_border_pe )  jn = jn + 1
+       il = mg_loc_ind(1,myid)
+       IF ( left_border_pe )   il = il - 1
+       ir = mg_loc_ind(2,myid)
+       IF ( right_border_pe )  ir = ir + 1
+       DO  i = il, ir
+          DO  j = js, jn
+             DO  k = nzb, nzt_mg(grid_level)+1
+                f2_l(k,j,i) = f2_sub(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Find out the number of array elements of the total array
+       nwords = SIZE( f2 )
+
+!
+!--    Gather subdomain data from all PEs
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
+                           f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),   &
+                           nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
+
+       CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
+#endif
+    
+    END SUBROUTINE mg_gather_fast
+
+
+
+    SUBROUTINE mg_scatter_fast( p2, p2_sub )
+!
+!--    TODO: It might be possible to improve the speed of this routine by using
+!--    non-blocking communication
+
+       USE control_parameters,                                                 &
+           ONLY:  grid_level
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  nwords  !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
+                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
+                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !:
+
+       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
+                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
+                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub  !:
+
+!
+!--    Find out the number of array elements of the subdomain array
+       nwords = SIZE( p2_sub )
+
+#if defined( __parallel )
+       CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
+
+       p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
+                     mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
+
+       CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
+#endif
+    
+    END SUBROUTINE mg_scatter_fast
+
+
+
+    RECURSIVE SUBROUTINE next_mg_level_fast( f_mg, p_mg, p3, r )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! This is where the multigrid technique takes place. V- and W- Cycle are
+! implemented and steered by the parameter "gamma". Parameter "nue" determines
+! the convergence of the multigrid iterative solution. There are nue times
+! RB-GS iterations. It should be set to "1" or "2", considering the time effort
+! one would like to invest. Last choice shows a very good converging factor,
+! but leads to an increase in computing time.
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir,      &
+                  gamma_mg, grid_level, grid_level_count, ibc_p_b, ibc_p_t,    &
+                  inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level,  &
+                  mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb, outflow_l,  &
+                  outflow_n, outflow_r, outflow_s
+
+       USE indices,                                                            &
+           ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,      &
+                  nyn_mg, nzb, nzt, nzt_mg
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i            !:
+       INTEGER(iwp) ::  j            !:
+       INTEGER(iwp) ::  k            !:
+       INTEGER(iwp) ::  nxl_mg_save  !:
+       INTEGER(iwp) ::  nxr_mg_save  !:
+       INTEGER(iwp) ::  nyn_mg_save  !:
+       INTEGER(iwp) ::  nys_mg_save  !:
+       INTEGER(iwp) ::  nzt_mg_save  !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg  !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg  !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3    !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r     !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
+                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
+                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2  !:
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
+                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
+                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !:
+
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub  !:
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p2_sub  !:
+
+!
+!--    Restriction to the coarsest grid
+    10 IF ( grid_level == 1 )  THEN
+
+!
+!--       Solution on the coarsest grid. Double the number of Gauss-Seidel
+!--       iterations in order to get a more accurate solution.
+          ngsrb = 2 * ngsrb
+
+          ind_even_odd = even_odd_level(grid_level)
+
+          CALL redblack_fast( f_mg, p_mg )
+
+          ngsrb = ngsrb / 2
+
+
+       ELSEIF ( grid_level /= 1 )  THEN
+
+          grid_level_count(grid_level) = grid_level_count(grid_level) + 1
+
+!
+!--       Solution on the actual grid level
+          ind_even_odd = even_odd_level(grid_level)
+
+          CALL redblack_fast( f_mg, p_mg )
+
+!
+!--       Determination of the actual residual
+          CALL resid_fast( f_mg, p_mg, r )
+
+!--       Restriction of the residual (finer grid values!) to the next coarser
+!--       grid. Therefore, the grid level has to be decremented now. nxl..nzt have
+!--       to be set to the coarse grid values, because these variables are needed
+!--       for the exchange of ghost points in routine exchange_horiz
+          grid_level = grid_level - 1
+
+          nxl = nxl_mg(grid_level)
+          nys = nys_mg(grid_level)
+          nxr = nxr_mg(grid_level)
+          nyn = nyn_mg(grid_level)
+          nzt = nzt_mg(grid_level)
+
+          IF ( grid_level == mg_switch_to_pe0_level )  THEN
+
+!
+!--          From this level on, calculations are done on PE0 only.
+!--          First, carry out restriction on the subdomain.
+!--          Therefore, indices of the level have to be changed to subdomain
+!--          values in between (otherwise, the restrict routine would expect
+!--          the gathered array)
+
+             nxl_mg_save = nxl_mg(grid_level)
+             nxr_mg_save = nxr_mg(grid_level)
+             nys_mg_save = nys_mg(grid_level)
+             nyn_mg_save = nyn_mg(grid_level)
+             nzt_mg_save = nzt_mg(grid_level)
+             nxl_mg(grid_level) = mg_loc_ind(1,myid)
+             nxr_mg(grid_level) = mg_loc_ind(2,myid)
+             nys_mg(grid_level) = mg_loc_ind(3,myid)
+             nyn_mg(grid_level) = mg_loc_ind(4,myid)
+             nzt_mg(grid_level) = mg_loc_ind(5,myid)
+             nxl = mg_loc_ind(1,myid)
+             nxr = mg_loc_ind(2,myid)
+             nys = mg_loc_ind(3,myid)
+             nyn = mg_loc_ind(4,myid)
+             nzt = mg_loc_ind(5,myid)
+
+             ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
+                              nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
+                              nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
+
+             CALL restrict_fast( f2_sub, r )
+
+!
+!--          Restore the correct indices of this level
+             nxl_mg(grid_level) = nxl_mg_save
+             nxr_mg(grid_level) = nxr_mg_save
+             nys_mg(grid_level) = nys_mg_save
+             nyn_mg(grid_level) = nyn_mg_save
+             nzt_mg(grid_level) = nzt_mg_save
+             nxl = nxl_mg(grid_level)
+             nxr = nxr_mg(grid_level)
+             nys = nys_mg(grid_level)
+             nyn = nyn_mg(grid_level)
+             nzt = nzt_mg(grid_level)
+!
+!--          Gather all arrays from the subdomains on PE0
+             CALL mg_gather_fast( f2, f2_sub )
+
+!
+!--          Set switch for routine exchange_horiz, that no ghostpoint exchange
+!--          has to be carried out from now on
+             mg_switch_to_pe0 = .TRUE.
+
+!
+!--          In case of non-cyclic lateral boundary conditions, both in- and
+!--          outflow conditions have to be used on all PEs after the switch,
+!--          because then they have the total domain.
+             IF ( bc_lr_dirrad )  THEN
+                inflow_l  = .TRUE.
+                inflow_r  = .FALSE.
+                outflow_l = .FALSE.
+                outflow_r = .TRUE.
+             ELSEIF ( bc_lr_raddir )  THEN
+                inflow_l  = .FALSE.
+                inflow_r  = .TRUE.
+                outflow_l = .TRUE.
+                outflow_r = .FALSE.
+             ENDIF
+
+             IF ( bc_ns_dirrad )  THEN
+                inflow_n  = .TRUE.
+                inflow_s  = .FALSE.
+                outflow_n = .FALSE.
+                outflow_s = .TRUE.
+             ELSEIF ( bc_ns_raddir )  THEN
+                inflow_n  = .FALSE.
+                inflow_s  = .TRUE.
+                outflow_n = .TRUE.
+                outflow_s = .FALSE.
+             ENDIF
+
+             DEALLOCATE( f2_sub )
+
+          ELSE
+
+             CALL restrict_fast( f2, r )
+
+             ind_even_odd = even_odd_level(grid_level)  ! must be after restrict
+
+          ENDIF
+
+          p2 = 0.0_wp
+
+!
+!--       Repeat the same procedure till the coarsest grid is reached
+          CALL next_mg_level_fast( f2, p2, p3, r )
+
+       ENDIF
+
+!
+!--    Now follows the prolongation
+       IF ( grid_level >= 2 )  THEN
+
+!
+!--       Prolongation of the new residual. The values are transferred
+!--       from the coarse to the next finer grid.
+          IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
+
+#if defined( __parallel )
+!
+!--          At this level, the new residual first has to be scattered from
+!--          PE0 to the other PEs
+             ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
+                       mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
+                       mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
+
+             CALL mg_scatter_fast( p2, p2_sub )
+
+!
+!--          Therefore, indices of the previous level have to be changed to
+!--          subdomain values in between (otherwise, the prolong routine would
+!--          expect the gathered array)
+             nxl_mg_save = nxl_mg(grid_level-1)
+             nxr_mg_save = nxr_mg(grid_level-1)
+             nys_mg_save = nys_mg(grid_level-1)
+             nyn_mg_save = nyn_mg(grid_level-1)
+             nzt_mg_save = nzt_mg(grid_level-1)
+             nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
+             nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
+             nys_mg(grid_level-1) = mg_loc_ind(3,myid)
+             nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
+             nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
+
+!
+!--          Set switch for routine exchange_horiz, that ghostpoint exchange
+!--          has to be carried again out from now on
+             mg_switch_to_pe0 = .FALSE.
+
+!
+!--          For non-cyclic lateral boundary conditions, restore the
+!--          in-/outflow conditions
+             inflow_l  = .FALSE.;  inflow_r  = .FALSE.
+             inflow_n  = .FALSE.;  inflow_s  = .FALSE.
+             outflow_l = .FALSE.;  outflow_r = .FALSE.
+             outflow_n = .FALSE.;  outflow_s = .FALSE.
+
+             IF ( pleft == MPI_PROC_NULL )  THEN
+                IF ( bc_lr_dirrad )  THEN
+                   inflow_l  = .TRUE.
+                ELSEIF ( bc_lr_raddir )  THEN
+                   outflow_l = .TRUE.
+                ENDIF
+             ENDIF
+
+             IF ( pright == MPI_PROC_NULL )  THEN
+                IF ( bc_lr_dirrad )  THEN
+                   outflow_r = .TRUE.
+                ELSEIF ( bc_lr_raddir )  THEN
+                   inflow_r  = .TRUE.
+                ENDIF
+             ENDIF
+
+             IF ( psouth == MPI_PROC_NULL )  THEN
+                IF ( bc_ns_dirrad )  THEN
+                   outflow_s = .TRUE.
+                ELSEIF ( bc_ns_raddir )  THEN
+                   inflow_s  = .TRUE.
+                ENDIF
+             ENDIF
+
+             IF ( pnorth == MPI_PROC_NULL )  THEN
+                IF ( bc_ns_dirrad )  THEN
+                   inflow_n  = .TRUE.
+                ELSEIF ( bc_ns_raddir )  THEN
+                   outflow_n = .TRUE.
+                ENDIF
+             ENDIF
+
+             CALL prolong_fast( p2_sub, p3 )
+
+!
+!--          Restore the correct indices of the previous level
+             nxl_mg(grid_level-1) = nxl_mg_save
+             nxr_mg(grid_level-1) = nxr_mg_save
+             nys_mg(grid_level-1) = nys_mg_save
+             nyn_mg(grid_level-1) = nyn_mg_save
+             nzt_mg(grid_level-1) = nzt_mg_save
+
+             DEALLOCATE( p2_sub )
+#endif
+
+          ELSE
+
+             CALL prolong_fast( p2, p3 )
+
+          ENDIF
+
+!
+!--       Computation of the new pressure correction. Therefore,
+!--       values from prior grids are added up automatically stage by stage.
+          DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
+             DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
+                DO  k = nzb, nzt_mg(grid_level)+1
+                   p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Relaxation of the new solution
+          CALL redblack_fast( f_mg, p_mg )
+
+       ENDIF
+
+
+!
+!--    The following few lines serve the steering of the multigrid scheme
+       IF ( grid_level == maximum_grid_level )  THEN
+
+          GOTO 20
+
+       ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
+                grid_level_count(grid_level) /= gamma_mg )  THEN
+
+          GOTO 10
+
+       ENDIF
+
+!
+!--    Reset counter for the next call of poismg_fast
+       grid_level_count(grid_level) = 0
+
+!
+!--    Continue with the next finer level. nxl..nzt have to be
+!--    set to the finer grid values, because these variables are needed for the
+!--    exchange of ghost points in routine exchange_horiz
+       grid_level = grid_level + 1
+       ind_even_odd = even_odd_level(grid_level)
+
+       nxl = nxl_mg(grid_level)
+       nxr = nxr_mg(grid_level)
+       nys = nys_mg(grid_level)
+       nyn = nyn_mg(grid_level)
+       nzt = nzt_mg(grid_level)
+
+    20 CONTINUE
+
+    END SUBROUTINE next_mg_level_fast
+
+
+
+    SUBROUTINE init_even_odd_blocks
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Initial settings for sorting k-dimension from sequential order (alternate
+! even/odd) into blocks of even and odd or vice versa
+!------------------------------------------------------------------------------!
+
+       USE arrays_3d,                                                          &
+           ONLY:  f1_mg, f2_mg, f3_mg
+
+       USE control_parameters,                                                 &
+           ONLY:  grid_level, masking_method, maximum_grid_level, topography
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzt, nzt_mg
+
+       USE indices,                                                            &
+           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
+                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
+                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
+                  nys_mg, nyn_mg, nzb, nzt_mg
+
+       IMPLICIT NONE
+!
+!--    Local variables
+       INTEGER(iwp) ::  i     !:
+       INTEGER(iwp) ::  l     !:
+
+       LOGICAL, SAVE ::  lfirst = .TRUE.
+
+
+       IF ( .NOT. lfirst )  RETURN
+
+       ALLOCATE( even_odd_level(maximum_grid_level) )
+
+       ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level),                        &
+                 f2_mg_b(nzb:nzt+1,maximum_grid_level),                        &
+                 f3_mg_b(nzb:nzt+1,maximum_grid_level) )
+
+!
+!--    Set border index between the even and odd block
+       DO  i = maximum_grid_level, 1, -1
+          even_odd_level(i) = nzt_mg(i) / 2
+       ENDDO
+
+!
+!--    Sort grid coefficients used in red/black scheme and for calculating the
+!--    residual to block (even/odd) structure
+       DO  l = maximum_grid_level, 1 , -1
+          CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l),   &
+                                          f1_mg_b(nzb:nzt_mg(grid_level)+1,l), &
+                                          l )
+          CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l),   &
+                                          f2_mg_b(nzb:nzt_mg(grid_level)+1,l), &
+                                          l )
+          CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l),   &
+                                          f3_mg_b(nzb:nzt_mg(grid_level)+1,l), &
+                                          l )
+       ENDDO
+
+!
+!--    Sort flags array for all levels to block (even/odd) structure
+       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
+
+          DO  l = maximum_grid_level, 1 , -1
+!
+!--          Assign the flag level to be calculated
+             SELECT CASE ( l )
+                CASE ( 1 )
+                   flags => wall_flags_1
+                CASE ( 2 )
+                   flags => wall_flags_2
+                CASE ( 3 )
+                   flags => wall_flags_3
+                CASE ( 4 )
+                   flags => wall_flags_4
+                CASE ( 5 )
+                   flags => wall_flags_5
+                CASE ( 6 )
+                   flags => wall_flags_6
+                CASE ( 7 )
+                   flags => wall_flags_7
+                CASE ( 8 )
+                   flags => wall_flags_8
+                CASE ( 9 )
+                   flags => wall_flags_9
+                CASE ( 10 )
+                   flags => wall_flags_10
+             END SELECT
+
+             CALL sort_k_to_even_odd_blocks( flags, l )
+
+          ENDDO
+
+       ENDIF
+
+       lfirst = .FALSE.
+
+     END SUBROUTINE init_even_odd_blocks
+
+
+
+     SUBROUTINE special_exchange_horiz ( p_mg, color )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Special exchange_horiz subroutine for use in redblack. Transfers only
+! "red" or "black" data points.
+!------------------------------------------------------------------------------!
+
+       USE control_parameters,                                                 &
+           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t,          &
+                  inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level,  &
+                  outflow_l, outflow_n, outflow_r, outflow_s,                  &
+                  synchronous_exchange
+
+       USE indices,                                                            &
+           ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,      &
+                  nyn_mg, nzb, nzt, nzt_mg
+
+       USE pegrid,                                                             &
+           ONLY:  ngp_xz, ngp_yz
+
+       IMPLICIT NONE
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
+                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
+
+       INTEGER(iwp), intent(IN) ::  color  !:
+!
+!--    Local variables
+       INTEGER(iwp) ::  i,i1,i2  !:
+       INTEGER(iwp) ::  j,j1,j2  !:
+       INTEGER(iwp) ::  k        !:
+       INTEGER(iwp) ::  l        !:
+       INTEGER(iwp) ::  jys      !:
+       INTEGER(iwp) ::  jyn      !:
+       INTEGER(iwp) ::  ixl      !:
+       INTEGER(iwp) ::  ixr      !:
+       logical      ::  sendrecv_in_background_save  !:
+       logical      ::  synchronous_exchange_save    !:
+
+       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
+                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
+                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  temp  !:
+
+
+       sendrecv_in_background_save = sendrecv_in_background
+       sendrecv_in_background      = .FALSE.
+       synchronous_exchange_save   = synchronous_exchange
+       synchronous_exchange        = .FALSE.
+
+       l = grid_level
+
+       ind_even_odd = even_odd_level(grid_level)
+
+       jys = nys_mg(grid_level-1)
+       jyn = nyn_mg(grid_level-1)
+       ixl = nxl_mg(grid_level-1)
+       ixr = nxr_mg(grid_level-1)
+
+!
+!--    Restricted transfer only on finer levels with enough data
+       IF ( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 )  THEN
+!
+!--       Handling the even k Values
+!--       Collecting data for the north - south exchange
+!--       Since only every second value has to be transfered, data are stored
+!--       on the next coarser grid level, because the arrays on that level
+!--       have just the required size
+          i1 = nxl_mg(grid_level-1)
+          i2 = nxl_mg(grid_level-1)
+
+          DO  i = nxl_mg(l), nxr_mg(l), 2
+             DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          grid_level = grid_level-1
+
+          nxl = nxl_mg(grid_level)
+          nys = nys_mg(grid_level)
+          nxr = nxr_mg(grid_level)
+          nyn = nyn_mg(grid_level)
+          nzt = nzt_mg(grid_level)
+
+          send_receive = 'ns'
+          CALL exchange_horiz( temp, 1 )
+
+          grid_level = grid_level+1
+
+          i1 = nxl_mg(grid_level-1)
+          i2 = nxl_mg(grid_level-1)
+
+          DO  i = nxl_mg(l), nxr_mg(l), 2
+             DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+!
+!--       Collecting data for the left - right exchange
+!--       Since only every second value has to be transfered, data are stored
+!--       on the next coarser grid level, because the arrays on that level
+!--       have just the required size
+          j1 = nys_mg(grid_level-1)
+          j2 = nys_mg(grid_level-1)
+
+          DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+             DO  i = nxl_mg(l), nxr_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
+             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          grid_level = grid_level-1
+          send_receive = 'lr'
+
+          CALL exchange_horiz( temp, 1 )
+
+          grid_level = grid_level+1
+
+          j1 = nys_mg(grid_level-1)
+          j2 = nys_mg(grid_level-1)
+
+          DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+             DO  i = nxl_mg(l), nxr_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,j,nxr_mg(l)+1)  = temp(k-ind_even_odd,j1,ixr+1)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,j,nxl_mg(l)-1)  = temp(k-ind_even_odd,j2,ixl-1)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
+             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,j,nxr_mg(l)+1)  = temp(k-ind_even_odd,j1,ixr+1)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = ind_even_odd+1, nzt_mg(l)
+                      p_mg(k,j,nxl_mg(l)-1)  = temp(k-ind_even_odd,j2,ixl-1)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+!
+!--       Now handling the even k values
+!--       Collecting data for the north - south exchange
+!--       Since only every second value has to be transfered, data are stored
+!--       on the next coarser grid level, because the arrays on that level
+!--       have just the required size
+          i1 = nxl_mg(grid_level-1)
+          i2 = nxl_mg(grid_level-1)
+
+          DO  i = nxl_mg(l), nxr_mg(l), 2
+             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,jys,i1) = p_mg(k,j,i)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,jyn,i2) = p_mg(k,j,i)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,jys,i1) = p_mg(k,j,i)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,jyn,i2) = p_mg(k,j,i)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          grid_level = grid_level-1
+
+          send_receive = 'ns'
+          CALL exchange_horiz( temp, 1 )
+
+          grid_level = grid_level+1
+
+          i1 = nxl_mg(grid_level-1)
+          i2 = nxl_mg(grid_level-1)
+
+          DO  i = nxl_mg(l), nxr_mg(l), 2
+             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+
+                IF ( j == nys_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
+                   ENDDO
+                   i1 = i1 + 1
+
+                ENDIF
+
+                IF ( j == nyn_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
+                   ENDDO
+                   i2 = i2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          j1 = nys_mg(grid_level-1)
+          j2 = nys_mg(grid_level-1)
+
+          DO  i = nxl_mg(l), nxr_mg(l), 2
+             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,j1,ixl) = p_mg(k,j,i)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,j2,ixr) = p_mg(k,j,i)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,j1,ixl) = p_mg(k,j,i)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      temp(k,j2,ixr) = p_mg(k,j,i)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          grid_level = grid_level-1
+
+          send_receive = 'lr'
+          CALL exchange_horiz( temp, 1 )
+
+          grid_level = grid_level+1
+
+          nxl = nxl_mg(grid_level)
+          nys = nys_mg(grid_level)
+          nxr = nxr_mg(grid_level)
+          nyn = nyn_mg(grid_level)
+          nzt = nzt_mg(grid_level)
+
+          j1 = nys_mg(grid_level-1)
+          j2 = nys_mg(grid_level-1)
+
+          DO  i = nxl_mg(l), nxr_mg(l), 2
+             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,j,nxr_mg(l)+1)  = temp(k,j1,ixr+1)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,j,nxl_mg(l)-1)  = temp(k,j2,ixl-1)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+
+                IF ( i == nxl_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,j,nxr_mg(l)+1)  = temp(k,j1,ixr+1)
+                   ENDDO
+                   j1 = j1 + 1
+
+                ENDIF
+
+                IF ( i == nxr_mg(l) )  THEN
+                   !DIR$ IVDEP
+                   DO  k = nzb+1, ind_even_odd
+                      p_mg(k,j,nxl_mg(l)-1)  = temp(k,j2,ixl-1)
+                   ENDDO
+                   j2 = j2 + 1
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+       ELSE
+!
+!--       Standard horizontal ghost boundary exchange for small coarse grid
+!--       levels, where the transfer time is latency bound
+          CALL exchange_horiz( p_mg, 1 )
+
+       ENDIF
+
+!
+!--    Reset values to default PALM setup
+       sendrecv_in_background = sendrecv_in_background_save
+       synchronous_exchange   = synchronous_exchange_save
+       send_receive = 'al'
+
+    END SUBROUTINE special_exchange_horiz
+
+ END MODULE poismg_mod
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 1574)
+++ /palm/trunk/SOURCE/pres.f90	(revision 1575)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! poismg_fast + respective module added, adjustments for psolver-queries
 ! 
 ! Former revisions:
@@ -114,4 +114,6 @@
         ONLY:  poisfft
 
+    USE poismg_mod
+
     USE statistics,                                                            &
         ONLY:  statistic_regions, sums_divnew_l, sums_divold_l, weight_pres,   &
@@ -144,5 +146,5 @@
 !-- Multigrid method expects array d to have one ghost layer.
 !-- 
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
      
        DEALLOCATE( d )
@@ -293,5 +295,5 @@
     CALL cpu_log( log_point_s(1), 'divergence', 'start' )
 
-    IF ( psolver == 'multigrid' )  THEN
+    IF ( psolver(1:9) == 'multigrid' )  THEN
        !$OMP PARALLEL DO SCHEDULE( STATIC )
        DO  i = nxl-1, nxr+1
@@ -388,13 +390,9 @@
 !
 !-- Compute the pressure perturbation solving the Poisson equation
-    IF ( psolver(1:7) == 'poisfft' )  THEN
+    IF ( psolver == 'poisfft' )  THEN
 
 !
 !--    Solve Poisson equation via FFT and solution of tridiagonal matrices
-       IF ( psolver == 'poisfft' )  THEN
-
-          CALL poisfft( d )
-
-       ENDIF
+       CALL poisfft( d )
 
 !
@@ -492,5 +490,5 @@
        tend = p_loc
 
-    ELSEIF ( psolver == 'multigrid' )  THEN
+    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
 
 !
@@ -511,5 +509,9 @@
        ENDIF
 
-       CALL poismg( tend )
+       IF ( psolver == 'multigrid' )  THEN
+          CALL poismg( tend )
+       ELSE
+          CALL poismg_fast( tend )
+       ENDIF
 
        IF ( gathered_size > subdomain_size )  THEN
