Index: /palm/trunk/SOURCE/cpulog_mod.f90
===================================================================
--- /palm/trunk/SOURCE/cpulog_mod.f90	(revision 4014)
+++ /palm/trunk/SOURCE/cpulog_mod.f90	(revision 4015)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! all reals changed to double precision in order to work with 32-bit working precision,
+! otherwise calculated time intervals would mostly give zero
+! 
+! 3885 2019-04-11 11:29:34Z kanani
 ! Changes related to global restructuring of location messages and introduction 
 ! of additional debug messages
@@ -151,9 +155,8 @@
  
 
-    USE control_parameters,                                                    &
-        ONLY: message_string, nr_timesteps_this_run, run_description_header,   &
-              synchronous_exchange
+    USE control_parameters,                                                                        &
+        ONLY: message_string, nr_timesteps_this_run, run_description_header, synchronous_exchange
                
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY: ngp_3d, nx, ny, nz
         
@@ -165,6 +168,6 @@
 
     PRIVATE
-    PUBLIC   cpu_log, cpu_log_barrierwait, cpu_log_nowait, cpu_statistics,     &
-             initial_wallclock_time, log_point, log_point_s
+    PUBLIC   cpu_log, cpu_log_barrierwait, cpu_log_nowait, cpu_statistics, initial_wallclock_time, &
+             log_point, log_point_s
 
     INTERFACE cpu_log
@@ -184,24 +187,22 @@
     LOGICAL, PARAMETER ::  cpu_log_nowait = .FALSE.       !< 
 
-    REAL(wp) ::  initial_wallclock_time  !< 
+    REAL(dp) ::  initial_wallclock_time  !<
 
     TYPE logpoint
-       REAL(wp)           ::  isum       !< 
-       REAL(wp)           ::  ivect      !< 
-       REAL(wp)           ::  mean       !< 
-       REAL(wp)           ::  mtime      !< 
-       REAL(wp)           ::  mtimevec   !< 
-       REAL(wp)           ::  sum        !< 
-       REAL(wp)           ::  vector     !< 
+       REAL(dp)           ::  isum       !<
+       REAL(dp)           ::  ivect      !<
+       REAL(dp)           ::  mean       !<
+       REAL(dp)           ::  mtime      !<
+       REAL(dp)           ::  mtimevec   !<
+       REAL(dp)           ::  sum        !<
+       REAL(dp)           ::  vector     !<
        INTEGER(iwp)       ::  counts     !< 
        CHARACTER (LEN=25) ::  place      !< 
     END TYPE logpoint
 
-    TYPE(logpoint), DIMENSION(100) ::  log_point = logpoint( 0.0_wp, 0.0_wp,   &
-                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
-                                       0, ' ' )
-    TYPE(logpoint), DIMENSION(100) ::  log_point_s = logpoint( 0.0_wp, 0.0_wp, &
-                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
-                                       0, ' ' )
+    TYPE(logpoint), DIMENSION(100) ::  log_point = logpoint( 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp,       &
+                                                             0.0_dp, 0.0_dp, 0.0_dp, 0, ' ' )
+    TYPE(logpoint), DIMENSION(100) ::  log_point_s = logpoint( 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp,     &
+                                                               0.0_dp, 0.0_dp, 0.0_dp, 0, ' ' )
 
     SAVE
@@ -225,6 +226,6 @@
        LOGICAL, SAVE     ::  first = .TRUE.     !< 
        
-       REAL(wp)          ::  mtime = 0.0_wp     !< 
-       REAL(wp)          ::  mtimevec = 0.0_wp  !< 
+       REAL(dp)          ::  mtime = 0.0_dp     !<
+       REAL(dp)          ::  mtimevec = 0.0_dp  !<
        TYPE(logpoint)    ::  log_event          !< 
 
@@ -238,6 +239,6 @@
           log_event%place = place
        ELSEIF ( log_event%place /= place )  THEN
