Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 4325)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 4327)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Setting of advection flags for vertical fluxes of w revised, air density for
+! vertical flux calculation of w at k=1 is considered now
+! 
+! 4325 2019-12-06 07:14:04Z Giersch
 ! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of 
 ! advection flags for fluxes in z-direction revised, comments extended
@@ -748,12 +752,9 @@
                 
                 flag_set = .FALSE.
-                IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)   .AND.            &
-                       .NOT. BTEST(wall_flags_0(k,j,i),3)     .AND.            &
+                IF ( ( .NOT. BTEST(wall_flags_0(k,j,i),3)     .AND.            &
                              BTEST(wall_flags_0(k+1,j,i),3) ) .OR.             &
-                     ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)   .AND.            &
-                             BTEST(wall_flags_0(k,j,i),3) )   .OR.             &
                      ( .NOT. BTEST(wall_flags_0(k+1,j,i),3)   .AND.            &
                              BTEST(wall_flags_0(k,j,i),3) )   .OR.             &        
-                     k == nzt )                                                &
+                     k == nzt -1 )                                             &
                 THEN
 !
@@ -765,20 +766,21 @@
                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
                    flag_set = .TRUE.
-                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3)    .OR.       &
-                           .NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND.      &
-                                 BTEST(wall_flags_0(k-1,j,i),3)     .AND.      &
-                                 BTEST(wall_flags_0(k,j,i),3)       .AND.      &
-                                 BTEST(wall_flags_0(k+1,j,i),3)     .AND.      &
-                           .NOT. flag_set                           .OR.       & 
-                         k == nzt - 1 )                                        &
+                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)    .AND.       &
+                                 BTEST(wall_flags_0(k,j,i),3)      .AND.       &
+                                 BTEST(wall_flags_0(k+1,j,i),3)    .AND.       &
+                                 BTEST(wall_flags_0(k_pp,j,i),3) ) .OR.        &
+                         ( .NOT. BTEST(wall_flags_0(k_pp,j,i),3)   .AND.       &
+                                 BTEST(wall_flags_0(k+1,j,i),3)    .AND.       &
+                                 BTEST(wall_flags_0(k,j,i),3)      .AND.       &
+                                 BTEST(wall_flags_0(k-1,j,i),3) )  .AND.       &
+                           .NOT. flag_set                          .OR.        & 
+                         k == nzt - 2 )                                        &
                 THEN                                                           
                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )      
                    flag_set = .TRUE.                                           
-                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3)  .AND.                &
-                         BTEST(wall_flags_0(k-1,j,i),3)   .AND.                &
+                ELSEIF ( BTEST(wall_flags_0(k-1,j,i),3)   .AND.                &
                          BTEST(wall_flags_0(k,j,i),3)     .AND.                &
                          BTEST(wall_flags_0(k+1,j,i),3)   .AND.                &
                          BTEST(wall_flags_0(k_pp,j,i),3)  .AND.                &
-                         BTEST(wall_flags_0(k_ppp,j,i),3) .AND.                &
                          .NOT. flag_set )                                      &
                 THEN
@@ -3355,6 +3357,8 @@
        k         = nzb + 1
        w_comp(k) = w(k,j,i) + w(k-1,j,i)
-       flux_t(0) = w_comp(k)       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
-       diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
+       flux_t(0) = w_comp(k)       * rho_air(k)                                &
+                 * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
+       diss_t(0) = -ABS(w_comp(k)) * rho_air(k)                                &
+                 * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
        
        DO  k = nzb+1, nzb+1
@@ -4574,5 +4578,5 @@
                                        ( 37.0_wp * ibit2 * adv_mom_5           &
                                     +     7.0_wp * ibit1 * adv_mom_3           &
-                                    +              ibit0 * adv_mom_1          &
+                                    +              ibit0 * adv_mom_1           &
                                        ) *                                     &
                                      ( u(k,j,i)   + u(k,j,i-1) )               &
@@ -5042,32 +5046,32 @@
 
                 w_comp = w(k,j,i) + w(k,j,i-1)
-                flux_t = w_comp * rho_air_zw(k) * (                          &
+                flux_t = w_comp * rho_air_zw(k) * (                            &
                           ( 37.0_wp * ibit8 * adv_mom_5                        &
                        +     7.0_wp * ibit7 * adv_mom_3                        &
                        +              ibit6 * adv_mom_1                        &
-                          ) *                                                &
-                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
+                          ) *                                                  &
+                             ( u(k+1,j,i)  + u(k,j,i)     )                    &
                    -      (  8.0_wp * ibit8 * adv_mom_5                        &
                        +              ibit7 * adv_mom_3                        &
-                          ) *                                                &
-                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
+                          ) *                                                  &
+                             ( u(k_pp,j,i) + u(k-1,j,i)   )                    &
                    +      (           ibit8 * adv_mom_5                        &
-                          ) *                                                &
-                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
+                          ) *                                                  &
+                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                    &
                                       )
 
-                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                   &
                           ( 10.0_wp * ibit8 * adv_mom_5                        &
                        +     3.0_wp * ibit7 * adv_mom_3                        &
                        +              ibit6 * adv_mom_1                        &
-                          ) *                                                &
-                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
+                          ) *                                                  &
+                             ( u(k+1,j,i)   - u(k,j,i)    )                    &
                    -      (  5.0_wp * ibit8 * adv_mom_5                        &
                        +              ibit7 * adv_mom_3                        &
-                          ) *                                                &
-                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
+                          ) *                                                  &
+                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                    &
                    +      (           ibit8 * adv_mom_5                        &
-                           ) *                                               &
-                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
+                           ) *                                                 &
+                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                    &
                                               )
 !
@@ -5075,9 +5079,9 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
-                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
-                     +  (   w_comp                      * rho_air_zw(k) -     &
-                          ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)     &
-                        ) * drho_air(k) * ddzw(k)                             &
+                div = ( ( u_comp      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
+                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
+                     +  (   w_comp                      * rho_air_zw(k) -      &
+                          ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)      &
+                        ) * drho_air(k) * ddzw(k)                              &
                       ) * 0.5_wp
 
