Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 2718)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 2731)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable loop vectorization by splitting the k-loop
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -1353,6 +1356,4 @@
        flux_t(0) = 0.0_wp
        diss_t(0) = 0.0_wp
-       flux_d    = 0.0_wp
-       diss_d    = 0.0_wp
 !        
 !--    Now compute the fluxes and tendency terms for the horizontal and
@@ -1473,4 +1474,10 @@
                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
                                          )
+       ENDDO
+
+       DO  k = nzb+1, nzb_max
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       Calculate the divergence of the velocity field. A respective
@@ -1513,6 +1520,4 @@
           swap_flux_x_local(k,j,tn) = flux_r(k)
           swap_diss_x_local(k,j,tn) = diss_r(k)
-          flux_d                    = flux_t(k)
-          diss_d                    = diss_t(k)
 
        ENDDO
@@ -1581,4 +1586,12 @@
                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
                                          )
+
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
+
 !
 !--       Calculate the divergence of the velocity field. A respective
@@ -1606,6 +1619,4 @@
           swap_flux_x_local(k,j,tn) = flux_r(k)
           swap_diss_x_local(k,j,tn) = diss_r(k)
-          flux_d                    = flux_t(k)
-          diss_d                    = diss_t(k)
 
        ENDDO
@@ -1765,6 +1776,4 @@
        REAL(wp)    ::  gv       !<
        REAL(wp)    ::  u_comp_l !<
-       REAL(wp)    ::  v_comp   !<
-       REAL(wp)    ::  w_comp   !<
        
        REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_n !<
@@ -1775,4 +1784,6 @@
        REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_t !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  u_comp !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_comp !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_comp !<
 
        gu = 2.0_wp * u_gtrans
@@ -1788,6 +1799,6 @@
              ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
 
-             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
-             flux_s_u(k,tn) = v_comp * (                                      &
+             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
+             flux_s_u(k,tn) = v_comp(k) * (                                   &
                             ( 37.0_wp * ibit14 * adv_mom_5                    &
                          +     7.0_wp * ibit13 * adv_mom_3                    &
@@ -1802,7 +1813,7 @@
                             ) *                                               &
                                         ( u(k,j+2,i) + u(k,j-3,i) )           &
-                                       )
-
-             diss_s_u(k,tn) = - ABS ( v_comp ) * (                            &
+                                          )
+
+             diss_s_u(k,tn) = - ABS ( v_comp(k) ) * (                         &
                             ( 10.0_wp * ibit14 * adv_mom_5                    &
                          +     3.0_wp * ibit13 * adv_mom_3                    &
@@ -1817,5 +1828,5 @@
                             ) *                                               &
                                         ( u(k,j+2,i) - u(k,j-3,i) )           &
-                                                 )
+                                                    )
 
           ENDDO
@@ -1823,10 +1834,10 @@
           DO  k = nzb_max+1, nzt
 
-             v_comp         = v(k,j,i) + v(k,j,i-1) - gv
-             flux_s_u(k,tn) = v_comp * (                                      &
+             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
+             flux_s_u(k,tn) = v_comp(k) * (                                   &
                            37.0_wp * ( u(k,j,i) + u(k,j-1,i)   )              &
                          -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )              &
                          +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
-             diss_s_u(k,tn) = - ABS(v_comp) * (                               &
+             diss_s_u(k,tn) = - ABS(v_comp(k)) * (                            &
                            10.0_wp * ( u(k,j,i) - u(k,j-1,i)   )              &
                          -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )              &
@@ -1897,6 +1908,5 @@
        flux_t(0) = 0.0_wp
        diss_t(0) = 0.0_wp
-       flux_d    = 0.0_wp
-       diss_d    = 0.0_wp
+       w_comp(0) = 0.0_wp
 !
 !--    Now compute the fluxes tendency terms for the horizontal and
@@ -1943,6 +1953,6 @@
           ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
 
-          v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-          flux_n(k) = v_comp * (                                              &
+          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
+          flux_n(k) = v_comp(k) * (                                           &
                      ( 37.0_wp * ibit14 * adv_mom_5                           &
                   +     7.0_wp * ibit13 * adv_mom_3                           &
@@ -1957,7 +1967,7 @@
                      ) *                                                      &
                                     ( u(k,j+3,i) + u(k,j-2,i) )               &
-                               )
-
-          diss_n(k) = - ABS ( v_comp ) * (                                    &
+                                  )
+
+          diss_n(k) = - ABS ( v_comp(k) ) * (                                 &
                      ( 10.0_wp * ibit14 * adv_mom_5                           &
                   +     3.0_wp * ibit13 * adv_mom_3                           &
@@ -1972,5 +1982,5 @@
                      ) *                                                      &
                                     ( u(k,j+3,i) - u(k,j-2,i) )               &
-                                         )
+                                            )
 !
 !--       k index has to be modified near bottom and top, else array
@@ -1984,6 +1994,6 @@
           k_mm  = k - 2 * ibit17
 
-          w_comp    = w(k,j,i) + w(k,j,i-1)
-          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
+          w_comp(k) = w(k,j,i) + w(k,j,i-1)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
                      ( 37.0_wp * ibit17 * adv_mom_5                           &
                   +     7.0_wp * ibit16 * adv_mom_3                           &
@@ -1998,7 +2008,7 @@
                      ) *                                                      &
                                 ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
-                                 )
-
-          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
                      ( 10.0_wp * ibit17 * adv_mom_5                           &
                   +     3.0_wp * ibit16 * adv_mom_3                           &
@@ -2013,5 +2023,11 @@
                      ) *                                                      &
                                 ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
-                                         )
+                                                           )
+       ENDDO
+
+       DO  k = nzb+1, nzb_max
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       Calculate the divergence of the velocity field. A respective
@@ -2025,5 +2041,5 @@
                                       )                                       &
                   ) * ddx                                                     &
