Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 1375)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 1557)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Enable monotone advection for scalars using monotonic limiter 
 ! 
 ! Former revisions:
@@ -194,6 +194,6 @@
        USE control_parameters,                                                 &
            ONLY:  cloud_physics, humidity, icloud_scheme, loop_optimization,   &
-                  passive_scalar, precipitation, ocean, ws_scheme_mom,         &
-                  ws_scheme_sca
+                  monotonic_adjustment, passive_scalar, precipitation, ocean,  &
+                  ws_scheme_mom, ws_scheme_sca
 
        USE indices,                                                            &
@@ -386,5 +386,6 @@
 
        USE control_parameters,                                                 &
-           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans 
+           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans, &
+                  v_gtrans 
 
        USE grid_variables,                                                     &
@@ -421,4 +422,5 @@
        INTEGER(iwp) ::  k     !:
        INTEGER(iwp) ::  k_mm  !:
+       INTEGER(iwp) ::  k_mmm !:
        INTEGER(iwp) ::  k_pp  !:
        INTEGER(iwp) ::  k_ppp !:
@@ -428,4 +430,22 @@
        REAL(wp)     ::  div    !:
        REAL(wp)     ::  flux_d !:
+       REAL(wp)     ::  fd_1   !:
+       REAL(wp)     ::  fl_1   !:
+       REAL(wp)     ::  fn_1   !:
+       REAL(wp)     ::  fr_1   !:
+       REAL(wp)     ::  fs_1   !:
+       REAL(wp)     ::  ft_1   !:
+       REAL(wp)     ::  phi_d  !:
+       REAL(wp)     ::  phi_l  !:
+       REAL(wp)     ::  phi_n  !:
+       REAL(wp)     ::  phi_r  !:
+       REAL(wp)     ::  phi_s  !:
+       REAL(wp)     ::  phi_t  !:
+       REAL(wp)     ::  rd     !:
+       REAL(wp)     ::  rl     !:
+       REAL(wp)     ::  rn     !:
+       REAL(wp)     ::  rr     !:
+       REAL(wp)     ::  rs     !:
+       REAL(wp)     ::  rt     !:
        REAL(wp)     ::  u_comp !:
        REAL(wp)     ::  v_comp !:
@@ -697,4 +717,142 @@
                                          )
 !
+!--       Apply monotonic adjustment.
+          IF ( monotonic_adjustment )  THEN
+!
+!--          At first, calculate first order fluxes.
+             u_comp = u(k,j,i) - u_gtrans
+             fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )               &
+                   -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )               &
+                       ) * adv_sca_1 
+
+             u_comp = u(k,j,i+1) - u_gtrans
+             fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )               &
+                   -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )               &
+                       ) * adv_sca_1 
+
+             v_comp = v(k,j,i) - v_gtrans
+             fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )               &
+                   -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )               &
+                       ) * adv_sca_1 
+
+             v_comp = v(k,j+1,i) - v_gtrans
+             fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )               &
+                   -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )               &
+                       ) * adv_sca_1 
+
+             fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
+                   -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
+                       ) * adv_sca_1 
+
+             ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
+                   -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
+                      ) * adv_sca_1 
+!
+!--          Calculate ratio of upwind gradients. Note, Min/Max is just to 
+!--          avoid if statements.
+             rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *                     & 
+                           ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, u(k,j,i) - u_gtrans ) *                     &
+                           ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /    &
+                                ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
+
+             rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *                   & 
+                           ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /    &
+                                ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, u(k,j,i+1) - u_gtrans ) *                   &
+                           ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
+
+             rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *                     & 
+                           ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, v(k,j,i) - v_gtrans ) *                     &
+                           ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /    &
+                                ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
+
+             rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *                   & 
+                           ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /    &
+                                ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )      &
+                              ) +                                             & 
+                       MIN( 0.0, v(k,j+1,i) - v_gtrans ) *                    &
+                           ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
+!   
+!--          Reuse k_mm and compute k_mmm for the vertical gradient ratios. 
+!--          Note, for vertical advection below the third grid point above 
+!--          surface ( or below the model top) rd and rt are set to 0, i.e. 
+!--          use of first order scheme is enforced.
+             k_mmm  = k - 3 * ibit8
+
+             rd     = ( MAX( 0.0, w(k-1,j,i) ) *                              & 
+                           ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) /  &
+                                ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )    &
+                              ) +                                             & 
+                        MIN( 0.0, w(k-1,j,i) ) *                              &
+                           ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /      &
+                                ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )   &
+                              )                                               &
+                      ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 
+ 
+             rt     = ( MAX( 0.0, w(k,j,i) ) *                                & 
+                           ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /    &
+                                ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, w(k,j,i) ) *                                &
+                           ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /    &
+                                ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )      &
+                              )                                               &
+                      ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
+!
+!--           Calculate empirical limiter function (van Albada2 limiter).
+              phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                          &
+                                         ( rl**2 + 1.0_wp ) ) 
+              phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                          &
+                                         ( rr**2 + 1.0_wp ) ) 
+              phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                          &
+                                         ( rs**2 + 1.0_wp ) ) 
+              phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                          &
+                                         ( rn**2 + 1.0_wp ) ) 
+              phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                          &
+                                         ( rd**2 + 1.0_wp ) ) 
+              phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                          &
+                                         ( rt**2 + 1.0_wp ) ) 
+!
+!--           Calculate the resulting monotone flux. 
+              swap_flux_x_local(k,j,tn) = fl_1 - phi_l *                      &
+                                        ( fl_1 - swap_flux_x_local(k,j,tn)  )
+              flux_r(k)                 = fr_1 - phi_r *                      &
+                                        ( fr_1 - flux_r(k)                  )
+              swap_flux_y_local(k,tn)   = fs_1 - phi_s *                      &
+                                        ( fs_1 - swap_flux_y_local(k,tn)    )
+              flux_n(k)                 = fn_1 - phi_n *                      &
+                                        ( fn_1 - flux_n(k)                  )
+              flux_d                    = fd_1  !- phi_d *                      &
+!                                         ( fd_1 - flux_d                     )
+              flux_t(k)                 = ft_1 !- phi_t *                      &
+!                                         ( ft_1 - flux_t(k)                  )
+!
+!--          Moreover, modify dissipation flux according to the limiter. 
+             swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l
+             diss_r(k)                 = diss_r(k)                 * phi_r
+             swap_diss_y_local(k,tn)   = swap_diss_y_local(k,tn)   * phi_s
+             diss_n(k)                 = diss_n(k)                 * phi_n
+             diss_d                    = 0.0 !diss_d                    * phi_d
+             diss_t(k)                 = 0.0 !diss_t(k)                 * phi_t
+
+          ENDIF
+!
 !--       Calculate the divergence of the velocity field. A respective
 !--       correction is needed to overcome numerical instabilities caused
@@ -756,5 +914,4 @@
           k_pp  = k + 2 * ( 1 - ibit6  )
           k_mm  = k - 2 * ibit8
-
 
           flux_t(k) = w(k,j,i) * (                                            &
@@ -786,4 +943,144 @@
                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
                                          )