@@ -5231,14 +5235,14 @@
              u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
              swap_flux_x_local_v(k,j) = u_comp * (                             &
-                                      ( 37.0_wp * ibit11 * adv_mom_5              &
-                                   +     7.0_wp * ibit10 * adv_mom_3              &
-                                   +              ibit9 * adv_mom_1              &
+                                      ( 37.0_wp * ibit11 * adv_mom_5           &
+                                   +     7.0_wp * ibit10 * adv_mom_3           &
+                                   +              ibit9 * adv_mom_1            &
                                       ) *                                      &
                                      ( v(k,j,i)   + v(k,j,i-1) )               &
-                               -      (  8.0_wp * ibit11 * adv_mom_5              &
-                                   +              ibit10 * adv_mom_3              &
+                               -      (  8.0_wp * ibit11 * adv_mom_5           &
+                                   +              ibit10 * adv_mom_3           &
                                       ) *                                      &
                                      ( v(k,j,i+1) + v(k,j,i-2) )               &
-                               +      (           ibit11 * adv_mom_5              &
+                               +      (           ibit11 * adv_mom_5           &
                                       ) *                                      &
                                      ( v(k,j,i+2) + v(k,j,i-3) )               &
@@ -5246,14 +5250,14 @@
 
               swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                   &
-                                      ( 10.0_wp * ibit11 * adv_mom_5              &
-                                   +     3.0_wp * ibit10 * adv_mom_3              &
-                                   +              ibit9 * adv_mom_1              &
+                                      ( 10.0_wp * ibit11 * adv_mom_5           &
+                                   +     3.0_wp * ibit10 * adv_mom_3           &
+                                   +              ibit9 * adv_mom_1            &
                                       ) *                                      &
                                      ( v(k,j,i)   - v(k,j,i-1) )               &
-                               -      (  5.0_wp * ibit11 * adv_mom_5              &
-                                   +              ibit10 * adv_mom_3              &
+                               -      (  5.0_wp * ibit11 * adv_mom_5           &
+                                   +              ibit10 * adv_mom_3           &
                                       ) *                                      &
                                      ( v(k,j,i+1) - v(k,j,i-2) )               &
-                               +      (           ibit11 * adv_mom_5              &
+                               +      (           ibit11 * adv_mom_5           &
                                       ) *                                      &
                                      ( v(k,j,i+2) - v(k,j,i-3) )               &
@@ -5265,11 +5269,11 @@
 
              u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
-             swap_flux_x_local_v(k,j) = u_comp * (                            &
-                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )               &
-                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )               &
+             swap_flux_x_local_v(k,j) = u_comp * (                             &
+                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )             &
+                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )             &
                            +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
-             swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                   &
-                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )               &
-                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )               &
+             swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                    &
+                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )             &
+                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )             &
                            +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
 
@@ -5307,32 +5311,32 @@
 
              v_comp                 = v(k,j,i) + v(k,j-1,i) - gv
-             swap_flux_y_local_v(k) = v_comp * (                              &
-                                   ( 37.0_wp * ibit14 * adv_mom_5                &
-                                +     7.0_wp * ibit13 * adv_mom_3                &
-                                +              ibit12 * adv_mom_1                &
-                                   ) *                                        &
-                                     ( v(k,j,i)   + v(k,j-1,i) )              &
-                            -      (  8.0_wp * ibit14 * adv_mom_5                &
-                                +              ibit13 * adv_mom_3                &
-                                   ) *                                        &
-                                     ( v(k,j+1,i) + v(k,j-2,i) )              &
-                            +      (           ibit14 * adv_mom_5                &
-                                   ) *                                        &
-                                     ( v(k,j+2,i) + v(k,j-3,i) )              &
+             swap_flux_y_local_v(k) = v_comp * (                               &
+                                   ( 37.0_wp * ibit14 * adv_mom_5              &
+                                +     7.0_wp * ibit13 * adv_mom_3              &
+                                +              ibit12 * adv_mom_1              &
+                                   ) *                                         &
+                                     ( v(k,j,i)   + v(k,j-1,i) )               &
+                            -      (  8.0_wp * ibit14 * adv_mom_5              &
+                                +              ibit13 * adv_mom_3              &
+                                   ) *                                         &
+                                     ( v(k,j+1,i) + v(k,j-2,i) )               &
+                            +      (           ibit14 * adv_mom_5              &
+                                   ) *                                         &
+                                     ( v(k,j+2,i) + v(k,j-3,i) )               &
                                                  )
 
-             swap_diss_y_local_v(k) = - ABS( v_comp ) * (                     &
-                                   ( 10.0_wp * ibit14 * adv_mom_5                &
-                                +     3.0_wp * ibit13 * adv_mom_3                &
-                                +              ibit12 * adv_mom_1                &
-                                   ) *                                        &
-                                     ( v(k,j,i)   - v(k,j-1,i) )              &
-                            -      (  5.0_wp * ibit14 * adv_mom_5                &
-                                +              ibit13 * adv_mom_3                &
-                                   ) *                                        &
-                                     ( v(k,j+1,i) - v(k,j-2,i) )              &
-                            +      (           ibit14 * adv_mom_5                &
-                                   ) *                                        &
-                                     ( v(k,j+2,i) - v(k,j-3,i) )              &
+             swap_diss_y_local_v(k) = - ABS( v_comp ) * (                      &
+                                   ( 10.0_wp * ibit14 * adv_mom_5              &
+                                +     3.0_wp * ibit13 * adv_mom_3              &
+                                +              ibit12 * adv_mom_1              &
+                                   ) *                                         &
+                                     ( v(k,j,i)   - v(k,j-1,i) )               &
+                            -      (  5.0_wp * ibit14 * adv_mom_5              &
+                                +              ibit13 * adv_mom_3              &
+                                   ) *                                         &
+                                     ( v(k,j+1,i) - v(k,j-2,i) )               &
+                            +      (           ibit14 * adv_mom_5              &
+                                   ) *                                         &
+                                     ( v(k,j+2,i) - v(k,j-3,i) )               &
                                                           )
 
@@ -5342,11 +5346,11 @@
 
              v_comp                 = v(k,j,i) + v(k,j-1,i) - gv
-             swap_flux_y_local_v(k) = v_comp * (                              &
-                           37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )                 &
-                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )                 &
+             swap_flux_y_local_v(k) = v_comp * (                               &
+                           37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )               &
+                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )               &
                          +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
-             swap_diss_y_local_v(k) = - ABS( v_comp ) * (                     &
-                           10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )                 &
-                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )                 &
+             swap_diss_y_local_v(k) = - ABS( v_comp ) * (                      &
+                           10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )               &
+                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )               &
                          +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
 
@@ -5366,32 +5370,32 @@
 
                 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