-               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
+               +  ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 )         &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
                                     * ( IBITS(advc_flags_1(k,j-1,i),12,1)     &
@@ -2032,6 +2048,6 @@
                                       )                                       &
                   ) * ddy                                                     &
-               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
-                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
+               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
+                -   w_comp(k-1) * rho_air_zw(k-1)                             &
                                     * ( IBITS(advc_flags_1(k-1,j,i),15,1)     &
                                       + IBITS(advc_flags_1(k-1,j,i),16,1)     &
@@ -2056,6 +2072,4 @@
            flux_s_u(k,tn)   = flux_n(k)
            diss_s_u(k,tn)   = diss_n(k)
-           flux_d           = flux_t(k)
-           diss_d           = diss_t(k)
 !
 !--        Statistical Evaluation of u'u'. The factor has to be applied for
@@ -2073,9 +2087,9 @@
            sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
                 + ( flux_t(k)                                                  &
-                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
-                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                )     &
+                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        )     &
                   + diss_t(k)                                                  &
-                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
-                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           )     &
+                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                    )     &
                   ) *   weight_substep(intermediate_timestep_count)
        ENDDO
@@ -2093,10 +2107,10 @@
                        +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
 
-          v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-          flux_n(k) = v_comp * (                                              &
+          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
+          flux_n(k) = v_comp(k) * (                                           &
                          37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
                        -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
                        +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
-          diss_n(k) = - ABS( v_comp ) * (                                     &
+          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
                          10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
                        -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
@@ -2113,6 +2127,6 @@
           k_mm  = k - 2 * ibit17
 
-          w_comp    = w(k,j,i) + w(k,j,i-1)
-          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
+          w_comp(k) = w(k,j,i) + w(k,j,i-1)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
                      ( 37.0_wp * ibit17 * adv_mom_5                           &
                   +     7.0_wp * ibit16 * adv_mom_3                           &
@@ -2127,7 +2141,7 @@
                      ) *                                                      &
                                 ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
-                                 )
-
-          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
                      ( 10.0_wp * ibit17 * adv_mom_5                           &
                   +     3.0_wp * ibit16 * adv_mom_3                           &
@@ -2142,13 +2156,20 @@
                      ) *                                                      &
                                 ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
-                                         )
+                                                           )
+
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       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_comp(k)   - ( 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)           &
+          div = ( ( u_comp(k)      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx    &
+               +  ( v_comp(k) + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy    &
+               +  ( w_comp(k)   * rho_air_zw(k)   -                           &
+                    w_comp(k-1) * rho_air_zw(k-1)                             &
                   ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
@@ -2168,6 +2189,4 @@
            flux_s_u(k,tn)   = flux_n(k)
            diss_s_u(k,tn)   = diss_n(k)
-           flux_d           = flux_t(k)
-           diss_d           = diss_t(k)
 !
 !--        Statistical Evaluation of u'u'. The factor has to be applied for
@@ -2185,9 +2204,9 @@
            sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
                 + ( flux_t(k)                                                  &
-                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
-                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
+                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
                   + diss_t(k)                                                  &
-                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
-                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
+                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
                   ) *   weight_substep(intermediate_timestep_count)
        ENDDO
@@ -2254,7 +2273,5 @@
        REAL(wp)     ::  gu       !<
        REAL(wp)     ::  gv       !<
-       REAL(wp)     ::  u_comp   !<
        REAL(wp)     ::  v_comp_l !<
-       REAL(wp)     ::  w_comp   !<
        
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !<
@@ -2264,5 +2281,7 @@
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_r !<
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_t !<
+       REAL(wp), DIMENSION(nzb:nzt+1)  ::  u_comp !<
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !<
+       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !<
 
        gu = 2.0_wp * u_gtrans
@@ -2279,6 +2298,6 @@
              ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
 
-             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
-             flux_l_v(k,j,tn) = u_comp * (                                    &
+             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
+             flux_l_v(k,j,tn) = u_comp(k) * (                                 &
                               ( 37.0_wp * ibit20 * adv_mom_5                  &
                            +     7.0_wp * ibit19 * adv_mom_3                  &
@@ -2293,7 +2312,7 @@
                               ) *                                             &
                                         ( v(k,j,i+2) + v(k,j,i-3) )           &
-                                         )
-
-              diss_l_v(k,j,tn) = - ABS( u_comp ) * (                          &
+                                            )
+
+              diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                       &
                               ( 10.0_wp * ibit20 * adv_mom_5                  &
                            +     3.0_wp * ibit19 * adv_mom_3                  &
@@ -2308,5 +2327,5 @@
                               ) *                                             &
                                         ( v(k,j,i+2) - v(k,j,i-3) )           &
-                                                   )
+                                                      )
 
           ENDDO
@@ -2314,10 +2333,10 @@
           DO  k = nzb_max+1, nzt
 
-             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
-             flux_l_v(k,j,tn) = u_comp * (                                    &
+             u_comp(k)           = u(k,j-1,i) + u(k,j,i) - gu
+             flux_l_v(k,j,tn) = u_comp(k) * (                                 &
                              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_v(k,j,tn) = - ABS( u_comp ) * (                           &
+             diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                        &
                              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) )            &
@@ -2388,6 +2407,5 @@
        flux_t(0) = 0.0_wp
        diss_t(0) = 0.0_wp
-       flux_d    = 0.0_wp
-       diss_d    = 0.0_wp
+       w_comp(0) = 0.0_wp
 !
 !--    Now compute the fluxes and tendency terms for the horizontal and
@@ -2399,6 +2417,6 @@
           ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
  
-          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp * (                                              &
+          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
                      ( 37.0_wp * ibit20 * adv_mom_5                           &
                   +     7.0_wp * ibit19 * adv_mom_3                           &
@@ -2413,7 +2431,7 @@
                      ) *                                                      &
                                     ( v(k,j,i+3) + v(k,j,i-2) )               &
-                               )
-
-          diss_r(k) = - ABS( u_comp ) * (                                     &
+                                  )
+
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
                      ( 10.0_wp * ibit20 * adv_mom_5                           &
                   +     3.0_wp * ibit19 * adv_mom_3                           &
@@ -2428,5 +2446,5 @@
                      ) *                                                      &
                                     ( v(k,j,i+3) - v(k,j,i-2) )               &
-                                        )
+                                           )
 
           ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
@@ -2476,6 +2494,6 @@
           k_mm  = k - 2 * ibit26
 
-          w_comp    = w(k,j-1,i) + w(k,j,i)
-          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
+          w_comp(k) = w(k,j-1,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
                      ( 37.0_wp * ibit26 * adv_mom_5                           &
                   +     7.0_wp * ibit25 * adv_mom_3                           &
@@ -2490,7 +2508,7 @@
                      ) *                                                      &
                                 ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
-                                 )
-
-          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
                      ( 10.0_wp * ibit26 * adv_mom_5                           &
                   +     3.0_wp * ibit25 * adv_mom_3                           &
@@ -2505,10 +2523,16 @@
                      ) *                                                      &
                                 ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
-                                         )
+                                                           )
+       ENDDO
+
+       DO  k = nzb+1, nzb_max
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       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_comp     + gu )                                       &
+          div = ( ( ( u_comp(k)     + gu )                                    &
                                        * ( ibit18 + ibit19 + ibit20 )         &
                   - ( u(k,j-1,i) + u(k,j,i) )                                 &
@@ -2526,6 +2550,6 @@
                                          )                                    &
                   ) * ddy                                                     &