+
+
+!
+!--       Apply monotonic adjustment.
+          IF ( monotonic_adjustment )  THEN
+!
+!--          At first, calculate first order fluxes.
+             u_comp = u(k,j,i) - u_gtrans
+             fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )               &
+                   -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )               &
+                       ) * adv_sca_1 
+
+             u_comp = u(k,j,i+1) - u_gtrans
+             fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )               &
+                   -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )               &
+                       ) * adv_sca_1 
+
+             v_comp = v(k,j,i) - v_gtrans
+             fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )               &
+                   -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )               &
+                       ) * adv_sca_1 
+
+             v_comp = v(k,j+1,i) - v_gtrans
+             fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )               &
+                   -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )               &
+                       ) * adv_sca_1 
+
+             fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
+                   -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
+                       ) * adv_sca_1 
+
+             ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
+                   -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
+                       ) * adv_sca_1 
+!
+!--          Calculate ratio of upwind gradients. Note, Min/Max is just to 
+!--          avoid if statements. 
+             rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *                     & 
+                           ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, u(k,j,i) - u_gtrans ) *                     &
+                           ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /    &
+                                ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
+
+             rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *                   & 
+                           ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /    &
+                                ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, u(k,j,i+1) - u_gtrans ) *                   &
+                           ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
+
+             rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *                     & 
+                           ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, v(k,j,i) - v_gtrans ) *                     &
+                           ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /    &
+                                ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
+
+             rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *                   & 
+                           ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /    &
+                                ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )      &
+                              ) +                                             & 
+                       MIN( 0.0, v(k,j+1,i) - v_gtrans ) *                    &
+                           ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /    &
+                                ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )      &
+                              )                                               &
+                      ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
+!
+!--          Reuse k_mm and compute k_mmm for the vertical gradient ratios. 
+!--          Note, for vertical advection below the third grid point above 
+!--          surface ( or below the model top) rd and rt are set to 0, i.e. 
+!--          use of first order scheme is enforced.
+             k_mmm  = k - 3 * ibit8
+
+             rd     = ( MAX( 0.0, w(k-1,j,i) ) *                              & 
+                           ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) /  &
+                                ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )    &
+                              ) +                                             & 
+                        MIN( 0.0, w(k-1,j,i) ) *                              &
+                           ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /      &
+                                ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )   &
+                              )                                               &
+                      ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 
+ 
+             rt     = ( MAX( 0.0, w(k,j,i) ) *                                & 
+                           ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /    &
+                                ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )      &
+                              ) +                                             & 
+                        MIN( 0.0, w(k,j,i) ) *                                &
+                           ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /    &
+                                ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )      &
+                              )                                               &
+                      ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 
+!
+!--           Calculate empirical limiter function (van Albada2 limiter).
+              phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                          &
+                                         ( rl**2 + 1.0_wp ) ) 
+              phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                          &
+                                         ( rr**2 + 1.0_wp ) ) 
+              phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                          &
+                                         ( rs**2 + 1.0_wp ) ) 
+              phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                          &
+                                         ( rn**2 + 1.0_wp ) ) 
+              phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                          &
+                                         ( rd**2 + 1.0_wp ) ) 
+              phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                          &
+                                         ( rt**2 + 1.0_wp ) ) 
+!
+!--           Calculate the resulting monotone flux. 
+              swap_flux_x_local(k,j,tn) = fl_1 - phi_l *                      &
+                                        ( fl_1 - swap_flux_x_local(k,j,tn)  )
+              flux_r(k)                 = fr_1 - phi_r *                      &
+                                        ( fr_1 - flux_r(k)                  )
+              swap_flux_y_local(k,tn)   = fs_1 - phi_s *                      &
+                                        ( fs_1 - swap_flux_y_local(k,tn)    )
+              flux_n(k)                 = fn_1 - phi_n *                      &
+                                        ( fn_1 - flux_n(k)                  )
+              flux_d                    = fd_1 - phi_d *                      &
+                                        ( fd_1 - flux_d                     )
+              flux_t(k)                 = ft_1 - phi_t *                      &
+                                        ( ft_1 - flux_t(k)                  )
+!
+!--          Moreover, modify dissipation flux according to the limiter. 
+             swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l
+             diss_r(k)                 = diss_r(k)                 * phi_r
+             swap_diss_y_local(k,tn)   = swap_diss_y_local(k,tn)   * phi_s
+             diss_n(k)                 = diss_n(k)                 * phi_n
+             diss_d                    = diss_d                    * phi_d
+             diss_t(k)                 = diss_t(k)                 * phi_t
+
+          ENDIF
 !
 !--       Calculate the divergence of the velocity field. A respective
@@ -803,4 +1100,5 @@
                                       ) + sk(k,j,i) * div
 
+
           swap_flux_y_local(k,tn)   = flux_n(k)
           swap_diss_y_local(k,tn)   = diss_n(k)
@@ -813,48 +1111,52 @@
 
 !
-!--    Evaluation of statistics
-       SELECT CASE ( sk_char )
-
-          CASE ( 'pt' )
-
-             DO  k = nzb, nzt
-                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +               &
-                                       ( flux_t(k) + diss_t(k) )              &
-                                 * weight_substep(intermediate_timestep_count)
+!--    Evaluation of statistics. Please note, in case of using monotone limiter
+!--    vertical fluxes will be not calculated by the advective fluxes.
+       IF ( .NOT. monotonic_adjustment )  THEN
+
+          SELECT CASE ( sk_char )
+
+             CASE ( 'pt' )
+
+                DO  k = nzb, nzt
+                   sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +            &
+                                          ( flux_t(k) + diss_t(k) )           &
+                                  * weight_substep(intermediate_timestep_count)
              ENDDO
              
-          CASE ( 'sa' )
-
-             DO  k = nzb, nzt
-                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +               &
-                                       ( flux_t(k) + diss_t(k) )              &
-                                 * weight_substep(intermediate_timestep_count)
-             ENDDO
+             CASE ( 'sa' )
+
+                DO  k = nzb, nzt
+                   sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +            &
+                                          ( flux_t(k) + diss_t(k) )           &
+                                  * weight_substep(intermediate_timestep_count)
+                ENDDO
              
-          CASE ( 'q' )
-
-             DO  k = nzb, nzt
-                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                &
-                                      ( flux_t(k) + diss_t(k) )               &
-                                 * weight_substep(intermediate_timestep_count)
-             ENDDO
-
-          CASE ( 'qr' )
-
-             DO  k = nzb, nzt
-                sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +              &
-                                      ( flux_t(k) + diss_t(k) )               &
-                                 * weight_substep(intermediate_timestep_count)
-             ENDDO
-
-          CASE ( 'nr' )
-
-             DO  k = nzb, nzt
-                sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +              &
-                                      ( flux_t(k) + diss_t(k) )               &
-                                 * weight_substep(intermediate_timestep_count)
-             ENDDO
-
-         END SELECT
+             CASE ( 'q' )
+
+                DO  k = nzb, nzt
+                   sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +             &
+                                         ( flux_t(k) + diss_t(k) )            &
+                                  * weight_substep(intermediate_timestep_count)
+                ENDDO
+
+             CASE ( 'qr' )
+
+                DO  k = nzb, nzt
+                   sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +           &
+                                         ( flux_t(k) + diss_t(k) )            &
+                                  * weight_substep(intermediate_timestep_count)
+                ENDDO
+
+             CASE ( 'nr' )
+
+                DO  k = nzb, nzt
+                   sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +           &
+                                         ( flux_t(k) + diss_t(k) )            &
+                                  * weight_substep(intermediate_timestep_count)
+                ENDDO
+
+            END SELECT
+         ENDIF
 
     END SUBROUTINE advec_s_ws_ij
@@ -2228,5 +2530,6 @@
 
        USE control_parameters,                                                &
-           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans 
+           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
+                  v_gtrans 
 
        USE grid_variables,                                                    &
@@ -2247,473 +2550,4 @@
        CHARACTER (LEN = *), INTENT(IN)    ::  sk_char !:
        