-          WRITE( message_string, * ) 'wrong argument expected: ',              &
-                            TRIM(log_event%place), ' given: ',  TRIM( place )
+          WRITE( message_string, * ) 'wrong argument expected: ',                                  &
+                                     TRIM(log_event%place), ' given: ',  TRIM( place )
           CALL message( 'cpu_log', 'PA0174', 1, 2, 0, 6, 0 )
        ENDIF
@@ -256,5 +257,5 @@
 !--    PEs that have not yet finished
 #if defined( __parallel )
-       IF ( cpu_log_barrierwait  .AND.  wait_allowed  .AND.                    &
+       IF ( cpu_log_barrierwait  .AND.  wait_allowed  .AND.                                        &
             ( modus == 'start'  .OR.  modus == 'continue' ) )  THEN
           CALL MPI_BARRIER( comm2d, ierr )
@@ -265,5 +266,5 @@
 !--    Take current time
        CALL SYSTEM_CLOCK( count, count_rate )
-       mtime = REAL( count, KIND=wp ) / REAL( count_rate, KIND=wp )
+       mtime = REAL( count, KIND=dp ) / REAL( count_rate, KIND=dp )
 
 !
@@ -274,7 +275,7 @@
        ELSEIF ( modus == 'pause' )  THEN
           IF ( ( mtime - log_event%mtime ) < 0.0  .AND.  first )  THEN
-             WRITE( message_string, * ) 'negative time interval occured',      &
-                         '&PE',myid,' L=PAUSE "',TRIM(log_event%place),        &
-                         '" new=', mtime,' last=',log_event%mtime
+             WRITE( message_string, * ) 'negative time interval occured',                          &
+                                        '&PE',myid,' L=PAUSE "',TRIM(log_event%place),             &
+                                        '" new=', mtime,' last=',log_event%mtime
              CALL message( 'cpu_log', 'PA0176', 0, 1, -1, 6, 0 )
              first = .FALSE.
@@ -283,9 +284,8 @@
           log_event%ivect    = log_event%ivect  + mtimevec - log_event%mtimevec
        ELSEIF ( modus == 'stop' )  THEN
-          IF ( ( mtime - log_event%mtime + log_event%isum ) < 0.0  .AND.       &
-               first )  THEN
-             WRITE( message_string, * ) 'negative time interval occured',      &
-                       '&PE',myid,' L=STOP "',TRIM(log_event%place),'" new=',  &
-                       mtime,' last=',log_event%mtime,' isum=',log_event%isum
+          IF ( ( mtime - log_event%mtime + log_event%isum ) < 0.0  .AND. first )  THEN
+             WRITE( message_string, * ) 'negative time interval occured',                          &
+                                        '&PE',myid,' L=STOP "',TRIM(log_event%place),'" new=',     &
+                                        mtime,' last=',log_event%mtime,' isum=',log_event%isum
              CALL message( 'cpu_log', 'PA0177', 0, 1, -1, 6, 0 )
              first = .FALSE.
@@ -295,7 +295,7 @@
           log_event%sum      = log_event%sum  + log_event%mtime
           IF ( log_event%sum < 0.0  .AND.  first )  THEN
-             WRITE( message_string, * ) 'negative time interval occured',      &
-                       '&PE',myid,' L=STOP "',TRIM(log_event%place),'" sum=',  &
-                       log_event%sum,' mtime=',log_event%mtime
+             WRITE( message_string, * ) 'negative time interval occured',                          &
+                                        '&PE',myid,' L=STOP "',TRIM(log_event%place),'" sum=',     &
+                                        log_event%sum,' mtime=',log_event%mtime
              CALL message( 'cpu_log', 'PA0178', 0, 1, -1, 6, 0 )
              first = .FALSE.
@@ -303,9 +303,8 @@
           log_event%vector   = log_event%vector + log_event%mtimevec
           log_event%counts   = log_event%counts + 1
-          log_event%isum     = 0.0_wp
-          log_event%ivect    = 0.0_wp
+          log_event%isum     = 0.0_dp
+          log_event%ivect    = 0.0_dp
        ELSE
-          message_string = 'unknown modus of time measurement: ' //            &
-                           TRIM( modus )
+          message_string = 'unknown modus of time measurement: ' // TRIM( modus )
           CALL message( 'cpu_log', 'PA0179', 0, 1, -1, 6, 0 )
        ENDIF
@@ -331,12 +330,12 @@
        INTEGER(iwp)    ::  iii             !< 
        INTEGER(iwp)    ::  sender          !< 
-       REAL(wp)       ::  average_cputime  !< 
-       REAL(wp), SAVE ::  norm = 1.0_wp    !< 
-       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  pe_max        !< 
-       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  pe_min        !< 
-       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  pe_rms        !< 
-       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  pe_tmp        !< 
-       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  sum           !< 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pe_log_points !< 
+       REAL(dp)       ::  average_cputime  !<
+       REAL(dp), SAVE ::  norm = 1.0_dp    !<
+       REAL(dp), DIMENSION(:),   ALLOCATABLE ::  pe_max        !<
+       REAL(dp), DIMENSION(:),   ALLOCATABLE ::  pe_min        !<
+       REAL(dp), DIMENSION(:),   ALLOCATABLE ::  pe_rms        !<
+       REAL(dp), DIMENSION(:),   ALLOCATABLE ::  pe_tmp        !<
+       REAL(dp), DIMENSION(:),   ALLOCATABLE ::  sum           !<
+       REAL(dp), DIMENSION(:,:), ALLOCATABLE ::  pe_log_points !<
 
 
@@ -358,11 +357,11 @@
 !
 !--       Allocate and initialize temporary arrays needed for statistics
-          ALLOCATE( pe_max( SIZE( log_point ) ), pe_min( SIZE( log_point ) ),  &
-                    pe_rms( SIZE( log_point ) ), pe_tmp( SIZE( log_point ) ),  &
+          ALLOCATE( pe_max( SIZE( log_point ) ), pe_min( SIZE( log_point ) ),                      &
+                    pe_rms( SIZE( log_point ) ), pe_tmp( SIZE( log_point ) ),                      &
                     pe_log_points( SIZE( log_point ), 0:numprocs-1 ) )
           pe_min = log_point%sum
           pe_max = log_point%sum    ! need to be set in case of 1 PE
-          pe_rms = 0.0_wp
-          pe_tmp = 0.0_wp
+          pe_rms = 0.0_dp
+          pe_tmp = 0.0_dp
 
 #if defined( __parallel )
@@ -370,6 +369,6 @@
 !--       Receive data from all PEs
           DO  i = 1, numprocs-1
-             CALL MPI_RECV( pe_tmp(1), SIZE( log_point ), MPI_REAL,            &
-                            i, i, comm2d, status, ierr )
+             CALL MPI_RECV( pe_tmp(1), SIZE( log_point ), MPI_DOUBLE_PRECISION, i, i, comm2d,      &
+                            status, ierr )
              sender = status(MPI_SOURCE)
              pe_log_points(:,sender) = pe_tmp
@@ -389,7 +388,5 @@
 !--          Calculate rms
              DO  i = 0, numprocs-1
-                pe_rms(iii) = pe_rms(iii) + (                                  &
-                                    pe_log_points(iii,i) - log_point(iii)%sum  &
-                                            )**2
+                pe_rms(iii) = pe_rms(iii) + ( pe_log_points(iii,i) - log_point(iii)%sum )**2
              ENDDO
              pe_rms(iii) = SQRT( pe_rms(iii) / numprocs )
@@ -401,6 +398,5 @@
           ALLOCATE( pe_max( SIZE( log_point ) ) )
           pe_max = log_point%sum
-          CALL MPI_SEND( pe_max(1), SIZE( log_point ), MPI_REAL, 0, myid,      &
-                         comm2d, ierr )
+          CALL MPI_SEND( pe_max(1), SIZE( log_point ), MPI_DOUBLE_PRECISION, 0, myid, comm2d, ierr )
 #endif
 
@@ -416,5 +412,5 @@
              sum = log_point%sum
           ELSEWHERE
-             sum = -1.0_wp
+             sum = -1.0_dp
           ENDWHERE
 
@@ -423,9 +419,8 @@
 !--       timestep
           IF ( nr_timesteps_this_run /= 0 )  THEN
-             average_cputime = log_point(1)%sum / REAL( ngp_3d(0), KIND=wp ) / &
-                               REAL( nr_timesteps_this_run, KIND=wp ) * 1E6_wp
-                               ! in micro-sec
+             average_cputime = log_point(1)%sum / REAL( ngp_3d(0), KIND=dp ) /                     &
+                               REAL( nr_timesteps_this_run, KIND=dp ) * 1E6_dp     ! in micro-sec
           ELSE
-             average_cputime = -1.0_wp
+             average_cputime = -1.0_dp
           ENDIF
 
@@ -434,15 +429,13 @@
           CALL check_open( 18 )
 #if defined( __parallel )
-          WRITE ( 18, 100 )  TRIM( run_description_header ),                   &
-                             numprocs * threads_per_task, pdims(1), pdims(2),  &
-                             threads_per_task, nx+1, ny+1, nz,                 &
+          WRITE ( 18, 100 )  TRIM( run_description_header ), numprocs * threads_per_task,          &
+                             pdims(1), pdims(2), threads_per_task, nx+1, ny+1, nz,                 &
                              nr_timesteps_this_run, average_cputime
 
           WRITE ( 18, 110 )
 #else
-          WRITE ( 18, 100 )  TRIM( run_description_header ),                   &
-                             numprocs * threads_per_task, 1, 1,                &
-                             threads_per_task, nx+1, ny+1, nz,                 &
-                             nr_timesteps_this_run, average_cputime
+          WRITE ( 18, 100 )  TRIM( run_description_header ), numprocs * threads_per_task, 1, 1,    &
+                             threads_per_task, nx+1, ny+1, nz, nr_timesteps_this_run,              &
+                             average_cputime
 
           WRITE ( 18, 110 )
@@ -451,10 +444,9 @@
              ii = MAXLOC( sum )
              i = ii(1)
-             IF ( sum(i) /= -1.0_wp )  THEN
-                WRITE ( 18, 102 ) &
-              log_point(i)%place, log_point(i)%sum,                            &
-                   log_point(i)%sum / log_point(1)%sum * 100.0_wp,             &
-                   log_point(i)%counts, pe_min(i), pe_max(i), pe_rms(i)
-                sum(i) = -1.0_wp
+             IF ( sum(i) /= -1.0_dp )  THEN
+                WRITE ( 18, 102 )  log_point(i)%place, log_point(i)%sum,                           &
+                                   log_point(i)%sum / log_point(1)%sum * 100.0_dp,                 &
+                                   log_point(i)%counts, pe_min(i), pe_max(i), pe_rms(i)
+                sum(i) = -1.0_dp
              ELSE
                 EXIT
@@ -488,5 +480,5 @@
           pe_min = log_point_s%sum
           pe_max = log_point_s%sum    ! need to be set in case of 1 PE
-          pe_rms = 0.0_wp
+          pe_rms = 0.0_dp
 
 #if defined( __parallel )
@@ -494,6 +486,6 @@
 !--       Receive data from all PEs
           DO  i = 1, numprocs-1
-             CALL MPI_RECV( pe_tmp(1), SIZE( log_point ), MPI_REAL, &
-                            MPI_ANY_SOURCE, MPI_ANY_TAG, comm2d, status, ierr )
+             CALL MPI_RECV( pe_tmp(1), SIZE( log_point ), MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE,    &
+                            MPI_ANY_TAG, comm2d, status, ierr )
              sender = status(MPI_SOURCE)
              pe_log_points(:,sender) = pe_tmp
@@ -505,6 +497,5 @@
           DO  iii = 1, SIZE( log_point )
              DO  i = 1, numprocs-1
-                log_point_s(iii)%sum = log_point_s(iii)%sum +                  &
-                                       pe_log_points(iii,i)
+                log_point_s(iii)%sum = log_point_s(iii)%sum + pe_log_points(iii,i)
                 pe_min(iii) = MIN( pe_min(iii), pe_log_points(iii,i) )
                 pe_max(iii) = MAX( pe_max(iii), pe_log_points(iii,i) )
@@ -514,7 +505,5 @@
 !--          Calculate rms
              DO  i = 0, numprocs-1
-                pe_rms(iii) = pe_rms(iii) + (                                  &
-                                   pe_log_points(iii,i) - log_point_s(iii)%sum &
-                                            )**2
+                pe_rms(iii) = pe_rms(iii) + ( pe_log_points(iii,i) - log_point_s(iii)%sum )**2
              ENDDO
              pe_rms(iii) = SQRT( pe_rms(iii) / numprocs )
@@ -525,6 +514,5 @@
 !--       the data in order to avoid sending the data type log)
           pe_max = log_point_s%sum
-          CALL MPI_SEND( pe_max(1), SIZE( log_point ), MPI_REAL, 0, 0, comm2d, &
-                         ierr )
+          CALL MPI_SEND( pe_max(1), SIZE( log_point ), MPI_DOUBLE_PRECISION, 0, 0, comm2d, ierr )
 #endif
 
@@ -539,5 +527,5 @@
              sum = log_point_s%sum
           ELSEWHERE
-             sum = -1.0_wp
+             sum = -1.0_dp
           ENDWHERE
 
@@ -548,10 +536,9 @@
              ii = MAXLOC( sum )
              i = ii(1)
-             IF ( sum(i) /= -1.0_wp )  THEN
-                WRITE ( 18, 102 )                                              &
-                   log_point_s(i)%place, log_point_s(i)%sum,                   &
-                   log_point_s(i)%sum / log_point(1)%sum * 100.0_wp,           &
-                   log_point_s(i)%counts, pe_min(i), pe_max(i), pe_rms(i)
-                sum(i) = -1.0_wp
+             IF ( sum(i) /= -1.0_dp )  THEN
+                WRITE ( 18, 102 )  log_point_s(i)%place, log_point_s(i)%sum,                       &
+                                   log_point_s(i)%sum / log_point(1)%sum * 100.0_dp,               &
+                                   log_point_s(i)%counts, pe_min(i), pe_max(i), pe_rms(i)
+                sum(i) = -1.0_dp
              ELSE
                 EXIT
@@ -587,12 +574,12 @@
 
 
-   100 FORMAT (A/11('-')//'CPU measures for ',I5,' PEs (',I5,'(x) * ',I5,'(y', &
-               &') tasks *',I5,' threads):'//                                  &
-               'gridpoints (x/y/z): ',20X,I5,' * ',I5,' * ',I5/                &
-               'nr of timesteps: ',22X,I6/                                     &
+   100 FORMAT (A/11('-')//'CPU measures for ',I5,' PEs (',I5,'(x) * ',I5,'(y',                     &
+               &') tasks *',I5,' threads):'//                                                      &
+               'gridpoints (x/y/z): ',20X,I5,' * ',I5,' * ',I5/                                    &
+               'nr of timesteps: ',22X,I6/                                                         &
                'cpu time per grid point and timestep: ',5X,F8.5,' * 10**-6 s')
 
    101 FORMAT (/'special measures:'/ &
-               &'-----------------------------------------------------------', &
+               &'-----------------------------------------------------------',                     &
                &'---------------------')
 
@@ -603,11 +590,11 @@
    106 FORMAT (/'Exchange of ghostpoints via MPI_ISEND/MPI_IRECV')
    107 FORMAT (//)
-   110 FORMAT ('------------------------------------------------------------', &
-               &'----------'//                                                 &
-               &'place:                              mean        counts     ', &
-               &' min       max    rms'/                                       &
-               &'                                sec.      %                ', &
-               &'sec.      sec.   sec.'/                                       &
-               &'-----------------------------------------------------------', &
+   110 FORMAT ('------------------------------------------------------------',                     &
+               &'----------'//                                                                     &
+               &'place:                              mean        counts     ',                     &
+               &' min       max    rms'/                                                           &
+               &'                                sec.      %                ',                     &
+               &'sec.      sec.   sec.'/                                                           &
+               &'-----------------------------------------------------------',                     &
                &'---------------------')
    111 FORMAT (/'Barriers are set at beginning (start/continue) of measurements')
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 4014)
+++ /palm/trunk/SOURCE/pres.f90	(revision 4015)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! variable child_domain_nvn eliminated
+! 
+! 3849 2019-04-01 16:35:16Z knoop
 ! OpenACC port for SPEC
 ! 
@@ -231,6 +234,4 @@
     REAL(wp), DIMENSION(1:nzt) ::  w_l_l               !<
 
-    LOGICAL :: child_domain_nvn      !<
-
 
     CALL cpu_log( log_point(8), 'pres', 'start' )
@@ -381,18 +382,19 @@
 
 !
-!-- Remove mean vertical velocity in case that Neumann conditions are
-!-- used both at bottom and top boundary, and if not a nested domain in a 
-!-- normal nesting run. In case of vertical nesting, this must be done.
-!-- Therefore an auxiliary logical variable child_domain_nvn is used here, and 
-!-- nvn stands for non-vertical nesting. 
-!-- This cannot be done before the first initial time step because ngp_2dh_outer
+!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
+!-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove
+!-- mean vertical velocities. They should be removed, because incompressibility requires that
+!-- the vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be
+!-- zero everywhere.
+!-- This must not be done in case of a 3d-nesting child domain, because a mean vertical velocity
+!-- can physically exist in such a domain.
+!-- Also in case of offline nesting, mean vertical velocities may exist (and must not be removed),
+!-- caused by horizontal divergence/convergence of the large scale flow that is prescribed at the
+!-- side boundaries.
+!-- The removal cannot be done before the first initial time step because ngp_2dh_outer
 !-- is not yet known then.
-    child_domain_nvn = child_domain
-    IF ( child_domain .AND. nesting_mode == 'vertical' )  THEN
-       child_domain_nvn = .FALSE.
-    ENDIF
-
-    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nesting_offline  .AND. &
-         .NOT. child_domain_nvn  .AND. intermediate_timestep_count /= 0 )       &
+    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nesting_offline                           &
+         .AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' )                             &
+         .AND. intermediate_timestep_count /= 0 )                                                  &
     THEN
        w_l = 0.0_wp;  w_l_l = 0.0_wp
@@ -400,8 +402,6 @@
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
-                w_l_l(k) = w_l_l(k) + w(k,j,i)                                 &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 3 )  &
-                                            )
+                w_l_l(k) = w_l_l(k) + w(k,j,i)                                                     &
+                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3 ) )
              ENDDO
           ENDDO
@@ -409,6 +409,5 @@
 #if defined( __parallel )   
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
-                           comm2d, ierr )
+       CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, ierr )
 #else
        w_l = w_l_l
@@ -420,8 +419,6 @@
           DO  j = nysg, nyng
              DO  k = nzb+1, nzt
-                w(k,j,i) = w(k,j,i) - w_l(k)                                   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 3 )  &
-                                            )
+                w(k,j,i) = w(k,j,i) - w_l(k)                                                       &
+                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3 ) )
              ENDDO
           ENDDO