-               +  ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )     &
-                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
+               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&
+                -   w_comp(k-1) * rho_air_zw(k-1)                             &
                                        * ( IBITS(advc_flags_1(k-1,j,i),24,1)  &
                                          + IBITS(advc_flags_1(k-1,j,i),25,1)  &
@@ -2550,6 +2574,4 @@
            flux_s_v(k,tn)   = flux_n(k)
            diss_s_v(k,tn)   = diss_n(k)
-           flux_d           = flux_t(k)
-           diss_d           = diss_t(k)
 
 !
@@ -2568,9 +2590,9 @@
            sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
                 + ( flux_t(k)                                                  &
-                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
-                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
+                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
                   + diss_t(k)                                                  &
-                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
-                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
+                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
                   ) *   weight_substep(intermediate_timestep_count)
 
@@ -2579,11 +2601,11 @@
        DO  k = nzb_max+1, nzt
 
-          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp * (                                              &
+          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
                       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(k) = - ABS( u_comp ) * (                                     &
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
                       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) )                   &
@@ -2612,6 +2634,6 @@
           k_mm  = k - 2 * ibit26
 
-          w_comp    = w(k,j-1,i) + w(k,j,i)
-          flux_t(k) = w_comp * rho_air_zw(k) * (                              &
+          w_comp(k) = w(k,j-1,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
                      ( 37.0_wp * ibit26 * adv_mom_5                           &
                   +     7.0_wp * ibit25 * adv_mom_3                           &
@@ -2626,7 +2648,7 @@
                      ) *                                                      &
                                 ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
-                                 )
-
-          diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
                      ( 10.0_wp * ibit26 * adv_mom_5                           &
                   +     3.0_wp * ibit25 * adv_mom_3                           &
@@ -2641,13 +2663,19 @@
                      ) *                                                      &
                                 ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
-                                         )
+                                                           )
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       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_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
-               +  ( v_comp(k)   - ( 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)           &
+          div = ( ( u_comp(k) + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx    &
+               +  ( v_comp(k)      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy    &
+               +  ( w_comp(k)   * rho_air_zw(k)   -                           &
+                    w_comp(k-1) * rho_air_zw(k-1)                             &
                   ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
@@ -2667,7 +2695,4 @@
            flux_s_v(k,tn)   = flux_n(k)
            diss_s_v(k,tn)   = diss_n(k)
-           flux_d           = flux_t(k)
-           diss_d           = diss_t(k)
-
 !
 !--        Statistical Evaluation of v'v'. The factor has to be applied for
@@ -2685,9 +2710,9 @@
            sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
                 + ( flux_t(k)                                                  &
-                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
-                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
+                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
                   + diss_t(k)                                                  &
-                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
-                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
+                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
                   ) *   weight_substep(intermediate_timestep_count)
 
@@ -2754,7 +2779,4 @@
        REAL(wp)    ::  gu      !<
        REAL(wp)    ::  gv      !<
-       REAL(wp)    ::  u_comp  !<
-       REAL(wp)    ::  v_comp  !<
-       REAL(wp)    ::  w_comp  !<
        
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !<
@@ -2764,4 +2786,7 @@
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_r !<
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_t !<
+       REAL(wp), DIMENSION(nzb:nzt+1)  ::  u_comp !<
+       REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !<
+       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !<
 
        gu = 2.0_wp * u_gtrans
@@ -2777,6 +2802,6 @@
              ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
 
-             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
-             flux_s_w(k,tn) = v_comp * (                                      &
+             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
+             flux_s_w(k,tn) = v_comp(k) * (                                   &
                             ( 37.0_wp * ibit32 * adv_mom_5                    &
                          +     7.0_wp * ibit31 * adv_mom_3                    &
@@ -2791,7 +2816,7 @@
                             ) *                                               &
                                         ( w(k,j+2,i) + w(k,j-3,i) )           &
-                                       )
-
-             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
+                                          )
+
+             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
                             ( 10.0_wp * ibit32 * adv_mom_5                    &
                          +     3.0_wp * ibit31 * adv_mom_3                    &
@@ -2806,5 +2831,5 @@
                             ) *                                               &
                                         ( w(k,j+2,i) - w(k,j-3,i) )           &
-                                                )
+                                                   )
 
           ENDDO
@@ -2812,10 +2837,10 @@
           DO  k = nzb_max+1, nzt
 
-             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
-             flux_s_w(k,tn) = v_comp * (                                      &
+             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
+             flux_s_w(k,tn) = v_comp(k) * (                                   &
                            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
-             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
+             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                         &
                            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) )              &
@@ -2835,6 +2860,6 @@
              ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
 
-             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
-             flux_l_w(k,j,tn) = u_comp * (                                    &
+             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
+             flux_l_w(k,j,tn) = u_comp(k) * (                                 &
                              ( 37.0_wp * ibit29 * adv_mom_5                   &
                           +     7.0_wp * ibit28 * adv_mom_3                   &
@@ -2849,7 +2874,7 @@
                              ) *                                              &
                                         ( w(k,j,i+2) + w(k,j,i-3) )           &
-                                         )
-
-               diss_l_w(k,j,tn) = - ABS( u_comp ) * (                         &
+                                            )
+
+               diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                      &
                              ( 10.0_wp * ibit29 * adv_mom_5                   &
                           +     3.0_wp * ibit28 * adv_mom_3                   &
@@ -2864,5 +2889,5 @@
                              ) *                                              &
                                         ( w(k,j,i+2) - w(k,j,i-3) )           &
-                                                    )
+                                                       )
 
           ENDDO
@@ -2870,10 +2895,10 @@
           DO  k = nzb_max+1, nzt
 
-             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
-             flux_l_w(k,j,tn) = u_comp * (                                    &
+             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
+             flux_l_w(k,j,tn) = u_comp(k) * (                                 &
                             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
-             diss_l_w(k,j,tn) = - ABS( u_comp ) * (                           &
+             diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                        &
                             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) )             &
@@ -2888,9 +2913,7 @@
 !--    advc_flags_1.
        k         = nzb + 1
-       w_comp    = w(k,j,i) + w(k-1,j,i)
-       flux_t(0) = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
-       diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
-       flux_d    = flux_t(0)
-       diss_d    = diss_t(0)
+       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
 !
 !--    Now compute the fluxes and tendency terms for the horizontal
@@ -2902,6 +2925,6 @@
           ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
 
-          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp * (                                              &
+          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
                      ( 37.0_wp * ibit29 * adv_mom_5                           &
                   +     7.0_wp * ibit28 * adv_mom_3                           &
@@ -2916,7 +2939,7 @@
                      ) *                                                      &
                                     ( w(k,j,i+3) + w(k,j,i-2) )               &
-                               )
-
-          diss_r(k) = - ABS( u_comp ) * (                                     &
+                                  )
+
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
                      ( 10.0_wp * ibit29 * adv_mom_5                           &
                   +     3.0_wp * ibit28 * adv_mom_3                           &
@@ -2931,5 +2954,5 @@
                      ) *                                                      &
                                     ( w(k,j,i+3) - w(k,j,i-2) )               &
-                                        )
+                                           )
 
           ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
@@ -2937,6 +2960,6 @@
           ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
 
-          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
-          flux_n(k) = v_comp * (                                              &
+          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
+          flux_n(k) = v_comp(k) * (                                           &
                      ( 37.0_wp * ibit32 * adv_mom_5                           &
                   +     7.0_wp * ibit31 * adv_mom_3                           &
@@ -2951,7 +2974,7 @@
                      ) *                                                      &
                                     ( w(k,j+3,i) + w(k,j-2,i) )               &
-                               )
-
-          diss_n(k) = - ABS( v_comp ) * (                                     &
+                                  )
+
+          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
                      ( 10.0_wp * ibit32 * adv_mom_5                           &
                   +     3.0_wp * ibit31 * adv_mom_3                           &
@@ -2966,5 +2989,5 @@
                      ) *                                                      &
                                    ( w(k,j+3,i) - w(k,j-2,i) )                &
-                                        )
+                                           )
 !
 !--       k index has to be modified near bottom and top, else array