-                flux_r = u_comp * (                                          &
-                          ( 37.0_wp * ibit11 * adv_mom_5                        &
-                       +     7.0_wp * ibit10 * adv_mom_3                        &
+                flux_r = u_comp * (                                            &
+                          ( 37.0_wp * ibit11 * adv_mom_5                       &
+                       +     7.0_wp * ibit10 * adv_mom_3                       &
                        +              ibit9 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit11 * adv_mom_5                        &
-                       +              ibit10 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j,i+2) + v(k,j,i-1) )                 &
-                   +      (           ibit11 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j,i+3) + v(k,j,i-2) )                 &
+                          ) *                                                  &
+                                 ( v(k,j,i+1) + v(k,j,i)   )                   &
+                   -      (  8.0_wp * ibit11 * adv_mom_5                       &
+                       +              ibit10 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( v(k,j,i+2) + v(k,j,i-1) )                   &
+                   +      (           ibit11 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( v(k,j,i+3) + v(k,j,i-2) )                   &
                                      )
 
-                diss_r = - ABS( u_comp ) * (                                 &
-                          ( 10.0_wp * ibit11 * adv_mom_5                        &
-                       +     3.0_wp * ibit10 * adv_mom_3                        &
+                diss_r = - ABS( u_comp ) * (                                   &
+                          ( 10.0_wp * ibit11 * adv_mom_5                       &
+                       +     3.0_wp * ibit10 * adv_mom_3                       &
                        +              ibit9 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit11 * adv_mom_5                        &
-                       +              ibit10 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j,i+2) - v(k,j,i-1) )                 &
-                   +      (           ibit11 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j,i+3) - v(k,j,i-2) )                 &
+                          ) *                                                  &
+                                 ( v(k,j,i+1) - v(k,j,i)  )                    &
+                   -      (  5.0_wp * ibit11 * adv_mom_5                       &
+                       +              ibit10 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( v(k,j,i+2) - v(k,j,i-1) )                   &
+                   +      (           ibit11 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( v(k,j,i+3) - v(k,j,i-2) )                   &
                                               )
 
@@ -5404,32 +5408,32 @@
 
                 u_comp_l  = u(k,j-1,i) + u(k,j,i) - gu
-                flux_l    = u_comp_l * (                                     &
-                                      ( 37.0_wp * ibit11_l * adv_mom_5          &
-                                   +     7.0_wp * ibit10_l * adv_mom_3          &
+                flux_l    = u_comp_l * (                                       &
+                                      ( 37.0_wp * ibit11_l * adv_mom_5         &
+                                   +     7.0_wp * ibit10_l * adv_mom_3         &
                                    +              ibit9_l * adv_mom_1          &
-                                      ) *                                    &
-                                     ( v(k,j,i)   + v(k,j,i-1) )             &
-                               -      (  8.0_wp * ibit11_l * adv_mom_5          &
-                                   +              ibit10_l * adv_mom_3          &
-                                      ) *                                    &
-                                     ( v(k,j,i+1) + v(k,j,i-2) )             &
-                               +      (           ibit11_l * adv_mom_5          &
-                                      ) *                                    &
-                                     ( v(k,j,i+2) + v(k,j,i-3) )             &
+                                      ) *                                      &
+                                     ( v(k,j,i)   + v(k,j,i-1) )               &
+                               -      (  8.0_wp * ibit11_l * adv_mom_5         &
+                                   +              ibit10_l * adv_mom_3         &
+                                      ) *                                      &
+                                     ( v(k,j,i+1) + v(k,j,i-2) )               &
+                               +      (           ibit11_l * adv_mom_5         &
+                                      ) *                                      &
+                                     ( v(k,j,i+2) + v(k,j,i-3) )               &
                                                  )
 
-                 diss_l   = - ABS( u_comp_l ) * (                            &
-                                      ( 10.0_wp * ibit11_l * adv_mom_5          &
-                                   +     3.0_wp * ibit10_l * adv_mom_3          &
+                 diss_l   = - ABS( u_comp_l ) * (                              &
+                                      ( 10.0_wp * ibit11_l * adv_mom_5         &
+                                   +     3.0_wp * ibit10_l * adv_mom_3         &
                                    +              ibit9_l * adv_mom_1          &
-                                      ) *                                    &
-                                     ( v(k,j,i)   - v(k,j,i-1) )             &
-                               -      (  5.0_wp * ibit11_l * adv_mom_5          &
-                                   +              ibit10_l * adv_mom_3          &
-                                      ) *                                    &
-                                     ( v(k,j,i+1) - v(k,j,i-2) )             &
-                               +      (           ibit11_l * adv_mom_5          &
-                                      ) *                                    &
-                                     ( v(k,j,i+2) - v(k,j,i-3) )             &
+                                      ) *                                      &
+                                     ( v(k,j,i)   - v(k,j,i-1) )               &
+                               -      (  5.0_wp * ibit11_l * adv_mom_5         &
+                                   +              ibit10_l * adv_mom_3         &
+                                      ) *                                      &
+                                     ( v(k,j,i+1) - v(k,j,i-2) )               &
+                               +      (           ibit11_l * adv_mom_5         &
+                                      ) *                                      &
+                                     ( v(k,j,i+2) - v(k,j,i-3) )               &
                                                            )
 #else
@@ -5443,32 +5447,32 @@
 
                 v_comp = v(k,j+1,i) + v(k,j,i)
-                flux_n = ( v_comp - gv ) * (                                 &
-                          ( 37.0_wp * ibit14 * adv_mom_5                        &
-                       +     7.0_wp * ibit13 * adv_mom_3                        &
-                       +              ibit12 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit14 * adv_mom_5                        &
-                       +              ibit13 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j+2,i) + v(k,j-1,i) )                 &
-                   +      (           ibit14 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j+3,i) + v(k,j-2,i) )                 &
+                flux_n = ( v_comp - gv ) * (                                   &
+                          ( 37.0_wp * ibit14 * adv_mom_5                       &
+                       +     7.0_wp * ibit13 * adv_mom_3                       &
+                       +              ibit12 * adv_mom_1                       &
+                          ) *                                                  &
+                                 ( v(k,j+1,i) + v(k,j,i)   )                   &
+                   -      (  8.0_wp * ibit14 * adv_mom_5                       &
+                       +              ibit13 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( v(k,j+2,i) + v(k,j-1,i) )                   &
+                   +      (           ibit14 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( v(k,j+3,i) + v(k,j-2,i) )                   &
                                      )
 
-                diss_n = - ABS( v_comp - gv ) * (                            &
-                          ( 10.0_wp * ibit14 * adv_mom_5                        &
-                       +     3.0_wp * ibit13 * adv_mom_3                        &
-                       +              ibit12 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit14 * adv_mom_5                        &
-                       +              ibit13 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j+2,i) - v(k,j-1,i) )                 &
-                   +      (           ibit14 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j+3,i) - v(k,j-2,i) )                 &
+                diss_n = - ABS( v_comp - gv ) * (                              &
+                          ( 10.0_wp * ibit14 * adv_mom_5                       &
+                       +     3.0_wp * ibit13 * adv_mom_3                       &
+                       +              ibit12 * adv_mom_1                       &
+                          ) *                                                  &
+                                 ( v(k,j+1,i) - v(k,j,i)  )                    &
+                   -      (  5.0_wp * ibit14 * adv_mom_5                       &
+                       +              ibit13 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( v(k,j+2,i) - v(k,j-1,i) )                   &
+                   +      (           ibit14 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( v(k,j+3,i) - v(k,j-2,i) )                   &
                                                       )
 
@@ -5481,32 +5485,32 @@
 
                 v_comp_s = v(k,j,i) + v(k,j-1,i) - gv
-                flux_s   = v_comp_s * (                                      &
-                                   ( 37.0_wp * ibit14_s * adv_mom_5             &
-                                +     7.0_wp * ibit13_s * adv_mom_3             &
-                                +              ibit12_s * adv_mom_1             &
-                                   ) *                                       &
-                                     ( v(k,j,i)   + v(k,j-1,i) )             &
-                            -      (  8.0_wp * ibit14_s * adv_mom_5             &
-                                +              ibit13_s * adv_mom_3             &
-                                   ) *                                       &
-                                     ( v(k,j+1,i) + v(k,j-2,i) )             &
-                            +      (           ibit14_s * adv_mom_5             &
-                                   ) *                                       &
-                                     ( v(k,j+2,i) + v(k,j-3,i) )             &
+                flux_s   = v_comp_s * (                                        &
+                                   ( 37.0_wp * ibit14_s * adv_mom_5            &
+                                +     7.0_wp * ibit13_s * adv_mom_3            &
+                                +              ibit12_s * adv_mom_1            &
+                                   ) *                                         &
+                                     ( v(k,j,i)   + v(k,j-1,i) )               &
+                            -      (  8.0_wp * ibit14_s * adv_mom_5            &
+                                +              ibit13_s * adv_mom_3            &
+                                   ) *                                         &
+                                     ( v(k,j+1,i) + v(k,j-2,i) )               &
+                            +      (           ibit14_s * adv_mom_5            &
+                                   ) *                                         &
+                                     ( v(k,j+2,i) + v(k,j-3,i) )               &
                                                  )
 
-                diss_s   = - ABS( v_comp_s ) * (                             &
-                                   ( 10.0_wp * ibit14_s * adv_mom_5             &
-                                +     3.0_wp * ibit13_s * adv_mom_3             &
-                                +              ibit12_s * adv_mom_1             &
-                                   ) *                                       &
-                                     ( v(k,j,i)   - v(k,j-1,i) )             &
-                            -      (  5.0_wp * ibit14_s * adv_mom_5             &
-                                +              ibit13_s * adv_mom_3             &
-                                   ) *                                       &
-                                     ( v(k,j+1,i) - v(k,j-2,i) )             &
-                            +      (           ibit14_s * adv_mom_5             &
-                                   ) *                                       &
-                                     ( v(k,j+2,i) - v(k,j-3,i) )             &
+                diss_s   = - ABS( v_comp_s ) * (                               &
+                                   ( 10.0_wp * ibit14_s * adv_mom_5            &
+                                +     3.0_wp * ibit13_s * adv_mom_3            &
+                                +              ibit12_s * adv_mom_1            &
+                                   ) *                                         &
+                                     ( v(k,j,i)   - v(k,j-1,i) )               &
+                            -      (  5.0_wp * ibit14_s * adv_mom_5            &
+                                +              ibit13_s * adv_mom_3            &
+                                   ) *                                         &
+                                     ( v(k,j+1,i) - v(k,j-2,i) )               &
+                            +      (           ibit14_s * adv_mom_5            &
+                                   ) *                                         &
+                                     ( v(k,j+2,i) - v(k,j-3,i) )               &
                                                           )
 #else
@@ -5527,32 +5531,32 @@
 
                 w_comp = w(k,j-1,i) + w(k,j,i)
-                flux_t = w_comp * rho_air_zw(k) * (                          &
-                          ( 37.0_wp * ibit17 * adv_mom_5                        &
-                       +     7.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
-                          ) *                                                &
-                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
-                   -      (  8.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
-                          ) *                                                &
-                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
-                   +      (           ibit17 * adv_mom_5                        &
-                          ) *                                                &
-                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
+                flux_t = w_comp * rho_air_zw(k) * (                            &
+                          ( 37.0_wp * ibit17 * adv_mom_5                       &
+                       +     7.0_wp * ibit16 * adv_mom_3                       &
+                       +              ibit15 * adv_mom_1                       &
+                          ) *                                                  &
+                             ( v(k+1,j,i)   + v(k,j,i)    )                    &
+                   -      (  8.0_wp * ibit17 * adv_mom_5                       &
+                       +              ibit16 * adv_mom_3                       &
+                          ) *                                                  &
+                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                    &
+                   +      (           ibit17 * adv_mom_5                       &
+                          ) *                                                  &
+                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                    &
                                       )
 
-                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
-                          ( 10.0_wp * ibit17 * adv_mom_5                        &
-                       +     3.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
-                          ) *                                                &
-                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
-                   -      (  5.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
-                          ) *                                                &
-                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
-                   +      (           ibit17 * adv_mom_5                        &
-                          ) *                                                &
-                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                   &
+                          ( 10.0_wp * ibit17 * adv_mom_5                       &
+                       +     3.0_wp * ibit16 * adv_mom_3                       &
+                       +              ibit15 * adv_mom_1                       &
+                          ) *                                                  &
+                             ( v(k+1,j,i)   - v(k,j,i)    )                    &
+                   -      (  5.0_wp * ibit17 * adv_mom_5                       &
+                       +              ibit16 * adv_mom_3                       &
+                          ) *                                                  &
+                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                    &
+                   +      (           ibit17 * adv_mom_5                       &
+                          ) *                                                  &
+                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                    &
                                                )
 !
@@ -5560,41 +5564,41 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( ( u_comp     + gu )                                 &
-                                       * ( ibit9 + ibit10 + ibit11 )         &
-                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
-                                       * (                                    &
-                        REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp )   &
-                      + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )  &
-                      + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )  &
-                                         )                                    &
-                  ) * ddx                                                     &
-               +  ( v_comp                                                    &
-                                       * ( ibit12 + ibit13 + ibit14 )         &
-                - ( v(k,j,i)     + v(k,j-1,i) )                               &
-                                       * (                                    &
-                        REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )  &
-                      + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )  &
-                      + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )  &
-                                         )                                    &
-                  ) * ddy                                                     &
-               +  ( w_comp * rho_air_zw(k)                                    &
-                                       * ( ibit15 + ibit16 + ibit17 )         &
-                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
-                                       * (                                    &
-                        REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp )  &
-                      + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp )  &
-                      + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp )  &
-                                         )                                    &
-                   ) * drho_air(k) * ddzw(k)                                  &
+                div = ( ( ( u_comp     + gu )                                  &
+                                       * ( ibit9 + ibit10 + ibit11 )           &
+                - ( u(k,j-1,i)   + u(k,j,i) )                                  &
+                                       * (                                     &
+                        REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )   &
+                                         )                                     &
+                  ) * ddx                                                      &
+               +  ( v_comp                                                     &
+                                       * ( ibit12 + ibit13 + ibit14 )          &
+                - ( v(k,j,i)     + v(k,j-1,i) )                                &
+                                       * (                                     &
+                        REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )   & 
+                      + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )   &
+                                         )                                     &
+                  ) * ddy                                                      & 
+               +  ( w_comp * rho_air_zw(k)                                     &
+                                       * ( ibit15 + ibit16 + ibit17 )          &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)              &
+                                       * (                                     &
+                        REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp )   &
+                                         )                                     &
+                   ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
 
 
-                tend(k,j,i) = tend(k,j,i) - (                                 &
-                       ( ( flux_r + diss_r )                                  &
-                     -   ( flux_l + diss_l ) ) * ddx                          &
-                     + ( ( flux_n + diss_n )                                  &
-                     -   ( flux_s + diss_s ) ) * ddy                          &
-                     + ( ( flux_t + diss_t )                                  &
-                     -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)        &
+                tend(k,j,i) = tend(k,j,i) - (                                  &
+                       ( ( flux_r + diss_r )                                   &
+                     -   ( flux_l + diss_l ) ) * ddx                           &
+                     + ( ( flux_n + diss_n )                                   &
+                     -   ( flux_s + diss_s ) ) * ddy                           &
+                     + ( ( flux_t + diss_t )                                   &
+                     -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)         &
                                             )  + v(k,j,i) * div
 
@@ -5637,12 +5641,12 @@
 
                 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
-                flux_r = u_comp * (                                           &
-                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                      &
-                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                      &
+                flux_r = u_comp * (                                            &
+                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                    &
+                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                    &
                     +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
 
-                diss_r = - ABS( u_comp ) * (                                  &
-                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                        &
-                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                      &
+                diss_r = - ABS( u_comp ) * (                                   &
+                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                      &
+                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                    &
                     +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
 
@@ -5651,11 +5655,11 @@
 !--             Recompute the left fluxes.
                 u_comp_l = u(k,j-1,i) + u(k,j,i) - gu
-                flux_l   = u_comp_l * (                                       &
-                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )               &
-                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )               &
+                flux_l   = u_comp_l * (                                        &
+                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )             &
+                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )             &
                            +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
-                diss_l   = - ABS( u_comp_l ) * (                              &
-                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )               &
-                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )               &
+                diss_l   = - ABS( u_comp_l ) * (                               &
+                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )             &
+                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )             &
                            +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
 #else
@@ -5665,12 +5669,12 @@
 
                 v_comp = v(k,j+1,i) + v(k,j,i)
-                flux_n = ( v_comp - gv ) * (                                  &
-                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                      &
-                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                      &
+                flux_n = ( v_comp - gv ) * (                                   &
+                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                    &
+                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                    &
                       +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
 
-                diss_n = - ABS( v_comp - gv ) * (                             &
-                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                      &
-                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                      &
+                diss_n = - ABS( v_comp - gv ) * (                              &
+                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                    &
+                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                    &
                     +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
 
@@ -5679,11 +5683,11 @@
 !--             Recompute the south fluxes.
                 v_comp_s = v(k,j,i) + v(k,j-1,i) - gv
-                flux_s   = v_comp_s * (                                       &
-                           37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )                 &
-                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )                 &
+                flux_s   = v_comp_s * (                                        &
+                           37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )               &
+                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )               &
                          +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
-                diss_s   = - ABS( v_comp_s ) * (                              &
-                           10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )                 &
-                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )                 &
+                diss_s   = - ABS( v_comp_s ) * (                               &
+                           10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )               &
+                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )               &
                          +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
 #else
@@ -5704,32 +5708,32 @@
 
                 w_comp = w(k,j-1,i) + w(k,j,i)
-                flux_t = w_comp * rho_air_zw(k) * (                          &
-                          ( 37.0_wp * ibit17 * adv_mom_5                        &
-                       +     7.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
-                          ) *                                                &
-                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
-                   -      (  8.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
-                          ) *                                                &
-                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
-                   +      (           ibit17 * adv_mom_5                        &
-                          ) *                                                &
-                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
+                flux_t = w_comp * rho_air_zw(k) * (                            &
+                          ( 37.0_wp * ibit17 * adv_mom_5                       &
+                       +     7.0_wp * ibit16 * adv_mom_3                       &
+                       +              ibit15 * adv_mom_1                       &
+                          ) *                                                  &
+                             ( v(k+1,j,i)   + v(k,j,i)    )                    &
+                   -      (  8.0_wp * ibit17 * adv_mom_5                       &
+                       +              ibit16 * adv_mom_3                       &
+                          ) *                                                  &
+                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                    &
+                   +      (           ibit17 * adv_mom_5                       &
+                          ) *                                                  &
+                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                    &
                                       )
 
-                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
-                          ( 10.0_wp * ibit17 * adv_mom_5                        &
-                       +     3.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
-                          ) *                                                &
-                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
-                   -      (  5.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
-                          ) *                                                &
-                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
-                   +      (           ibit17 * adv_mom_5                        &
-                          ) *                                                &
-                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                   &
+                          ( 10.0_wp * ibit17 * adv_mom_5                       &
+                       +     3.0_wp * ibit16 * adv_mom_3                       &
+                       +              ibit15 * adv_mom_1                       &
+                          ) *                                                  &
+                             ( v(k+1,j,i)   - v(k,j,i)    )                    &
+                   -      (  5.0_wp * ibit17 * adv_mom_5                       &
+                       +              ibit16 * adv_mom_3                       &
+                          ) *                                                  &
+                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                    &
+                   +      (           ibit17 * adv_mom_5                       &
+                          ) *                                                  &
+                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                    &
                                                )
 !
@@ -5737,18 +5741,18 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
-                     +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
-                     +  (   w_comp                      * rho_air_zw(k) -     &
-                          ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)     &
-                        ) * drho_air(k) * ddzw(k)                             &
+                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx  &
+                     +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy  &
+                     +  (   w_comp                      * rho_air_zw(k) -      &
+                          ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)      &
+                        ) * drho_air(k) * ddzw(k)                              &
                       ) * 0.5_wp
  