-       INTEGER(iwp) ::  i      !:
-       INTEGER(iwp) ::  ibit0  !:
-       INTEGER(iwp) ::  ibit1  !:
-       INTEGER(iwp) ::  ibit2  !:
-       INTEGER(iwp) ::  ibit3  !:
-       INTEGER(iwp) ::  ibit4  !:
-       INTEGER(iwp) ::  ibit5  !:
-       INTEGER(iwp) ::  ibit6  !:
-       INTEGER(iwp) ::  ibit7  !:
-       INTEGER(iwp) ::  ibit8  !:
-       INTEGER(iwp) ::  j      !:
-       INTEGER(iwp) ::  k      !:
-       INTEGER(iwp) ::  k_mm   !:
-       INTEGER(iwp) ::  k_pp   !:
-       INTEGER(iwp) ::  k_ppp  !:
-       INTEGER(iwp) ::  tn = 0 !:
-       
-#if defined( __nopointer )
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !:
-#else
-       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !:
-#endif
-
-       REAL(wp) ::  diss_d !:
-       REAL(wp) ::  div    !:
-       REAL(wp) ::  flux_d !:
-       REAL(wp) ::  u_comp !:
-       REAL(wp) ::  v_comp !:
-       
-       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_n !:
-       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_r !:
-       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_t !:
-       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_n !:
-       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_r !:
-       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_t !:
-       
-       REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local !:
-       REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local !:
-       
-       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !:
-       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !:
-        
-
-!
-!--    Compute the fluxes for the whole left boundary of the processor domain.
-       i = nxl
-       DO  j = nys, nyn
-
-          DO  k = nzb+1, nzb_max
-
-             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
-             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
-             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
-
-             u_comp                 = u(k,j,i) - u_gtrans
-             swap_flux_x_local(k,j) = u_comp * (                              &
-                                             ( 37.0_wp * ibit2 * adv_sca_5    &
-                                          +     7.0_wp * ibit1 * adv_sca_3    &
-                                          +              ibit0 * adv_sca_1    &
-                                             ) *                              &
-                                          ( sk(k,j,i)   + sk(k,j,i-1)    )    &
-                                      -      (  8.0_wp * ibit2 * adv_sca_5    &
-                                          +              ibit1 * adv_sca_3    &
-                                             ) *                              &
-                                          ( sk(k,j,i+1) + sk(k,j,i-2)    )    &
-                                      +      (           ibit2 * adv_sca_5    & 
-                                             ) *                              &
-                                          ( sk(k,j,i+2) + sk(k,j,i-3)    )    &
-                                               )
-
-              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
-                                             ( 10.0_wp * ibit2 * adv_sca_5    &
-                                          +     3.0_wp * ibit1 * adv_sca_3    &
-                                          +              ibit0 * adv_sca_1    &
-                                             ) *                              &
-                                          ( sk(k,j,i)   - sk(k,j,i-1) )       &
-                                      -      (  5.0_wp * ibit2 * adv_sca_5    &
-                                          +              ibit1 * adv_sca_3    &
-                                             ) *                              &
-                                         ( sk(k,j,i+1) - sk(k,j,i-2)  )       &
-                                      +      (           ibit2 * adv_sca_5    &
-                                             ) *                              &
-                                          ( sk(k,j,i+2) - sk(k,j,i-3) )       &
-                                                        )
-
-          ENDDO
-
-          DO  k = nzb_max+1, nzt
-
-             u_comp                 = u(k,j,i) - u_gtrans
-             swap_flux_x_local(k,j) = u_comp * (                              &
-                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
-                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
-                                    +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
-                                               ) * adv_sca_5
-
-             swap_diss_x_local(k,j) = -ABS( u_comp ) * (                      &
-                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
-                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
-                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
-                                                       ) * adv_sca_5
-
-          ENDDO
-
-       ENDDO
-
-       DO  i = nxl, nxr
-
-          j = nys
-          DO  k = nzb+1, nzb_max
-
-             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
-             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
-             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
-
-             v_comp               = v(k,j,i) - v_gtrans
-             swap_flux_y_local(k) = v_comp * (                                &
-                                             ( 37.0_wp * ibit5 * adv_sca_5    &
-                                          +     7.0_wp * ibit4 * adv_sca_3    &
-                                          +              ibit3 * adv_sca_1    &
-                                             ) *                              &
-                                         ( sk(k,j,i)  + sk(k,j-1,i)     )     &
-                                       -     (  8.0_wp * ibit5 * adv_sca_5    &
-                                          +              ibit4 * adv_sca_3    &
-                                              ) *                             &
-                                         ( sk(k,j+1,i) + sk(k,j-2,i)    )     &
-                                      +      (           ibit5 * adv_sca_5    &
-                                             ) *                              &
-                                        ( sk(k,j+2,i) + sk(k,j-3,i)     )     &
-                                             )
-
-             swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
-                                             ( 10.0_wp * ibit5 * adv_sca_5    &
-                                          +     3.0_wp * ibit4 * adv_sca_3    &
-                                          +              ibit3 * adv_sca_1    &
-                                             ) *                              &
-                                          ( sk(k,j,i)   - sk(k,j-1,i)    )    &
-                                      -      (  5.0_wp * ibit5 * adv_sca_5    &
-                                          +              ibit4 * adv_sca_3    &
-                                             ) *                              &
-                                          ( sk(k,j+1,i) - sk(k,j-2,i)    )    &
-                                      +      (           ibit5 * adv_sca_5    &
-                                             ) *                              &
-                                          ( sk(k,j+2,i) - sk(k,j-3,i)    )    &
-                                                     )
-
-          ENDDO
-!
-!--       Above to the top of the highest topography. No degradation necessary.
-          DO  k = nzb_max+1, nzt
-
-             v_comp               = v(k,j,i) - v_gtrans
-             swap_flux_y_local(k) = v_comp * (                               &
-                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
-                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )  &
-                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
-                                             ) * adv_sca_5
-              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
-                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
-                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
-                                  +             sk(k,j+2,i) - sk(k,j-3,i)    &
-                                                      ) * adv_sca_5
-
-          ENDDO
-
-          DO  j = nys, nyn
-
-             flux_t(0) = 0.0_wp
-             diss_t(0) = 0.0_wp
-             flux_d    = 0.0_wp
-             diss_d    = 0.0_wp
-
-             DO  k = nzb+1, nzb_max
-
-                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
-                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
-                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
-
-                u_comp    = u(k,j,i+1) - u_gtrans
-                flux_r(k) = u_comp * (                                        &
-                          ( 37.0_wp * ibit2 * adv_sca_5                       &
-                      +      7.0_wp * ibit1 * adv_sca_3                       &
-                      +               ibit0 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+1) + sk(k,j,i)   )                    &
-                   -      (  8.0_wp * ibit2 * adv_sca_5                       &
-                       +              ibit1 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
-                   +      (           ibit2 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
-                                     )
-
-                diss_r(k) = -ABS( u_comp ) * (                                &
-                          ( 10.0_wp * ibit2 * adv_sca_5                       &
-                       +     3.0_wp * ibit1 * adv_sca_3                       &
-                       +              ibit0 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+1) - sk(k,j,i)   )                    &
-                   -      (  5.0_wp * ibit2 * adv_sca_5                       &
-                       +              ibit1 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
-                   +      (           ibit2 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
-                                             )
-
-                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
-                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
-                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
-
-                v_comp    = v(k,j+1,i) - v_gtrans
-                flux_n(k) = v_comp * (                                        &
-                          ( 37.0_wp * ibit5 * adv_sca_5                       &
-                       +     7.0_wp * ibit4 * adv_sca_3                       &
-                       +              ibit3 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j+1,i) + sk(k,j,i)   )                    &
-                   -      (  8.0_wp * ibit5 * adv_sca_5                       &
-                       +              ibit4 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
-                   +      (           ibit5 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
-                                     )
-
-                diss_n(k) = -ABS( v_comp ) * (                                &
-                          ( 10.0_wp * ibit5 * adv_sca_5                       &
-                       +     3.0_wp * ibit4 * adv_sca_3                       &
-                       +              ibit3 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j+1,i) - sk(k,j,i)    )                   &
-                   -      (  5.0_wp * ibit5 * adv_sca_5                       &
-                       +              ibit4 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                   &
-                   +      (           ibit5 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
-                                             )
-!
-!--             k index has to be modified near bottom and top, else array
-!--             subscripts will be exceeded.
-                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
-                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
-                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
-
-                k_ppp = k + 3 * ibit8
-                k_pp  = k + 2 * ( 1 - ibit6  )
-                k_mm  = k - 2 * ibit8
-
-
-                flux_t(k) = w(k,j,i) * (                                      &
-                           ( 37.0_wp * ibit8 * adv_sca_5                      &
-                        +     7.0_wp * ibit7 * adv_sca_3                      &
-                        +           ibit6 * adv_sca_1                         &
-                           ) *                                                &
-                                   ( sk(k+1,j,i)  + sk(k,j,i)    )            &
-                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
-                        +              ibit7 * adv_sca_3                      &
-                           ) *                                                &
-                                   ( sk(k_pp,j,i) + sk(k-1,j,i)  )            &
-                    +      (           ibit8 * adv_sca_5                      &
-                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
-                                       )
-
-                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
-                           ( 10.0_wp * ibit8 * adv_sca_5                      &
-                        +     3.0_wp * ibit7 * adv_sca_3                      &
-                        +              ibit6 * adv_sca_1                      &
-                           ) *                                                &
-                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
-                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
-                        +              ibit7 * adv_sca_3                      &
-                           ) *                                                &
-                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
-                    +      (           ibit8 * adv_sca_5                      &
-                           ) *                                                &
-                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
-                                               )
-!
-!--             Calculate the divergence of the velocity field. A respective
-!--             correction is needed to overcome numerical instabilities caused
-!--             by a not sufficient reduction of divergences near topography.
-                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
-                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
-                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
-
-                tend(k,j,i) = tend(k,j,i) - (                                 &
-                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
-                          swap_diss_x_local(k,j)            ) * ddx           &
-                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
-                          swap_diss_y_local(k)              ) * ddy           &
-                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
-                                                               ) * ddzw(k)    &
-                                            ) + sk(k,j,i) * div
-
-                swap_flux_y_local(k)   = flux_n(k)
-                swap_diss_y_local(k)   = diss_n(k)
-                swap_flux_x_local(k,j) = flux_r(k)
-                swap_diss_x_local(k,j) = diss_r(k)
-                flux_d                 = flux_t(k)
-                diss_d                 = diss_t(k)
-
-             ENDDO
-
-             DO  k = nzb_max+1, nzt
-
-                u_comp    = u(k,j,i+1) - u_gtrans
-                flux_r(k) = u_comp * (                                        &
-                      37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
-                    -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
-                    +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
-                diss_r(k) = -ABS( u_comp ) * (                                &
-                      10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
-                    -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
-                    +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
-
-                v_comp    = v(k,j+1,i) - v_gtrans
-                flux_n(k) = v_comp * (                                        &
-                      37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
-                    -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
-                    +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
-                diss_n(k) = -ABS( v_comp ) * (                                &
-                      10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
-                    -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
-                    +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
-!
-!--             k index has to be modified near bottom and top, else array
-!--             subscripts will be exceeded.
-                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
-                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
-                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
-
-                k_ppp = k + 3 * ibit8
-                k_pp  = k + 2 * ( 1 - ibit6  )
-                k_mm  = k - 2 * ibit8
-
-
-                flux_t(k) = w(k,j,i) * (                                      &
-                           ( 37.0_wp * ibit8 * adv_sca_5                      &
-                        +     7.0_wp * ibit7 * adv_sca_3                      &
-                        +              ibit6 * adv_sca_1                      &
-                           ) *                                                &
-                                   ( sk(k+1,j,i)  + sk(k,j,i)     )           &
-                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
-                        +              ibit7 * adv_sca_3                      &
-                           ) *                                                &
-                                   ( sk(k_pp,j,i) + sk(k-1,j,i)   )           &
-                    +      (           ibit8 * adv_sca_5                      &
-                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i)  )           &
-                                       )
-
-                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
-                           ( 10.0_wp * ibit8 * adv_sca_5                      &
-                        +     3.0_wp * ibit7 * adv_sca_3                      &
-                        +              ibit6 * adv_sca_1                      &
-                           ) *                                                &
-                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
-                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
-                        +              ibit7 * adv_sca_3                      &
-                           ) *                                                &
-                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
-                    +      (           ibit8 * adv_sca_5                      &
-                           ) *                                                &
-                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
-                                               )
-!
-!--             Calculate the divergence of the velocity field. A respective
-!--             correction is needed to overcome numerical instabilities introduced
-!--             by a not sufficient reduction of divergences near topography.
-                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
-                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
-                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
-
-                tend(k,j,i) = tend(k,j,i) - (                                 &
-                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
-                          swap_diss_x_local(k,j)            ) * ddx           &
-                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
-                          swap_diss_y_local(k)              ) * ddy           &
-                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
-                                                               ) * ddzw(k)    &
-                                            ) + sk(k,j,i) * div
-
-                swap_flux_y_local(k)   = flux_n(k)
-                swap_diss_y_local(k)   = diss_n(k)
-                swap_flux_x_local(k,j) = flux_r(k)
-                swap_diss_x_local(k,j) = diss_r(k)
-                flux_d                 = flux_t(k)
-                diss_d                 = diss_t(k)
-
-             ENDDO
-!
-!--          evaluation of statistics
-             SELECT CASE ( sk_char )
-
-                 CASE ( 'pt' )
-                    DO  k = nzb, nzt
-                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
-                        + ( flux_t(k) + diss_t(k) )                          &
-                        *   weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'sa' )
-                    DO  k = nzb, nzt
-                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
-                        + ( flux_t(k) + diss_t(k) )                          &
-                        *   weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'q' )
-                    DO  k = nzb, nzt
-                       sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
-                        + ( flux_t(k) + diss_t(k) )                          &
-                        *   weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'qr' )
-                    DO  k = nzb, nzt
-                       sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
-                        + ( flux_t(k) + diss_t(k) )                          &
-                        *   weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'nr' )
-                    DO  k = nzb, nzt
-                       sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
-                        + ( flux_t(k) + diss_t(k) )                          &
-                        *   weight_substep(intermediate_timestep_count)
-                    ENDDO
-
-              END SELECT
-
-         ENDDO
-      ENDDO
-
-    END SUBROUTINE advec_s_ws
-
-
-!------------------------------------------------------------------------------!
-! Scalar advection - Call for all grid points - accelerator version
-!------------------------------------------------------------------------------!
-    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
-
-       USE arrays_3d,                                                         &
-           ONLY:  ddzw, tend, u, v, w
-
-       USE constants,                                                         &
-           ONLY:  adv_sca_1, adv_sca_3, adv_sca_5
-
-       USE control_parameters,                                                &
-           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
-
-       USE grid_variables,                                                    &
-           ONLY:  ddx, ddy
-
-       USE indices,                                                           &
-           ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,  &
-                  nzb, nzb_max, nzt, wall_flags_0
-
-       USE kinds
-       
-!        USE statistics,                                                       &
-!            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
-!                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
-
-       IMPLICIT NONE
-
-       CHARACTER (LEN = *), INTENT(IN)    :: sk_char !:
-
        INTEGER(iwp) ::  i      !:
        INTEGER(iwp) ::  ibit0  !:
@@ -2733,4 +2567,775 @@
        INTEGER(iwp) ::  k_ppp  !:
        INTEGER(iwp) ::  tn = 0 !:
+       
+#if defined( __nopointer )
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !:
+#else
+       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !:
+#endif
+
+       REAL(wp) ::  diss_d !:
+       REAL(wp) ::  div    !:
+       REAL(wp) ::  flux_d !:
+       REAL(wp) ::  fd_1   !:
+       REAL(wp) ::  fl_1   !:
+       REAL(wp) ::  fn_1   !:
+       REAL(wp) ::  fr_1   !:
+       REAL(wp) ::  fs_1   !:
+       REAL(wp) ::  ft_1   !:
+       REAL(wp) ::  phi_d  !:
+       REAL(wp) ::  phi_l  !:
+       REAL(wp) ::  phi_n  !:
+       REAL(wp) ::  phi_r  !:
+       REAL(wp) ::  phi_s  !:
+       REAL(wp) ::  phi_t  !:
+       REAL(wp) ::  rd     !:
+       REAL(wp) ::  rl     !:
+       REAL(wp) ::  rn     !:
+       REAL(wp) ::  rr     !:
+       REAL(wp) ::  rs     !:
+       REAL(wp) ::  rt     !:
+       REAL(wp) ::  u_comp !:
+       REAL(wp) ::  v_comp !:
+       
+       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_n !:
+       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_r !:
+       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_t !:
+       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_n !:
+       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_r !:
+       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_t !:
+       
+       REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local !:
+       REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local !:
+       
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !:
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !:
+        
+
+!
+!--    Compute the fluxes for the whole left boundary of the processor domain.
+       i = nxl
+       DO  j = nys, nyn
+
+          DO  k = nzb+1, nzb_max
+
+             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
+             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
+             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
+
+             u_comp                 = u(k,j,i) - u_gtrans
+             swap_flux_x_local(k,j) = u_comp * (                              &
+                                             ( 37.0_wp * ibit2 * adv_sca_5    &
+                                          +     7.0_wp * ibit1 * adv_sca_3    &
+                                          +              ibit0 * adv_sca_1    &
+                                             ) *                              &
+                                          ( sk(k,j,i)   + sk(k,j,i-1)    )    &
+                                      -      (  8.0_wp * ibit2 * adv_sca_5    &
+                                          +              ibit1 * adv_sca_3    &
+                                             ) *                              &
+                                          ( sk(k,j,i+1) + sk(k,j,i-2)    )    &
+                                      +      (           ibit2 * adv_sca_5    & 
+                                             ) *                              &
+                                          ( sk(k,j,i+2) + sk(k,j,i-3)    )    &
+                                               )
+
+              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
+                                             ( 10.0_wp * ibit2 * adv_sca_5    &
+                                          +     3.0_wp * ibit1 * adv_sca_3    &
+                                          +              ibit0 * adv_sca_1    &
+                                             ) *                              &
+                                          ( sk(k,j,i)   - sk(k,j,i-1) )       &
+                                      -      (  5.0_wp * ibit2 * adv_sca_5    &
+                                          +              ibit1 * adv_sca_3    &
+                                             ) *                              &
+                                         ( sk(k,j,i+1) - sk(k,j,i-2)  )       &
+                                      +      (           ibit2 * adv_sca_5    &
+                                             ) *                              &
+                                          ( sk(k,j,i+2) - sk(k,j,i-3) )       &
+                                                        )
+
+          ENDDO
+
+          DO  k = nzb_max+1, nzt
+
+             u_comp                 = u(k,j,i) - u_gtrans
+             swap_flux_x_local(k,j) = u_comp * (                              &
+                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
+                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
+                                    +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
+                                               ) * adv_sca_5
+
+             swap_diss_x_local(k,j) = -ABS( u_comp ) * (                      &
+                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
+                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
+                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
+                                                       ) * adv_sca_5
+
+          ENDDO
+
+       ENDDO
+
+       DO  i = nxl, nxr
+
+          j = nys
+          DO  k = nzb+1, nzb_max
+
+             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
+             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
+             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
+
+             v_comp               = v(k,j,i) - v_gtrans
+             swap_flux_y_local(k) = v_comp * (                                &
+                                             ( 37.0_wp * ibit5 * adv_sca_5    &
+                                          +     7.0_wp * ibit4 * adv_sca_3    &
+                                          +              ibit3 * adv_sca_1    &
+                                             ) *                              &
+                                         ( sk(k,j,i)  + sk(k,j-1,i)     )     &
+                                       -     (  8.0_wp * ibit5 * adv_sca_5    &
+                                          +              ibit4 * adv_sca_3    &
+                                              ) *                             &
+                                         ( sk(k,j+1,i) + sk(k,j-2,i)    )     &
+                                      +      (           ibit5 * adv_sca_5    &
+                                             ) *                              &
+                                        ( sk(k,j+2,i) + sk(k,j-3,i)     )     &
+                                             )
+
+             swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
+                                             ( 10.0_wp * ibit5 * adv_sca_5    &
+                                          +     3.0_wp * ibit4 * adv_sca_3    &
+                                          +              ibit3 * adv_sca_1    &
+                                             ) *                              &
+                                          ( sk(k,j,i)   - sk(k,j-1,i)    )    &
+                                      -      (  5.0_wp * ibit5 * adv_sca_5    &
+                                          +              ibit4 * adv_sca_3    &
+                                             ) *                              &
+                                          ( sk(k,j+1,i) - sk(k,j-2,i)    )    &
+                                      +      (           ibit5 * adv_sca_5    &
+                                             ) *                              &
+                                          ( sk(k,j+2,i) - sk(k,j-3,i)    )    &
+                                                     )
+
+          ENDDO
+!
+!--       Above to the top of the highest topography. No degradation necessary.
+          DO  k = nzb_max+1, nzt
+
+             v_comp               = v(k,j,i) - v_gtrans
+             swap_flux_y_local(k) = v_comp * (                               &
+                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
+                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )  &
+                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
+                                             ) * adv_sca_5
+              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
+                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
+                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
+                                  +             sk(k,j+2,i) - sk(k,j-3,i)    &
+                                                      ) * adv_sca_5
+
+          ENDDO
+
+          DO  j = nys, nyn
+
+             flux_t(0) = 0.0_wp
+             diss_t(0) = 0.0_wp
+             flux_d    = 0.0_wp
+             diss_d    = 0.0_wp
+
+             DO  k = nzb+1, nzb_max
+
+                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
+                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
+                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
+
+                u_comp    = u(k,j,i+1) - u_gtrans
+                flux_r(k) = u_comp * (                                        &
+                          ( 37.0_wp * ibit2 * adv_sca_5                       &
+                      +      7.0_wp * ibit1 * adv_sca_3                       &
+                      +               ibit0 * adv_sca_1                       &
+                          ) *                                                 &
+                             ( sk(k,j,i+1) + sk(k,j,i)   )                    &
+                   -      (  8.0_wp * ibit2 * adv_sca_5                       &
+                       +              ibit1 * adv_sca_3                       &
+                          ) *                                                 &
+                             ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
+                   +      (           ibit2 * adv_sca_5                       &
+                          ) *                                                 &
+                             ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
+                                     )
+
+                diss_r(k) = -ABS( u_comp ) * (                                &
+                          ( 10.0_wp * ibit2 * adv_sca_5                       &
+                       +     3.0_wp * ibit1 * adv_sca_3                       &
+                       +              ibit0 * adv_sca_1                       &
+                          ) *                                                 &
+                             ( sk(k,j,i+1) - sk(k,j,i)   )                    &
+                   -      (  5.0_wp * ibit2 * adv_sca_5                       &
+                       +              ibit1 * adv_sca_3                       &
+                          ) *                                                 &
+                             ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
+                   +      (           ibit2 * adv_sca_5                       &
+                          ) *                                                 &
+                             ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
+                                             )
+
+                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
+                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
+                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
+
+                v_comp    = v(k,j+1,i) - v_gtrans
+                flux_n(k) = v_comp * (                                        &
+                          ( 37.0_wp * ibit5 * adv_sca_5                       &
+                       +     7.0_wp * ibit4 * adv_sca_3                       &
+                       +              ibit3 * adv_sca_1                       &
+                          ) *                                                 &
+                             ( sk(k,j+1,i) + sk(k,j,i)   )                    &
+                   -      (  8.0_wp * ibit5 * adv_sca_5                       &
+                       +              ibit4 * adv_sca_3                       &
+                          ) *                                                 &
+                             ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
+                   +      (           ibit5 * adv_sca_5                       &
+                          ) *                                                 &
+                             ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
+                                     )
+
+                diss_n(k) = -ABS( v_comp ) * (                                &
+                          ( 10.0_wp * ibit5 * adv_sca_5                       &
+                       +     3.0_wp * ibit4 * adv_sca_3                       &
+                       +              ibit3 * adv_sca_1                       &
+                          ) *                                                 &
+                             ( sk(k,j+1,i) - sk(k,j,i)    )                   &
+                   -      (  5.0_wp * ibit5 * adv_sca_5                       &
+                       +              ibit4 * adv_sca_3                       &
+                          ) *                                                 &
+                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                   &
+                   +      (           ibit5 * adv_sca_5                       &
+                          ) *                                                 &
+                             ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
+                                             )
+!
+!--             k index has to be modified near bottom and top, else array
+!--             subscripts will be exceeded.
+                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
+                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
+                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
+
+                k_ppp = k + 3 * ibit8
+                k_pp  = k + 2 * ( 1 - ibit6  )
+                k_mm  = k - 2 * ibit8
+
+
+                flux_t(k) = w(k,j,i) * (                                      &
+                           ( 37.0_wp * ibit8 * adv_sca_5                      &
+                        +     7.0_wp * ibit7 * adv_sca_3                      &
+                        +           ibit6 * adv_sca_1                         &
+                           ) *                                                &
+                                   ( sk(k+1,j,i)  + sk(k,j,i)    )            &
+                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
+                        +              ibit7 * adv_sca_3                      &
+                           ) *                                                &
+                                   ( sk(k_pp,j,i) + sk(k-1,j,i)  )            &
+                    +      (           ibit8 * adv_sca_5                      &
+                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
+                                       )
+
+                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
+                           ( 10.0_wp * ibit8 * adv_sca_5                      &
+                        +     3.0_wp * ibit7 * adv_sca_3                      &
+                        +              ibit6 * adv_sca_1                      &
+                           ) *                                                &
+                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
+                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
+                        +              ibit7 * adv_sca_3                      &
+                           ) *                                                &
+                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
+                    +      (           ibit8 * adv_sca_5                      &
+                           ) *                                                &
+                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
+                                               )
+!
+!--             Apply monotonic adjustment.
+                IF ( monotonic_adjustment )  THEN
+!
+!--                At first, calculate first order fluxes.
+                   u_comp = u(k,j,i) - u_gtrans
+                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
+                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
+                             ) * adv_sca_1 
+
+                   u_comp = u(k,j,i+1) - u_gtrans
+                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
+                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
+                             ) * adv_sca_1 
+
+                   v_comp = v(k,j,i) - v_gtrans
+                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
+                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
+                             ) * adv_sca_1 
+
+                   v_comp = v(k,j+1,i) - v_gtrans
+                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
+                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
+                             ) * adv_sca_1 
+
+                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
+                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
+                            ) * adv_sca_1 
+
+                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
+                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
+                            ) * adv_sca_1 
+!
+!--                Calculate ratio of upwind gradients. Note, Min/Max is just
+!--                to avoid if statements.
+                   rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *               & 
+                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, u(k,j,i) - u_gtrans ) *               &
+                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
+                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
+
+                   rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *             & 
+                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
+                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, u(k,j,i+1) - u_gtrans ) *             &
+                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
+
+                   rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *               & 
+                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, v(k,j,i) - v_gtrans ) *               &
+                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
+                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
+
+                   rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *             & 
+                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
+                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, v(k,j+1,i) - v_gtrans ) *             &
+                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
+!   
+!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios. 
+!--                Note, for vertical advection below the third grid point above 
+!--                surface ( or below the model top) rd and rt are set to 0, i.e. 
+!--                use of first order scheme is enforced.
+                   k_mmm  = k - 3 * ibit8
+
+                   rd     = ( MAX( 0.0, w(k-1,j,i) ) *                        & 
+                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
+                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
+                               ) +                                            & 
+                              MIN( 0.0, w(k-1,j,i) ) *                        &
+                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
+                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
+                               )                                              &
+                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 
+ 
+                   rt     = ( MAX( 0.0, w(k,j,i) ) *                          & 
+                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
+                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
+                               ) +                                            & 
+                              MIN( 0.0, w(k,j,i) ) *                          &
+                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
+                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
+                               )                                              &
+                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
+!
+!--                Calculate empirical limiter function (van Albada2 limiter).
+                   phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                     &
+                                        ( rl**2 + 1.0_wp ) ) 
+                   phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                     &
+                                        ( rr**2 + 1.0_wp ) ) 
+                   phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                     &
+                                        ( rs**2 + 1.0_wp ) ) 
+                   phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                     &
+                                        ( rn**2 + 1.0_wp ) ) 
+                   phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                     &
+                                        ( rd**2 + 1.0_wp ) ) 
+                   phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                     &
+                                        ( rt**2 + 1.0_wp ) ) 
+!
+!--                Calculate the resulting monotone flux. 
+                   swap_flux_x_local(k,j) = fl_1 - phi_l *                    &
+                                          ( fl_1 - swap_flux_x_local(k,j) )
+                   flux_r(k)              = fr_1 - phi_r *                    &
+                                          ( fr_1 - flux_r(k)              )
+                   swap_flux_y_local(k)   = fs_1 - phi_s *                    &
+                                          ( fs_1 - swap_flux_y_local(k)   )
+                   flux_n(k)              = fn_1 - phi_n *                    &
+                                          ( fn_1 - flux_n(k)              )
+                   flux_d                 = fd_1 - phi_d *                    &
+                                          ( fd_1 - flux_d                 )
+                   flux_t(k)              = ft_1 - phi_t *                    &
+                                          ( ft_1 - flux_t(k)              )
+!
+!--                Moreover, modify dissipation flux according to the limiter. 
+                   swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l
+                   diss_r(k)              = diss_r(k)              * phi_r
+                   swap_diss_y_local(k)   = swap_diss_y_local(k)   * phi_s
+                   diss_n(k)              = diss_n(k)              * phi_n
+                   diss_d                 = diss_d                 * phi_d
+                   diss_t(k)              = diss_t(k)              * phi_t
+
+                ENDIF
+!
+!--             Calculate the divergence of the velocity field. A respective
+!--             correction is needed to overcome numerical instabilities caused
+!--             by a not sufficient reduction of divergences near topography.
+                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
+                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
+                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+
+                tend(k,j,i) = tend(k,j,i) - (                                 &
+                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
+                          swap_diss_x_local(k,j)            ) * ddx           &
+                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
+                          swap_diss_y_local(k)              ) * ddy           &
+                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
+                                                               ) * ddzw(k)    &
+                                            ) + sk(k,j,i) * div
+
+                swap_flux_y_local(k)   = flux_n(k)
+                swap_diss_y_local(k)   = diss_n(k)
+                swap_flux_x_local(k,j) = flux_r(k)
+                swap_diss_x_local(k,j) = diss_r(k)
+                flux_d                 = flux_t(k)
+                diss_d                 = diss_t(k)
+
+             ENDDO
+
+             DO  k = nzb_max+1, nzt
+
+                u_comp    = u(k,j,i+1) - u_gtrans
+                flux_r(k) = u_comp * (                                        &
+                      37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
+                    -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
+                    +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
+                diss_r(k) = -ABS( u_comp ) * (                                &
+                      10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
+                    -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
+                    +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
+
+                v_comp    = v(k,j+1,i) - v_gtrans
+                flux_n(k) = v_comp * (                                        &
+                      37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
+                    -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
+                    +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
+                diss_n(k) = -ABS( v_comp ) * (                                &
+                      10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
+                    -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
+                    +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
+!
+!--             k index has to be modified near bottom and top, else array
+!--             subscripts will be exceeded.
+                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
+                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
+                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
+
+                k_ppp = k + 3 * ibit8
+                k_pp  = k + 2 * ( 1 - ibit6  )
+                k_mm  = k - 2 * ibit8
+
+
+                flux_t(k) = w(k,j,i) * (                                      &
+                           ( 37.0_wp * ibit8 * adv_sca_5                      &
+                        +     7.0_wp * ibit7 * adv_sca_3                      &
+                        +              ibit6 * adv_sca_1                      &
+                           ) *                                                &
+                                   ( sk(k+1,j,i)  + sk(k,j,i)     )           &
+                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
+                        +              ibit7 * adv_sca_3                      &
+                           ) *                                                &
+                                   ( sk(k_pp,j,i) + sk(k-1,j,i)   )           &
+                    +      (           ibit8 * adv_sca_5                      &
+                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i)  )           &
+                                       )
+
+                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
+                           ( 10.0_wp * ibit8 * adv_sca_5                      &
+                        +     3.0_wp * ibit7 * adv_sca_3                      &
+                        +              ibit6 * adv_sca_1                      &
+                           ) *                                                &
+                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
+                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
+                        +              ibit7 * adv_sca_3                      &
+                           ) *                                                &
+                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
+                    +      (           ibit8 * adv_sca_5                      &
+                           ) *                                                &
+                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
+                                               )
+!
+!--             Apply monotonic adjustment.
+                IF ( monotonic_adjustment )  THEN
+!
+!--                At first, calculate first order fluxes.
+                   u_comp = u(k,j,i) - u_gtrans
+                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
+                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
+                             ) * adv_sca_1 
+
+                   u_comp = u(k,j,i+1) - u_gtrans
+                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
+                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
+                             ) * adv_sca_1 
+
+                   v_comp = v(k,j,i) - v_gtrans
+                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
+                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
+                             ) * adv_sca_1 
+
+                   v_comp = v(k,j+1,i) - v_gtrans
+                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
+                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
+                             ) * adv_sca_1 
+
+                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
+                         -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
+                            ) * adv_sca_1 
+
+                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
+                         -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
+                            ) * adv_sca_1 
+!
+!--                Calculate ratio of upwind gradients. Note, Min/Max is just
+!--                to avoid if statements.
+                   rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *               & 
+                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, u(k,j,i) - u_gtrans ) *               &
+                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
+                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
+
+                   rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *             & 
+                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
+                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, u(k,j,i+1) - u_gtrans ) *             &
+                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
+
+                   rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *               & 
+                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, v(k,j,i) - v_gtrans ) *               &
+                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
+                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
+
+                   rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *             & 
+                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
+                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, v(k,j+1,i) - v_gtrans ) *             &
+                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
+!   
+!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios. 
+!--                Note, for vertical advection below the third grid point above 
+!--                surface ( or below the model top) rd and rt are set to 0, i.e. 
+!--                use of first order scheme is enforced.
+                   k_mmm  = k - 3 * ibit8
+
+                   rd     = ( MAX( 0.0, w(k-1,j,i) ) *                        & 
+                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
+                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
+                               ) +                                            & 
+                              MIN( 0.0, w(k-1,j,i) ) *                        &
+                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
+                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
+                               )                                              &
+                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 
+ 
+                   rt     = ( MAX( 0.0, w(k,j,i) ) *                          & 
+                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
+                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
+                               ) +                                            & 
+                              MIN( 0.0, w(k,j,i) ) *                          &
+                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
+                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
+                               )                                              &
+                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 
+!
+!--                Calculate empirical limiter function (van Albada2 limiter).
+                   phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                     &
+                                        ( rl**2 + 1.0_wp ) ) 
+                   phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                     &
+                                        ( rr**2 + 1.0_wp ) ) 
+                   phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                     &
+                                        ( rs**2 + 1.0_wp ) ) 
+                   phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                     &
+                                        ( rn**2 + 1.0_wp ) ) 
+                   phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                     &
+                                        ( rd**2 + 1.0_wp ) ) 
+                   phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                     &
+                                        ( rt**2 + 1.0_wp ) ) 
+!
+!--                Calculate the resulting monotone flux. 
+                   swap_flux_x_local(k,j) = fl_1 - phi_l *                    &
+                                          ( fl_1 - swap_flux_x_local(k,j) )
+                   flux_r(k)              = fr_1 - phi_r *                    &
+                                          ( fr_1 - flux_r(k)              )
+                   swap_flux_y_local(k)   = fs_1 - phi_s *                    &
+                                          ( fs_1 - swap_flux_y_local(k)   )
+                   flux_n(k)              = fn_1 - phi_n *                    &
+                                          ( fn_1 - flux_n(k)              )
+                   flux_d                 = fd_1 - phi_d *                    &
+                                          ( fd_1 - flux_d                 )
+                   flux_t(k)              = ft_1 - phi_t *                    &
+                                          ( ft_1 - flux_t(k)              )
+!
+!--                Moreover, modify dissipation flux according to the limiter. 
+                   swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l
+                   diss_r(k)              = diss_r(k)              * phi_r
+                   swap_diss_y_local(k)   = swap_diss_y_local(k)   * phi_s
+                   diss_n(k)              = diss_n(k)              * phi_n
+                   diss_d                 = diss_d                 * phi_d
+                   diss_t(k)              = diss_t(k)              * phi_t
+
+                ENDIF
+!
+!--             Calculate the divergence of the velocity field. A respective
+!--             correction is needed to overcome numerical instabilities introduced
+!--             by a not sufficient reduction of divergences near topography.
+                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
+                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
+                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+
+                tend(k,j,i) = tend(k,j,i) - (                                 &
+                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
+                          swap_diss_x_local(k,j)            ) * ddx           &
+                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
+                          swap_diss_y_local(k)              ) * ddy           &
+                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
+                                                               ) * ddzw(k)    &
+                                            ) + sk(k,j,i) * div
+
+                swap_flux_y_local(k)   = flux_n(k)
+                swap_diss_y_local(k)   = diss_n(k)
+                swap_flux_x_local(k,j) = flux_r(k)
+                swap_diss_x_local(k,j) = diss_r(k)
+                flux_d                 = flux_t(k)
+                diss_d                 = diss_t(k)
+
+             ENDDO
+!
+!--          Evaluation of statistics. Please note, in case of using monotone
+!--          limiter vertical fluxes will be not calculated by the advective 
+!--          fluxes.
+             IF ( .NOT. monotonic_adjustment )  THEN
+
+                SELECT CASE ( sk_char )
+
+                    CASE ( 'pt' )
+                       DO  k = nzb, nzt
+                          sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)       &
+                           + ( flux_t(k) + diss_t(k) )                        &
+                           *   weight_substep(intermediate_timestep_count)
+                       ENDDO
+                    CASE ( 'sa' )
+                       DO  k = nzb, nzt
+                          sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)       &
+                           + ( flux_t(k) + diss_t(k) )                        &
+                           *   weight_substep(intermediate_timestep_count)
+                       ENDDO
+                    CASE ( 'q' )
+                       DO  k = nzb, nzt
+                          sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)         &
+                           + ( flux_t(k) + diss_t(k) )                        &
+                           *   weight_substep(intermediate_timestep_count)
+                       ENDDO
+                    CASE ( 'qr' )
+                       DO  k = nzb, nzt
+                          sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)       &
+                           + ( flux_t(k) + diss_t(k) )                        &
+                           *   weight_substep(intermediate_timestep_count)
+                       ENDDO
+                    CASE ( 'nr' )
+                       DO  k = nzb, nzt
+                          sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)       &
+                           + ( flux_t(k) + diss_t(k) )                        &
+                           *   weight_substep(intermediate_timestep_count)
+                       ENDDO
+
+                 END SELECT
+
+              ENDIF
+
+         ENDDO
+      ENDDO
+
+    END SUBROUTINE advec_s_ws
+
+
+!------------------------------------------------------------------------------!
+! Scalar advection - Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
+
+       USE arrays_3d,                                                         &
+           ONLY:  ddzw, tend, u, v, w
+
+       USE constants,                                                         &
+           ONLY:  adv_sca_1, adv_sca_3, adv_sca_5
+
+       USE control_parameters,                                                &
+           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
+                  v_gtrans
+
+       USE grid_variables,                                                    &
+           ONLY:  ddx, ddy
+
+       USE indices,                                                           &
+           ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,  &
+                  nzb, nzb_max, nzt, wall_flags_0
+
+       USE kinds
+       
+!        USE statistics,                                                       &
+!            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
+!                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN = *), INTENT(IN)    :: sk_char !:
+
+       INTEGER(iwp) ::  i      !:
+       INTEGER(iwp) ::  ibit0  !:
+       INTEGER(iwp) ::  ibit1  !:
+       INTEGER(iwp) ::  ibit2  !:
+       INTEGER(iwp) ::  ibit3  !:
+       INTEGER(iwp) ::  ibit4  !:
+       INTEGER(iwp) ::  ibit5  !:
+       INTEGER(iwp) ::  ibit6  !:
+       INTEGER(iwp) ::  ibit7  !:
+       INTEGER(iwp) ::  ibit8  !:
+       INTEGER(iwp) ::  j      !:
+       INTEGER(iwp) ::  k      !:
+       INTEGER(iwp) ::  k_mm   !:
+       INTEGER(iwp) ::  k_mmm  !:
+       INTEGER(iwp) ::  k_pp   !:
+       INTEGER(iwp) ::  k_ppp  !:
+       INTEGER(iwp) ::  tn = 0 !:
 
        REAL(wp)    ::  diss_d !:
@@ -2747,4 +3352,22 @@
        REAL(wp)    ::  flux_s !:
        REAL(wp)    ::  flux_t !:
+       REAL(wp)    ::  fd_1   !:
+       REAL(wp)    ::  fl_1   !:
+       REAL(wp)    ::  fn_1   !:
+       REAL(wp)    ::  fr_1   !:
+       REAL(wp)    ::  fs_1   !:
+       REAL(wp)    ::  ft_1   !:
+       REAL(wp)    ::  phi_d  !:
+       REAL(wp)    ::  phi_l  !:
+       REAL(wp)    ::  phi_n  !:
+       REAL(wp)    ::  phi_r  !:
+       REAL(wp)    ::  phi_s  !:
+       REAL(wp)    ::  phi_t  !:
+       REAL(wp)    ::  rd     !:
+       REAL(wp)    ::  rl     !:
+       REAL(wp)    ::  rn     !:
+       REAL(wp)    ::  rr     !:
+       REAL(wp)    ::  rs     !:
+       REAL(wp)    ::  rt     !:
        REAL(wp)    ::  u_comp !:
        REAL(wp)    ::  v_comp !:
@@ -2967,4 +3590,136 @@
                                    ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
                                          )
+!
+!--             Apply monotonic adjustment.
+                IF ( monotonic_adjustment )  THEN
+!
+!--                At first, calculate first order fluxes.
+                   u_comp = u(k,j,i) - u_gtrans
+                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
+                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
+                             ) * adv_sca_1 
+
+                   u_comp = u(k,j,i+1) - u_gtrans
+                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
+                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
+                             ) * adv_sca_1 
+
+                   v_comp = v(k,j,i) - v_gtrans
+                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
+                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
+                             ) * adv_sca_1 
+
+                   v_comp = v(k,j+1,i) - v_gtrans
+                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
+                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
+                             ) * adv_sca_1 
+
+                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
+                        -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )      &
+                            ) * adv_sca_1 
+
+                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
+                        -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )        &
+                            ) * adv_sca_1 
+!
+!--                Calculate ratio of upwind gradients. Note, Min/Max is just
+!--                to avoid if statements.
+                   rl     = ( MAX( 0.0, u(k,j,i) - u_gtrans ) *               & 
+                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, u(k,j,i) - u_gtrans ) *               &
+                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
+                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
+
+                   rr     = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) *             & 
+                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
+                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, u(k,j,i+1) - u_gtrans ) *             &
+                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
+
+                   rs     = ( MAX( 0.0, v(k,j,i) - v_gtrans ) *               & 
+                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, v(k,j,i) - v_gtrans ) *               &
+                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
+                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
+
+                   rn     = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) *             & 
+                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
+                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
+                                  ) +                                         & 
+                              MIN( 0.0, v(k,j+1,i) - v_gtrans ) *             &
+                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
+                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
+                                  )                                           &
+                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
+!   
+!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios. 
+!--                Note, for vertical advection below the third grid point above 
+!--                surface ( or below the model top) rd and rt are set to 0, i.e. 
+!--                use of first order scheme is enforced.
+                   k_mmm  = k - 3 * ibit8
+
+                   rd     = ( MAX( 0.0, w(k-1,j,i) ) *                        & 
+                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
+                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
+                               ) +                                            & 
+                              MIN( 0.0, w(k-1,j,i) ) *                        &
+                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
+                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
+                               )                                              &
+                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 
+ 
+                   rt     = ( MAX( 0.0, w(k,j,i) ) *                          & 
+                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
+                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
+                               ) +                                            & 
+                              MIN( 0.0, w(k,j,i) ) *                          &
+                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
+                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
+                               )                                              &
+                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
+!
+!--                Calculate empirical limiter function (van Albada2 limiter).
+                   phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) /                     &
+                                        ( rl**2 + 1.0_wp ) ) 
+                   phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) /                     &
+                                        ( rr**2 + 1.0_wp ) ) 
+                   phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) /                     &
+                                        ( rs**2 + 1.0_wp ) ) 
+                   phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) /                     &
+                                        ( rn**2 + 1.0_wp ) ) 
+                   phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) /                     &
+                                        ( rd**2 + 1.0_wp ) ) 
+                   phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) /                     &
+                                        ( rt**2 + 1.0_wp ) ) 
+!
+!--                Calculate the resulting monotone flux. 
+                   flux_l = fl_1 - phi_l * ( fl_1 - flux_l )
+                   flux_r = fr_1 - phi_r * ( fr_1 - flux_r )
+                   flux_s = fs_1 - phi_s * ( fs_1 - flux_s )
+                   flux_n = fn_1 - phi_n * ( fn_1 - flux_n )
+                   flux_d = fd_1 - phi_d * ( fd_1 - flux_d )
+                   flux_t = ft_1 - phi_t * ( ft_1 - flux_t )
+!
+!--                Moreover, modify dissipation flux according to the limiter. 
+                   diss_l = diss_l * phi_l
+                   diss_r = diss_r * phi_r
+                   diss_s = diss_s * phi_s
+                   diss_n = diss_n * phi_n
+                   diss_d = diss_d * phi_d
+                   diss_t = diss_t * phi_t
+
+                ENDIF
 !
 !--             Calculate the divergence of the velocity field. A respective