@@ -2978,6 +3001,6 @@
           k_mm  = k - 2 * ibit35
 
-          w_comp    = w(k+1,j,i) + w(k,j,i)
-          flux_t(k) = w_comp * rho_air(k+1) * (                               &
+          w_comp(k) = w(k+1,j,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
                      ( 37.0_wp * ibit35 * adv_mom_5                           &
                   +     7.0_wp * ibit34 * adv_mom_3                           &
@@ -2992,7 +3015,7 @@
                      ) *                                                      &
                                 ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
-                                )
-
-          diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
+                                                 )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
                      ( 10.0_wp * ibit35 * adv_mom_5                           &
                   +     3.0_wp * ibit34 * adv_mom_3                           &
@@ -3007,11 +3030,16 @@
                      ) *                                                      &
                                 ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
-                                        )
-
+                                                          )
+       ENDDO
+
+       DO  k = nzb+1, nzb_max
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       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_comp + gu ) * ( ibit27 + ibit28 + ibit29 )            &
+          div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 )         &
                   - ( u(k+1,j,i) + u(k,j,i)   )                               & 
                                     * ( IBITS(advc_flags_1(k,j,i-1),27,1)     &
@@ -3020,5 +3048,5 @@
                                       )                                       &
                   ) * ddx                                                     &