-                tend(k,j,i) = tend(k,j,i) - (                                 &
-                       ( ( flux_r + diss_r )                                  &
-                     -   ( flux_l + diss_l ) ) * ddx                          &
-                     + ( ( flux_n + diss_n )                                  &
-                     -   ( flux_s + diss_s ) ) * ddy                          &
-                     + ( ( flux_t + diss_t )                                  &
-                     -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)        &
+                tend(k,j,i) = tend(k,j,i) - (                                  &
+                       ( ( flux_r + diss_r )                                   &
+                     -   ( flux_l + diss_l ) ) * ddx                           &
+                     + ( ( flux_n + diss_n )                                   &
+                     -   ( flux_s + diss_s ) ) * ddy                           &
+                     + ( ( flux_t + diss_t )                                   &
+                     -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)         &
                                             )  + v(k,j,i) * div
 
@@ -5896,14 +5900,14 @@
              u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
              swap_flux_x_local_w(k,j) = u_comp * (                             &
-                                      ( 37.0_wp * ibit20 * adv_mom_5              &
-                                   +     7.0_wp * ibit19 * adv_mom_3              &
-                                   +              ibit18 * adv_mom_1              &
+                                      ( 37.0_wp * ibit20 * adv_mom_5           &
+                                   +     7.0_wp * ibit19 * adv_mom_3           &
+                                   +              ibit18 * adv_mom_1           &
                                       ) *                                      &
                                      ( w(k,j,i)   + w(k,j,i-1) )               &
-                               -      (  8.0_wp * ibit20 * adv_mom_5              &
-                                   +              ibit19 * adv_mom_3              &
+                               -      (  8.0_wp * ibit20 * adv_mom_5           &
+                                   +              ibit19 * adv_mom_3           &
                                       ) *                                      &
                                      ( w(k,j,i+1) + w(k,j,i-2) )               &
-                               +      (           ibit20 * adv_mom_5              &
+                               +      (           ibit20 * adv_mom_5           &
                                       ) *                                      &
                                      ( w(k,j,i+2) + w(k,j,i-3) )               &
@@ -5911,14 +5915,14 @@
 
                swap_diss_x_local_w(k,j) = - ABS( u_comp ) * (                  &
-                                        ( 10.0_wp * ibit20 * adv_mom_5            &
-                                     +     3.0_wp * ibit19 * adv_mom_3            &
-                                     +              ibit18 * adv_mom_1            &
+                                        ( 10.0_wp * ibit20 * adv_mom_5         &
+                                     +     3.0_wp * ibit19 * adv_mom_3         &
+                                     +              ibit18 * adv_mom_1         &
                                         ) *                                    &
                                      ( w(k,j,i)   - w(k,j,i-1) )               &
-                                 -      (  5.0_wp * ibit20 * adv_mom_5            &
-                                     +              ibit19 * adv_mom_3            &
+                                 -      (  5.0_wp * ibit20 * adv_mom_5         &
+                                     +              ibit19 * adv_mom_3         &
                                         ) *                                    &
                                      ( w(k,j,i+1) - w(k,j,i-2) )               &
-                                 +      (           ibit20 * adv_mom_5            &
+                                 +      (           ibit20 * adv_mom_5         &
                                         ) *                                    &
                                      ( w(k,j,i+2) - w(k,j,i-3) )               &
@@ -5931,10 +5935,10 @@
              u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
              swap_flux_x_local_w(k,j) = u_comp * (                             &
-                            37.0_wp * ( w(k,j,i) + w(k,j,i-1)   )                 &
-                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )                 &
+                            37.0_wp * ( w(k,j,i) + w(k,j,i-1)   )              &
+                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )              &
                           +           ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
              swap_diss_x_local_w(k,j) = - ABS( u_comp ) * (                    &
-                            10.0_wp * ( w(k,j,i) - w(k,j,i-1)   )                 &
-                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )                 &
+                            10.0_wp * ( w(k,j,i) - w(k,j,i-1)   )              &
+                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )              &
                           +           ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
 
@@ -5972,32 +5976,32 @@
 
              v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
-             swap_flux_y_local_w(k) = v_comp * (                              &
-                                    ( 37.0_wp * ibit23 * adv_mom_5               &
-                                 +     7.0_wp * ibit22 * adv_mom_3               &
-                                 +              ibit21 * adv_mom_1               &
+             swap_flux_y_local_w(k) = v_comp * (                               &
+                                    ( 37.0_wp * ibit23 * adv_mom_5             &
+                                 +     7.0_wp * ibit22 * adv_mom_3             &
+                                 +              ibit21 * adv_mom_1             &
                                     ) *                                        &
-                                     ( w(k,j,i)   + w(k,j-1,i) )              &
-                             -      (  8.0_wp * ibit23 * adv_mom_5               &
-                                 +              ibit22 * adv_mom_3               &
-                                    ) *                                       &
-                                     ( w(k,j+1,i) + w(k,j-2,i) )              &
-                             +      (           ibit23 * adv_mom_5               &
-                                    ) *                                       &
-                                     ( w(k,j+2,i) + w(k,j-3,i) )              &
+                                     ( w(k,j,i)   + w(k,j-1,i) )               &
+                             -      (  8.0_wp * ibit23 * adv_mom_5             &
+                                 +              ibit22 * adv_mom_3             &
+                                    ) *                                        &
+                                     ( w(k,j+1,i) + w(k,j-2,i) )               &
+                             +      (           ibit23 * adv_mom_5             &
+                                    ) *                                        &
+                                     ( w(k,j+2,i) + w(k,j-3,i) )               &
                                                )
 
-             swap_diss_y_local_w(k) = - ABS( v_comp ) * (                     &
-                                    ( 10.0_wp * ibit23 * adv_mom_5               &
-                                 +     3.0_wp * ibit22 * adv_mom_3               &
-                                 +              ibit21 * adv_mom_1               &
-                                    ) *                                       &
-                                     ( w(k,j,i)   - w(k,j-1,i) )              &
-                             -      (  5.0_wp * ibit23 * adv_mom_5               &
-                                 +              ibit22 * adv_mom_3               &
-                                    ) *                                       &
-                                     ( w(k,j+1,i) - w(k,j-2,i) )              &
-                             +      (           ibit23 * adv_mom_5               &
-                                    ) *                                       &
-                                     ( w(k,j+2,i) - w(k,j-3,i) )              &
+             swap_diss_y_local_w(k) = - ABS( v_comp ) * (                      &
+                                    ( 10.0_wp * ibit23 * adv_mom_5             &
+                                 +     3.0_wp * ibit22 * adv_mom_3             &
+                                 +              ibit21 * adv_mom_1             &
+                                    ) *                                        &
+                                     ( w(k,j,i)   - w(k,j-1,i) )               &
+                             -      (  5.0_wp * ibit23 * adv_mom_5             &
+                                 +              ibit22 * adv_mom_3             &
+                                    ) *                                        &
+                                     ( w(k,j+1,i) - w(k,j-2,i) )               &
+                             +      (           ibit23 * adv_mom_5             &
+                                    ) *                                        &
+                                     ( w(k,j+2,i) - w(k,j-3,i) )               &
                                                         )
 
@@ -6007,11 +6011,11 @@
 
              v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
-             swap_flux_y_local_w(k) = v_comp * (                              &
-                           37.0_wp * ( w(k,j,i) + w(k,j-1,i)   )                 &
-                         -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
+             swap_flux_y_local_w(k) = v_comp * (                               &
+                           37.0_wp * ( w(k,j,i) + w(k,j-1,i)   )               &
+                         -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )               &
                          +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
-             swap_diss_y_local_w(k) = - ABS( v_comp ) * (                     &
-                           10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )                 &
-                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )                 &
+             swap_diss_y_local_w(k) = - ABS( v_comp ) * (                      &
+                           10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )               &
+                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )               &
                          +           ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
 
@@ -6027,6 +6031,8 @@
              k      = nzb + 1
              w_comp = w(k,j,i) + w(k-1,j,i)
-             flux_d = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
-             diss_d = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
+             flux_d = w_comp       * rho_air(k)                                &
+                    * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
+             diss_d = -ABS(w_comp) * rho_air(k)                                &
+                    * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
 
              DO  k = nzb+1, nzb_max_l
@@ -6037,32 +6043,32 @@
 
                 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
-                flux_r = u_comp * (                                          &
-                          ( 37.0_wp * ibit20 * adv_mom_5                        &
-                       +     7.0_wp * ibit19 * adv_mom_3                        &
-                       +              ibit18 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit20 * adv_mom_5                        &
-                       +              ibit19 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( w(k,j,i+2) + w(k,j,i-1) )                 &
-                   +      (           ibit20 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( w(k,j,i+3) + w(k,j,i-2) )                 &
+                flux_r = u_comp * (                                            &
+                          ( 37.0_wp * ibit20 * adv_mom_5                       &
+                       +     7.0_wp * ibit19 * adv_mom_3                       &
+                       +              ibit18 * adv_mom_1                       &
+                          ) *                                                  &
+                                 ( w(k,j,i+1) + w(k,j,i)   )                   &
+                   -      (  8.0_wp * ibit20 * adv_mom_5                       &
+                       +              ibit19 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( w(k,j,i+2) + w(k,j,i-1) )                   &
+                   +      (           ibit20 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( w(k,j,i+3) + w(k,j,i-2) )                   &
                                      )
 
-                diss_r = - ABS( u_comp ) * (                                 &
-                          ( 10.0_wp * ibit20 * adv_mom_5                        &
-                       +     3.0_wp * ibit19 * adv_mom_3                        &
-                       +              ibit18 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit20 * adv_mom_5                        &
-                       +              ibit19 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( w(k,j,i+2) - w(k,j,i-1) )                 &
-                   +      (           ibit20 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
+                diss_r = - ABS( u_comp ) * (                                   &
+                          ( 10.0_wp * ibit20 * adv_mom_5                       &
+                       +     3.0_wp * ibit19 * adv_mom_3                       &
+                       +              ibit18 * adv_mom_1                       &
+                          ) *                                                  &
+                                 ( w(k,j,i+1) - w(k,j,i)  )                    &
+                   -      (  5.0_wp * ibit20 * adv_mom_5                       &
+                       +              ibit19 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( w(k,j,i+2) - w(k,j,i-1) )                   &
+                   +      (           ibit20 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( w(k,j,i+3) - w(k,j,i-2) )                   &
                                               )
 
@@ -6075,32 +6081,32 @@
 
                 u_comp_l = u(k+1,j,i) + u(k,j,i) - gu
-                flux_l   = u_comp_l * (                                      &
-                                      ( 37.0_wp * ibit20_l * adv_mom_5          &
-                                   +     7.0_wp * ibit19_l * adv_mom_3          &
-                                   +              ibit18_l * adv_mom_1          &
-                                      ) *                                    &
-                                     ( w(k,j,i)   + w(k,j,i-1) )             &
-                               -      (  8.0_wp * ibit20_l * adv_mom_5          &
-                                   +              ibit19_l * adv_mom_3          &
-                                      ) *                                    &
-                                     ( w(k,j,i+1) + w(k,j,i-2) )             &
-                               +      (           ibit20_l * adv_mom_5          &
-                                      ) *                                    &
-                                     ( w(k,j,i+2) + w(k,j,i-3) )             &
+                flux_l   = u_comp_l * (                                        &
+                                      ( 37.0_wp * ibit20_l * adv_mom_5         &
+                                   +     7.0_wp * ibit19_l * adv_mom_3         &
+                                   +              ibit18_l * adv_mom_1         &
+                                      ) *                                      &
+                                     ( w(k,j,i)   + w(k,j,i-1) )               &
+                               -      (  8.0_wp * ibit20_l * adv_mom_5         &
+                                   +              ibit19_l * adv_mom_3         &
+                                      ) *                                      &
+                                     ( w(k,j,i+1) + w(k,j,i-2) )               &
+                               +      (           ibit20_l * adv_mom_5         &
+                                      ) *                                      &
+                                     ( w(k,j,i+2) + w(k,j,i-3) )               &
                                                  )
 
-                diss_l   = - ABS( u_comp_l ) * (                             &
-                                        ( 10.0_wp * ibit20_l * adv_mom_5        &
-                                     +     3.0_wp * ibit19_l * adv_mom_3        &
-                                     +              ibit18_l * adv_mom_1        &
-                                        ) *                                  &
-                                     ( w(k,j,i)   - w(k,j,i-1) )             &
-                                 -      (  5.0_wp * ibit20_l * adv_mom_5        &
-                                     +              ibit19_l * adv_mom_3        &
-                                        ) *                                  &
-                                     ( w(k,j,i+1) - w(k,j,i-2) )             &
-                                 +      (           ibit20_l * adv_mom_5        &
-                                        ) *                                  &
-                                     ( w(k,j,i+2) - w(k,j,i-3) )             &
+                diss_l   = - ABS( u_comp_l ) * (                               &
+                                        ( 10.0_wp * ibit20_l * adv_mom_5       &
+                                     +     3.0_wp * ibit19_l * adv_mom_3       &
+                                     +              ibit18_l * adv_mom_1       &
+                                        ) *                                    &
+                                     ( w(k,j,i)   - w(k,j,i-1) )               &
+                                 -      (  5.0_wp * ibit20_l * adv_mom_5       &
+                                     +              ibit19_l * adv_mom_3       &
+                                        ) *                                    &
+                                     ( w(k,j,i+1) - w(k,j,i-2) )               &
+                                 +      (           ibit20_l * adv_mom_5       &
+                                        ) *                                    &
+                                     ( w(k,j,i+2) - w(k,j,i-3) )               &
                                                             )
 #else
@@ -6115,32 +6121,32 @@
 
                 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
-                flux_n = v_comp * (                                          &
-                          ( 37.0_wp * ibit23 * adv_mom_5                        &
-                       +     7.0_wp * ibit22 * adv_mom_3                        &
-                       +              ibit21 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit23 * adv_mom_5                        &
-                       +              ibit22 * adv_mom_3                        &
-                          ) *                                                 &
-                                 ( w(k,j+2,i) + w(k,j-1,i) )                 &
-                   +      (           ibit23 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( w(k,j+3,i) + w(k,j-2,i) )                 &
+                flux_n = v_comp * (                                            &
+                          ( 37.0_wp * ibit23 * adv_mom_5                       &
+                       +     7.0_wp * ibit22 * adv_mom_3                       &
+                       +              ibit21 * adv_mom_1                       &
+                          ) *                                                  &
+                                 ( w(k,j+1,i) + w(k,j,i)   )                   &
+                   -      (  8.0_wp * ibit23 * adv_mom_5                       &
+                       +              ibit22 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( w(k,j+2,i) + w(k,j-1,i) )                   &
+                   +      (           ibit23 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( w(k,j+3,i) + w(k,j-2,i) )                   &
                                      )
 
-                diss_n = - ABS( v_comp ) * (                                 &
-                          ( 10.0_wp * ibit23 * adv_mom_5                        &
-                       +     3.0_wp * ibit22 * adv_mom_3                        &
-                       +              ibit21 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit23 * adv_mom_5                        &
-                       +              ibit22 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( w(k,j+2,i) - w(k,j-1,i) )                 &
-                   +      (           ibit23 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( w(k,j+3,i) - w(k,j-2,i) )                 &
+                diss_n = - ABS( v_comp ) * (                                   &
+                          ( 10.0_wp * ibit23 * adv_mom_5                       &
+                       +     3.0_wp * ibit22 * adv_mom_3                       &
+                       +              ibit21 * adv_mom_1                       &
+                          ) *                                                  &
+                                 ( w(k,j+1,i) - w(k,j,i)  )                    &
+                   -      (  5.0_wp * ibit23 * adv_mom_5                       &
+                       +              ibit22 * adv_mom_3                       &
+                          ) *                                                  &
+                                 ( w(k,j+2,i) - w(k,j-1,i) )                   &
+                   +      (           ibit23 * adv_mom_5                       &
+                          ) *                                                  &
+                                 ( w(k,j+3,i) - w(k,j-2,i) )                   &
                                               )
 
@@ -6153,32 +6159,32 @@
 
                 v_comp_s = v(k+1,j,i) + v(k,j,i) - gv
-                flux_s   = v_comp_s * (                                      &
-                                    ( 37.0_wp * ibit23_s * adv_mom_5            &
-                                 +     7.0_wp * ibit22_s * adv_mom_3            &
-                                 +              ibit21_s * adv_mom_1            &
-                                    ) *                                      &
-                                     ( w(k,j,i)   + w(k,j-1,i) )             &
-                             -      (  8.0_wp * ibit23_s * adv_mom_5            &
-                                 +              ibit22_s * adv_mom_3            &
-                                    ) *                                      &
-                                     ( w(k,j+1,i) + w(k,j-2,i) )             &
-                             +      (           ibit23_s * adv_mom_5            &
-                                    ) *                                      &
-                                     ( w(k,j+2,i) + w(k,j-3,i) )             &
+                flux_s   = v_comp_s * (                                        &
+                                    ( 37.0_wp * ibit23_s * adv_mom_5           &
+                                 +     7.0_wp * ibit22_s * adv_mom_3           &
+                                 +              ibit21_s * adv_mom_1           &
+                                    ) *                                        &
+                                     ( w(k,j,i)   + w(k,j-1,i) )               &
+                             -      (  8.0_wp * ibit23_s * adv_mom_5           &
+                                 +              ibit22_s * adv_mom_3           &
+                                    ) *                                        &
+                                     ( w(k,j+1,i) + w(k,j-2,i) )               &
+                             +      (           ibit23_s * adv_mom_5           &
+                                    ) *                                        &
+                                     ( w(k,j+2,i) + w(k,j-3,i) )               &
                                                )
 
-                diss_s   = - ABS( v_comp_s ) * (                             &
-                                    ( 10.0_wp * ibit23_s * adv_mom_5            &
-                                 +     3.0_wp * ibit22_s * adv_mom_3            &
-                                 +              ibit21_s * adv_mom_1            &
-                                    ) *                                      &
-                                     ( w(k,j,i)   - w(k,j-1,i) )             &
-                             -      (  5.0_wp * ibit23_s * adv_mom_5            &
-                                 +              ibit22_s * adv_mom_3            &
-                                    ) *                                      &
-                                     ( w(k,j+1,i) - w(k,j-2,i) )             &
-                             +      (           ibit23_s * adv_mom_5            &
-                                    ) *                                      &
-                                     ( w(k,j+2,i) - w(k,j-3,i) )             &
+                diss_s   = - ABS( v_comp_s ) * (                               &
+                                    ( 10.0_wp * ibit23_s * adv_mom_5           &
+                                 +     3.0_wp * ibit22_s * adv_mom_3           &
+                                 +              ibit21_s * adv_mom_1           &
+                                    ) *                                        &
+                                     ( w(k,j,i)   - w(k,j-1,i) )               &
+                             -      (  5.0_wp * ibit23_s * adv_mom_5           &
+                                 +              ibit22_s * adv_mom_3           &
+                                    ) *                                        &
+                                     ( w(k,j+1,i) - w(k,j-2,i) )               &
+                             +      (           ibit23_s * adv_mom_5           &
+                                    ) *                                        &
+                                     ( w(k,j+2,i) - w(k,j-3,i) )               &
                                                         )
 #else
@@ -6199,32 +6205,32 @@
 
                 w_comp = w(k+1,j,i) + w(k,j,i)
-                flux_t = w_comp * rho_air(k+1) * (                           &
-                          ( 37.0_wp * ibit26 * adv_mom_5                        &
-                       +     7.0_wp * ibit25 * adv_mom_3                        &
-                       +              ibit24 * adv_mom_1                        &
-                          ) *                                                &
-                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
-                   -      (  8.0_wp * ibit26 * adv_mom_5                        &
-                       +              ibit25 * adv_mom_3                        &
-                          ) *                                                &
-                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
-                   +      (           ibit26 * adv_mom_5                        &
-                          ) *                                                &
-                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
+                flux_t = w_comp * rho_air(k+1) * (                             &
+                          ( 37.0_wp * ibit26 * adv_mom_5                       &
+                       +     7.0_wp * ibit25 * adv_mom_3                       &
+                       +              ibit24 * adv_mom_1                       &
+                          ) *                                                  &
+                             ( w(k+1,j,i)  + w(k,j,i)     )                    &
+                   -      (  8.0_wp * ibit26 * adv_mom_5                       &
+                       +              ibit25 * adv_mom_3                       &
+                          ) *                                                  &
+                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                    &
+                   +      (           ibit26 * adv_mom_5                       &
+                          ) *                                                  &
+                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                    &
                                        )
 
-                diss_t = - ABS( w_comp ) * rho_air(k+1) * (                  &
-                          ( 10.0_wp * ibit26 * adv_mom_5                        &
-                       +     3.0_wp * ibit25 * adv_mom_3                        &
-                       +              ibit24 * adv_mom_1                        &
-                          ) *                                                &
-                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
-                   -      (  5.0_wp * ibit26 * adv_mom_5                        &
-                       +              ibit25 * adv_mom_3                        &
-                          ) *                                                &
-                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
-                   +      (           ibit26 * adv_mom_5                        &
-                          ) *                                                &
-                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
+                diss_t = - ABS( w_comp ) * rho_air(k+1) * (                    &
+                          ( 10.0_wp * ibit26 * adv_mom_5                       &
+                       +     3.0_wp * ibit25 * adv_mom_3                       &
+                       +              ibit24 * adv_mom_1                       &
+                          ) *                                                  &
+                             ( w(k+1,j,i)   - w(k,j,i)    )                    &
+                   -      (  5.0_wp * ibit26 * adv_mom_5                       &
+                       +              ibit25 * adv_mom_3                       &
+                          ) *                                                  &
+                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                    &
+                   +      (           ibit26 * adv_mom_5                       &
+                          ) *                                                  &
+                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                    &
                                                )
 !