-              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
+              +   ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 )         & 
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
                                     * ( IBITS(advc_flags_1(k,j-1,i),30,1)     &
@@ -3027,6 +3055,7 @@
                                       )                                       &
                   ) * ddy                                                     &
-              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
-                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
+              +   ( w_comp(k)               * rho_air(k+1)                    &
+                                            * ( ibit33 + ibit34 + ibit35 )    &
+                - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
                                     * ( IBITS(advc_flags_2(k-1,j,i),1,1)      &
                                       + IBITS(advc_flags_2(k-1,j,i),2,1)      &
@@ -3051,15 +3080,13 @@
           flux_s_w(k,tn)   = flux_n(k)
           diss_s_w(k,tn)   = diss_n(k)
-          flux_d           = flux_t(k)
-          diss_d           = diss_t(k)
 !
 !--       Statistical Evaluation of w'w'.
           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
                       + ( flux_t(k)                                           &
-                       * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
-                       / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
+                       * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                ) &
+                       / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        ) &
                         + diss_t(k)                                           &
-                       *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
-                       / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
+                       *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           ) &
+                       / ( ABS( w_comp(k) ) + 1.0E-20_wp                    ) &
                         ) *   weight_substep(intermediate_timestep_count)
 
@@ -3068,22 +3095,22 @@
        DO  k = nzb_max+1, nzt
 
-          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp * (                                              &
+          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
                       37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
                     -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
                     +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
 
-          diss_r(k) = - ABS( u_comp ) * (                                     &
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
                       10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
                     -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
                     +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
 
-          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
-          flux_n(k) = v_comp * (                                              &
+          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
+          flux_n(k) = v_comp(k) * (                                           &
                       37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
                     -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
                     +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
 
-          diss_n(k) = - ABS( v_comp ) * (                                     &
+          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
                       10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
                     -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
@@ -3100,6 +3127,6 @@
           k_mm  = k - 2 * ibit35
 
-          w_comp    = w(k+1,j,i) + w(k,j,i)
-          flux_t(k) = w_comp * rho_air(k+1) * (                               &
+          w_comp(k) = w(k+1,j,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
                      ( 37.0_wp * ibit35 * adv_mom_5                           &
                   +     7.0_wp * ibit34 * adv_mom_3                           &
@@ -3114,7 +3141,7 @@
                      ) *                                                      &
                                 ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
-                                )
-
-          diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
+                                                 )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
                      ( 10.0_wp * ibit35 * adv_mom_5                           &
                   +     3.0_wp * ibit34 * adv_mom_3                           &
@@ -3129,13 +3156,19 @@
                      ) *                                                      &
                                 ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
-                                        )
+                                                          )
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          flux_d    = flux_t(k-1)
+          diss_d    = diss_t(k-1)
 !
 !--       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_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
-              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
-              +   (   w_comp                    * rho_air(k+1) -              &
-                    ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)                  &
+          div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx      &
+              +   ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy      &
+              +   ( w_comp(k)               * rho_air(k+1) -                  &
+                  ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
                   ) * drho_air_zw(k) * ddzu(k+1)                              &
                 ) * 0.5_wp
@@ -3155,15 +3188,13 @@
           flux_s_w(k,tn)   = flux_n(k)
           diss_s_w(k,tn)   = diss_n(k)
-          flux_d           = flux_t(k)
-          diss_d           = diss_t(k)
 !
 !--       Statistical Evaluation of w'w'.
           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
                       + ( flux_t(k)                                           &
-                       * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
-                       / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
+                       * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                ) &
+                       / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        ) &
                         + diss_t(k)                                           &
-                       *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
-                       / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
+                       *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           ) &
+                       / ( ABS( w_comp(k) ) + 1.0E-20_wp                    ) &
                         ) *   weight_substep(intermediate_timestep_count)
 
