Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 712)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 713)
@@ -4,5 +4,9 @@
 ! Current revisions:
 ! -----------------
-! 
+! File reformatted.
+! Bugfix in vertical advection of w concerning the optimized version for   
+! vector architecture.
+! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
+! broken numbers.
 !
 ! Former revisions:
@@ -13,5 +17,5 @@
 ! formatting adjustments
 !
-! 705 2011-03-25 11:21:43 Z suehring
+! 705 2011-03-25 11:21:43 Z suehring $
 ! Bugfix in declaration of logicals concerning outflow boundaries.
 ! 
@@ -111,8 +115,8 @@
 !
 !--    Set the appropriate factors for scalar and momentum advection.
-       adv_sca_5 = 0.016666666666666
-       adv_sca_3 = 0.083333333333333
-       adv_mom_5 = 0.0083333333333333
-       adv_mom_3 = 0.041666666666666
+       adv_sca_5 = 1./60.
+       adv_sca_3 = 1./12.
+       adv_mom_5 = 1./120.
+       adv_mom_3 = 1./24.
 
 !         
@@ -317,5 +321,5 @@
              CASE ( 3 )
 
-                DO  k = nzb_s_inner(j,i) + 1, nzt
+                DO  k = nzb_s_inner(j,i)+1, nzt
                    v_comp = v(k,j+1,i) - v_gtrans
                    flux_n(k) = v_comp * (                                      &
@@ -473,5 +477,5 @@
        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
        
-           DO  k = nzb_s_inner(j,i) + 1, nzt
+           DO  k = nzb_s_inner(j,i)+1, nzt
               v_comp               = v(k,j,i) - v_gtrans
               swap_flux_y_local(k) = v_comp * (                               &
@@ -534,5 +538,5 @@
 !       
 !--    2nd-order scheme (bottom)
-       k         = nzb_s_inner(j,i) + 1
+       k         = nzb_s_inner(j,i)+1
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
@@ -684,18 +688,19 @@
          
             CASE ( 1 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
                   flux_r(k) = ( u_comp(k) - gu ) * (                          &
-                              7. * (u(k,j,i+1) + u(k,j,i)    )                &
-                              -    ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
-                  diss_r(k) = - abs(u_comp(k) - gu) * (                       &
-                              3. * (u(k,j,i+1) - u(k,j,i)    )                & 
-                              -    ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
+                              7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
+                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
+                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                     &
+                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               & 
+                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
               ENDDO
               
             CASE ( 2 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                  flux_r(k) = (u_comp(k) - gu) * ( u(k,j,i+1) + u(k,j,i) )*0.25 
+                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
+                                u(k,j,i+1) + u(k,j,i) ) * 0.25 
                   diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i),     &
                                         u(k,j,i-1), u(k,j,i-2), u_comp(k),    &
@@ -704,16 +709,16 @@
                
             CASE ( 3 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                   flux_n(k) = v_comp * (                                      &
-                              7. * ( u(k,j+1,i) + u(k,j,i) )                  &
-                              -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( u(k,j+1,i) - u(k,j,i)   )                & 
-                              -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
+                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
+                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
+                  diss_n(k) = - ABS( v_comp ) * (                             &
+                              3.0 * ( u(k,j+1,i) - u(k,j,i)     )             & 
+                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
                ENDDO
                
             CASE ( 4 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                   flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
@@ -724,36 +729,36 @@
                
             CASE ( 5 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
 !               
 !--               Compute leftside fluxes for the left boundary of PE domain
                   u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
                   flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
-                  diss_l_u(k,j) = diss_2nd(u(k,j,i+2), u(k,j,i+1), u(k,j,i),  &
-                                           u(k,j,i-1), u(k,j,i-1), u_comp(k), &
-                                           0.25, ddx )
+                  diss_l_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
+                                            u(k,j,i-1), u(k,j,i-1), u_comp(k),&
+                                            0.25, ddx )
                  
                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
                   flux_r(k) = ( u_comp(k) - gu ) * (                          &
-                              7. * (u(k,j,i+1) + u(k,j,i)    )                &
-                              -    ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
-                  diss_r(k) = - abs( u_comp(k) -gu ) * (                      &
-                              3. * ( u(k,j,i+1) - u(k,j,i)     )              & 
-                              -    ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
+                              7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
+                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
+                  diss_r(k) = - ABS( u_comp(k) -gu ) * (                      &
+                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               & 
+                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
                ENDDO
                degraded_l = .TRUE.
                
             CASE ( 7 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt                           
+               DO  k = nzb_u_inner(j,i)+1, nzt                           
                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                   flux_n(k) = v_comp * (                                      &
-                              7. * ( u(k,j+1,i) + u(k,j,i)   )                &
-                              -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( u(k,j+1,i) - u(k,j,i)   )                & 
-                              -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
+                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
+                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
+                  diss_n(k) = - ABS( v_comp ) * (                             &
+                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               & 
+                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
                ENDDO 
                
             CASE ( 8 )
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
 !               
 !--              Compute southside fluxes for the south boundary of PE domain           
@@ -766,9 +771,9 @@
                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                   flux_n(k) = v_comp * (                                      &
-                              7. * ( u(k,j+1,i) + u(k,j,i)   )                &
-                              -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( u(k,j+1,i) - u(k,j,i)   )                & 
-                              -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
+                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
+                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
+                  diss_n(k) = - ABS( v_comp ) * (                             &
+                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               & 
+                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
                ENDDO  
                degraded_s = .TRUE.
@@ -779,31 +784,31 @@
 !         
 !--      Compute the crosswise 5th order fluxes at the outflow
-         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
-         .OR. boundary_flags(j,i) == 5 )  THEN
+         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
+              boundary_flags(j,i) == 5 )  THEN
          
              DO  k = nzb_u_inner(j,i)+1, nzt
                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                flux_n(k) = v_comp * (                                         &
-                           37.  * ( u(k,j+1,i) + u(k,j,i)   )                 &
-                           - 8. * ( 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) * (                                  &
-                           10.  * ( u(k,j+1,i) - u(k,j,i)   )                 &
-                           - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                 &
-                           +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
+                           37.0 * ( u(k,j+1,i) + u(k,j,i)   )                 &
+                         -  8.0 * ( 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 ) * (                               &
+                           10.0 * ( u(k,j+1,i) - u(k,j,i)   )                 &
+                         -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                 &
+                         +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
              ENDDO
              
          ELSE
          
-            DO  k = nzb_u_inner(j,i) + 1, nzt
+            DO  k = nzb_u_inner(j,i)+1, nzt
                u_comp(k) = u(k,j,i+1) + u(k,j,i)
                flux_r(k) = ( u_comp(k) - gu ) * (                             &
-                           37.  * ( u(k,j,i+1) + u(k,j,i)   )                 &
-                           - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                 &
-                           +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-               diss_r(k) = - abs(u_comp(k) - gu) * (                          &
-                           10.  * ( u(k,j,i+1) - u(k,j,i) )                   &
-                           - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                 &
-                           +      ( u(k,j,i+3) - u(k,j,i-2) ) ) *adv_mom_5
+                           37.0 * ( u(k,j,i+1) + u(k,j,i)   )                 &
+                         -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                 &
+                         +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
+               diss_r(k) = - ABS( u_comp(k) - gu ) * (                        &
+                           10.0 * ( u(k,j,i+1) - u(k,j,i) )                   &
+                         -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                 &
+                         +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
             ENDDO
             
@@ -813,24 +818,24 @@
 !       
 !--       Compute the fifth order fluxes for the interior of PE domain.                  
-          DO  k = nzb_u_inner(j,i) + 1, nzt
+          DO  k = nzb_u_inner(j,i)+1, nzt
              u_comp(k) = u(k,j,i+1) + u(k,j,i)
              flux_r(k) = ( u_comp(k) - gu ) * (                               &
-                         37.  * ( u(k,j,i+1) + u(k,j,i)   )                   &
-                         - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                   &
-                         +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-             diss_r(k) = - abs(u_comp(k) - gu) * (                            &
-                         10.  * ( u(k,j,i+1) - u(k,j,i)   )                   &
-                         - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                   &
-                         +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
+                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
+                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
+                       +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
+             diss_r(k) = - ABS( u_comp(k) - gu ) * (                          &
+                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
+                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
+                       +        ( 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 * (                                           &
-                         37.  * ( u(k,j+1,i) + u(k,j,i)   )                   &
-                         - 8. * ( 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) * (                                    &
-                         10.  * ( u(k,j+1,i) - u(k,j,i)   )                   &
-                         - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                   &
-                         +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
+                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
+                       -  8.0 * ( 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 ) * (                                  &
+                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
+                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
+                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
           ENDDO 
          
@@ -840,14 +845,14 @@
        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
        
-          DO  k = nzb_u_inner(j,i) + 1, nzt
+          DO  k = nzb_u_inner(j,i)+1, nzt
              v_comp      = v(k,j,i) + v(k,j,i-1) - gv
              flux_s_u(k) = v_comp * (                                         &
-                           37.  * ( u(k,j,i) + u(k,j-1,i)   )                 &
-                           - 8. * ( 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) = - abs(v_comp) * (                                  &
-                           10.  * ( u(k,j,i) - u(k,j-1,i)    )                &
-                           - 5. * ( u(k,j+1,i) - u(k,j-2,i)  )                &
-                           +      (  u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
+                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
+                         -  8.0 * ( 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) = - ABS(v_comp) * (                                  &
+                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
+                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
+                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
           ENDDO
           
@@ -859,11 +864,11 @@
              u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
              flux_l_u(k,j) = u_comp_l * (                                     &
-                             37.  * ( u(k,j,i) + u(k,j,i-1)   )               &
-                             - 8. * ( u(k,j,i+1) + u(k,j,i-2) )               &
-                             +      ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
-             diss_l_u(k,j) = - abs(u_comp_l) * (                              &
-                             10.  * ( u(k,j,i) - u(k,j,i-1)   )               &
-                             - 5. * ( u(k,j,i+1) - u(k,j,i-2) )               &
-                             +      ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
+                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )               &
+                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )               &
+                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
+             diss_l_u(k,j) = - ABS(u_comp_l) * (                              &
+                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )               &
+                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )               &
+                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
           ENDDO
           
@@ -871,10 +876,10 @@
 !       
 !--    Now compute the tendency terms for the horizontal parts.
-       DO  k = nzb_u_inner(j,i) + 1, nzt                   
+       DO  k = nzb_u_inner(j,i)+1, nzt                   
           tend(k,j,i) = tend(k,j,i) - (                                       &
                             ( flux_r(k) + diss_r(k)                           &
-                          - ( flux_l_u(k,j) + diss_l_u(k,j) ) ) * ddx         &
+                          -   flux_l_u(k,j) - diss_l_u(k,j) ) * ddx         &
                           + ( flux_n(k) + diss_n(k)                           &
-                          - ( flux_s_u(k) + diss_s_u(k) )     ) * ddy  )
+                          -   flux_s_u(k) - diss_s_u(k)     ) * ddy  )
 
            flux_l_u(k,j) = flux_r(k)
@@ -887,12 +892,13 @@
            sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:)                            &
              + ( flux_r(k) *                                                  &
-             ( u_comp(k) - 2. * hom(k,1,1,:) ) / ( u_comp(k) - gu + 1.0E-20 ) &
-             + diss_r(k)                                                      &
-             * abs(u_comp(k) - 2. * hom(k,1,1,:) )                            & 
-             / (abs(u_comp(k) - gu) + 1.0E-20)         )                      &
-             * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+               ( u_comp(k) - 2.0 * hom(k,1,1,:) )                             &
+             / ( u_comp(k) - gu + 1.0E-20      )                              &
+             +   diss_r(k) *                                                  &
+                 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                        & 
+             / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )                          &
+             *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
        ENDDO 
-       sums_us2_ws_l(nzb_u_inner(j,i),:) =                                   &
-                                           sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
+       sums_us2_ws_l(nzb_u_inner(j,i),:) = sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
+                                           
 
 !
@@ -903,6 +909,6 @@
        diss_t(nzb_u_inner(j,i)) = 0.
 !
-!--    2nd order scheme
-       k         = nzb_u_inner(j,i) + 1
+!--    2nd order scheme (bottom)
+       k         = nzb_u_inner(j,i)+1
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
@@ -913,39 +919,40 @@
              
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!--    WS3 as an intermediate step
-       k         = nzb_u_inner(j,i) + 2
+                                 -   flux_d    - diss_d      ) * ddzw(k)
+!
+!--    WS3 as an intermediate step (bottom)
+       k         = nzb_u_inner(j,i)+2
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
        w_comp    = w(k,j,i) + w(k,j,i-1)
-       flux_t(k) = w_comp*(                                                   &
-                   7. * ( u(k+1,j,i) + u(k,j,i)   )                           &
-                   -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
-       diss_t(k) = -abs(w_comp)*(                                             &
-                   3. * ( u(k+1,j,i) - u(k,j,i)   )                           & 
-                   -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
+       flux_t(k) = w_comp * (                                                 &
+                   7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
+                 -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
+       diss_t(k) = -ABS( w_comp) * (                                          &
+                   3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          & 
+                 -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
 
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
-
-       DO  k = nzb_u_inner(j,i) + 3, nzt - 2
+                                 -   flux_d    - diss_d      ) * ddzw(k)
+!
+!--    WS5
+       DO  k = nzb_u_inner(j,i)+3, nzt-2
           flux_d    = flux_t(k-1)
           diss_d    = diss_t(k-1)
           w_comp    = w(k,j,i) + w(k,j,i-1)
-          flux_t(k) = w_comp*(                                                &
-                      37.*   ( u(k+1,j,i) + u(k,j,i)   )                      &
-                      - 8. * ( u(k+2,j,i) + u(k-1,j,i) )                      &
-                      +      ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
-          diss_t(k) = - abs(w_comp) * (                                       &
-                      10. *  ( u(k+1,j,i) - u(k,j,i)   )                      & 
-                      - 5. * ( u(k+2,j,i) - u(k-1,j,i) )                      &
-                      +      ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
+          flux_t(k) = w_comp * (                                              &
+                      37.0 * ( u(k+1,j,i) + u(k,j,i)   )                      &
+                    -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                      &
+                    +        ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
+          diss_t(k) = - ABS( w_comp ) * (                                     &
+                      10.0 *  ( u(k+1,j,i) - u(k,j,i)   )                     & 
+                    -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                      &
+                    +        ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
 
           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                    -   flux_d    - diss_d      ) * ddzw(k)
        ENDDO
 !
-!--    WS3 as an intermediate step
+!--    WS3 as an intermediate step (top)
        k         = nzt - 1
        flux_d    = flux_t(k-1)
@@ -953,15 +960,15 @@
        w_comp    = w(k,j,i) + w(k,j,i-1)
        flux_t(k) = w_comp * (                                                 &
-                   7. * ( u(k+1,j,i) + u(k,j,i)   )                           &
-                   -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
-       diss_t(k) = - abs(w_comp) * (                                          &
-                   3. * ( u(k+1,j,i) - u(k,j,i)   )                           & 
-                   -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
+                   7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
+                 -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
+       diss_t(k) = - ABS( w_comp ) * (                                        &
+                   3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          & 
+                 -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
                    
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                 -   flux_d    - diss_d      ) * ddzw(k)
        
 !
-!--    2nd order scheme
+!--    2nd order scheme (top)
        k         = nzt
        flux_d    = flux_t(k-1)
@@ -973,5 +980,5 @@
 
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                 -   flux_d    - diss_d      ) * ddzw(k)
 
 !
@@ -980,5 +987,5 @@
           sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:)                           &
               + ( flux_t(k) + diss_t(k) )                                     &
-              * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+              *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
        ENDDO
 
@@ -1021,16 +1028,16 @@
           
              CASE ( 1 )
-                DO  k = nzb_v_inner(j,i) + 1, nzt
+                DO  k = nzb_v_inner(j,i)+1, nzt
                    u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                    flux_r(k) = u_comp * (                                     &
-                               7. * (v(k,j,i+1) + v(k,j,i)    )               &
-                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
-                   diss_r(k) = - abs(u_comp) * (                              &
-                               3. * ( v(k,j,i+1) - v(k,j,i)   )               & 
-                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
+                               7.0 * (v(k,j,i+1) + v(k,j,i)    )              &
+                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
+                   diss_r(k) = - ABS( u_comp ) * (                            &
+                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              & 
+                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
                 ENDDO
                 
              CASE ( 2 )
-                DO  k = nzb_v_inner(j,i) + 1, nzt
+                DO  k = nzb_v_inner(j,i)+1, nzt
                    u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                    flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 
@@ -1041,16 +1048,16 @@
                 
              CASE ( 3 )
-                DO  k = nzb_v_inner(j,i) + 1, nzt
+                DO  k = nzb_v_inner(j,i)+1, nzt
                    v_comp(k) = v(k,j+1,i) + v(k,j,i)
                    flux_n(k) = ( v_comp(k)- gv ) * (                          &
-                               7. * ( v(k,j+1,i) + v(k,j,i)   )               &
-                               -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
-                   diss_n(k) = - abs(v_comp(k) - gv) * (                      &
-                               3. * ( v(k,j+1,i) - v(k,j,i)   )               & 
-                               -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
+                               7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
+                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
+                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
+                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              & 
+                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
                 ENDDO
                 
              CASE ( 4 )
-                DO  k = nzb_v_inner(j,i) + 1, nzt
+                DO  k = nzb_v_inner(j,i)+1, nzt
                    v_comp(k) = v(k,j+1,i) + v(k,j,i)
                    flux_n(k) = ( v_comp(k) - gv ) *                           &
@@ -1062,16 +1069,16 @@
                 
              CASE ( 5 )
-                DO  k = nzb_w_inner(j,i) + 1, nzt
+                DO  k = nzb_w_inner(j,i)+1, nzt
                    u_comp    = u(k,j-1,i) + u(k,j,i) - gu 
                    flux_r(k) = u_comp * (                                     &
-                               7. * ( v(k,j,i+1) + v(k,j,i)   )               &
-                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
-                   diss_r(k) = - abs(u_comp) * (                              &
-                               3. * ( w(k,j,i+1) - w(k,j,i)   )               & 
-                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
+                               7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
+                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
+                   diss_r(k) = - ABS( u_comp ) * (                            &
+                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )              & 
+                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
                 ENDDO
                 
              CASE ( 6 )
-                DO  k = nzb_v_inner(j,i) + 1, nzt
+                DO  k = nzb_v_inner(j,i)+1, nzt
                    u_comp        = u(k,j-1,i) + u(k,j,i) - gu
                    flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
@@ -1082,14 +1089,14 @@
                    u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
                    flux_r(k) = u_comp * (                                     &
-                               7. * ( v(k,j,i+1) + v(k,j,i)   )               &
-                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
-                   diss_r(k) = - abs(u_comp) * (                              &
-                               3. * ( v(k,j,i+1) - v(k,j,i)   )               & 
-                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 
+                               7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
+                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
+                   diss_r(k) = - ABS( u_comp ) * (                            &
+                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              & 
+                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 
                 ENDDO
                 degraded_l = .TRUE.
                 
              CASE ( 7 )
-                DO  k = nzb_v_inner(j,i) + 1, nzt
+                DO  k = nzb_v_inner(j,i)+1, nzt
                    v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
                    flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
@@ -1100,9 +1107,9 @@
                    v_comp(k) = v(k,j+1,i) + v(k,j,i)
                    flux_n(k) = ( v_comp(k) - gv ) * (                         &
-                               7. * ( v(k,j+1,i) + v(k,j,i)   )               &
-                               -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
-                   diss_n(k) = - abs(v_comp(k) - gv) * (                      &
-                               3. * ( v(k,j+1,i) - v(k,j,i)   )               & 
-                               -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
+                               7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
+                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
+                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
+                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              & 
+                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
                 ENDDO
                 degraded_s = .TRUE.
@@ -1113,31 +1120,31 @@
 !          
 !--       Compute the crosswise 5th order fluxes at the outflow
-          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2         &
-          .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
+          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.   &
+               boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
           
-             DO  k = nzb_v_inner(j,i) + 1, nzt
+             DO  k = nzb_v_inner(j,i)+1, nzt
                 v_comp(k) = v(k,j+1,i) + v(k,j,i)
                 flux_n(k) = ( v_comp(k) - gv ) * (                            &
-                            37.  * ( v(k,j+1,i) + v(k,j,i)   )                &
-                            - 8. * ( 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(k) = - abs(v_comp(k) - gv ) * (                        &
-                             10.  * ( v(k,j+1,i) - v(k,j,i)   )               &
-                             - 5. * ( v(k,j+2,i) - v(k,j-1,i) )               &
-                             +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
+                            37.0 * ( v(k,j+1,i) + v(k,j,i)   )                &
+                          -  8.0 * ( 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(k) = - ABS( v_comp(k) - gv ) * (                       &
+                             10.0 * ( v(k,j+1,i) - v(k,j,i)   )               &
+                           -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )               &
+                           +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
               ENDDO
               
           ELSE
           
-             DO  k = nzb_v_inner(j,i) + 1, nzt
+             DO  k = nzb_v_inner(j,i)+1, nzt
                 u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                 flux_r(k) = u_comp * (                                        &
-                            37.  * ( v(k,j,i+1) + v(k,j,i)   )                &
-                            - 8. * ( 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) * (                                 &
-                            10. * ( v(k,j,i+1) - v(k,j,i)   )                 &
-                            -5. * ( v(k,j,i+2) - v(k,j,i-1) )                 &
-                            +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
+                            37.0 * ( v(k,j,i+1) + v(k,j,i)   )                &
+                          -  8.0 * ( 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 ) * (                               &
+                            10.0 * ( v(k,j,i+1) - v(k,j,i)   )                &
+                          -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                &
+                          +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
              ENDDO
              
@@ -1147,24 +1154,24 @@
 !       
 !--       Compute the fifth order fluxes for the interior of PE domain.                 
-          DO  k = nzb_v_inner(j,i) + 1, nzt
+          DO  k = nzb_v_inner(j,i)+1, nzt
              u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
              flux_r(k) = u_comp * (                                           &
-                         37.  * ( v(k,j,i+1) + v(k,j,i)   )                   &
-                         - 8. * ( 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) * (                                    &
-                         10.  * ( v(k,j,i+1) - v(k,j,i) )                     &
-                         - 5. * ( v(k,j,i+2) - v(k,j,i-1) )                   &
-                         +      ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
+                         37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
+                       -  8.0 * ( 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 ) * (                                  &
+                         10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
+                       -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
+                       +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
 
              v_comp(k) = v(k,j+1,i) + v(k,j,i)
              flux_n(k) = ( v_comp(k) - gv ) * (                               &
-                         37.  * ( v(k,j+1,i) + v(k,j,i)   )                   &
-                         - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                   &
+                         37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
+                       -  8.0 * ( 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(k) = - abs(v_comp(k) - gv) * (                            &
-                         10.  * ( v(k,j+1,i) - v(k,j,i)   )                   &
-                         - 5. * ( 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(k) = - ABS( v_comp(k) - gv ) * (                          &
+                         10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
+                       -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
+                       +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
           ENDDO 
           
@@ -1173,14 +1180,14 @@
 !--    Compute left- and southside fluxes for the respective boundary       
        IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
-          DO  k = nzb_v_inner(j,i) + 1, nzt
+          DO  k = nzb_v_inner(j,i)+1, nzt
              u_comp        = u(k,j-1,i) + u(k,j,i) - gu
              flux_l_v(k,j) = u_comp * (                                       &
-                             37.  * ( v(k,j,i) + v(k,j,i-1)   )               &
-                             - 8. * ( 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) = - abs(u_comp) * (                                &
-                             10.  * ( v(k,j,i) - v(k,j,i-1)   )               &
-                             - 5. * ( v(k,j,i+1) - v(k,j,i-2) )               &
-                             +      ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
+                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
+                           -  8.0 * ( 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) = - ABS( u_comp ) * (                              &
+                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
+                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
+                           +        ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
           ENDDO
           
@@ -1189,14 +1196,14 @@
        IF ( j == nysv .AND. ( .NOT. degraded_s ) )  THEN
        
-          DO  k = nzb_v_inner(j,i) + 1, nzt
+          DO  k = nzb_v_inner(j,i)+1, nzt
              v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
              flux_s_v(k) = v_comp_l * (                                       &
-                           37.  * ( v(k,j,i) + v(k,j-1,i)   )                 &
-                           - 8. * ( 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_v(k) = - abs(v_comp_l) * (                                &
-                           10.  * ( v(k,j,i) - v(k,j-1,i)   )                 &
-                           - 5. * ( v(k,j+1,i) - v(k,j-2,i) )                 &
-                           +      ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
+                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
+                         -  8.0 * ( 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_v(k) = - ABS( v_comp_l ) * (                              &
+                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
+                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
+                         +        ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
           ENDDO
           
@@ -1204,10 +1211,10 @@
 !       
 !--    Now compute the tendency terms for the horizontal parts.         
-       DO  k = nzb_v_inner(j,i) + 1, nzt                 
+       DO  k = nzb_v_inner(j,i)+1, nzt                 
           tend(k,j,i) = tend(k,j,i) - (                                       &
                          ( flux_r(k) + diss_r(k)                              &
-                       - ( flux_l_v(k,j) + diss_l_v(k,j) ) ) * ddx            &
+                       -   flux_l_v(k,j) - diss_l_v(k,j)   ) * ddx            &
                        + ( flux_n(k) + diss_n(k)                              &
-                       - ( flux_s_v(k) + diss_s_v(k)     ) ) * ddy     )
+                       -   flux_s_v(k) - diss_s_v(k)       ) * ddy     )
         
            flux_l_v(k,j) = flux_r(k)
@@ -1222,67 +1229,68 @@
            sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                            &
              + ( flux_n(k)                                                    &
-             *  ( v_comp(k) - 2. * hom(k,1,2,:) )                             &
+             * ( v_comp(k) - 2.0 * hom(k,1,2,:) )                             &
              / ( v_comp(k) - gv + 1.0E-20 )                                   &
-             + diss_n(k)                                                      &
-             * abs( v_comp(k) - 2. * hom(k,1,2,:) )                           &
-             / ( abs( v_comp(k) - gv ) +1.0E-20 ) )                           &
-             * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+             +   diss_n(k)                                                    &
+             *   ABS( v_comp(k) - 2.0 * hom(k,1,2,:) )                        &
+             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
+             *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
 
        ENDDO
-       sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                   &  
-                                           sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
+       sums_vs2_ws_l(nzb_v_inner(j,i),:) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)  
+                                           
 !
 !--    Vertical advection, degradation of order near surface and top. 
 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
 !--    statistical evaluation the top flux at the surface should be 0
-       flux_t(nzb_v_inner(j,i)) = 0.  !statistical reasons
-       diss_t(nzb_v_inner(j,i)) = 0.
-!
-!--    2nd order scheme       
-       k         = nzb_v_inner(j,i) + 1
+       flux_t(nzb_v_inner(j,i)) = 0.0  !statistical reasons
+       diss_t(nzb_v_inner(j,i)) = 0.0
+!
+!--    2nd order scheme (bottom)      
+       k         = nzb_v_inner(j,i)+1
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1) 
        w_comp    = w(k,j-1,i) + w(k,j,i)
        flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
-       diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0., 0., w_comp,&
-                             0.25, ddzw(k) ) 
+       diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0.0, 0.0,      &
+                             w_comp, 0.25, ddzw(k) ) 
 
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                 -   flux_d    - diss_d       ) * ddzw(k)
         
 !
-!--    WS3 as an intermediate step
-       k         = nzb_v_inner(j,i) + 2
+!--    WS3 as an intermediate step (bottom)
+       k         = nzb_v_inner(j,i)+2
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
        w_comp    = w(k,j-1,i) + w(k,j,i)
        flux_t(k) = w_comp * (                                                 &
-                   7. * ( v(k+1,j,i) + v(k,j,i)   )                           &
-                   -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
-       diss_t(k) = - abs(w_comp) * (                                          &
-                   3. * ( v(k+1,j,i) - v(k,j,i)   )                           & 
-                   -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
+                   7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &
+                 -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
+       diss_t(k) = - ABS( w_comp ) * (                                        &
+                   3.0 * ( v(k+1,j,i) - v(k,j,i)   )                          & 
+                 -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
 
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                 -   flux_d    - diss_d       ) * ddzw(k)
+!                                 
 !--    WS5
-       DO  k = nzb_v_inner(j,i) + 3, nzt - 2
+       DO  k = nzb_v_inner(j,i)+3, nzt-2
           flux_d    = flux_t(k-1)
           diss_d    = diss_t(k-1)
           w_comp    = w(k,j-1,i) + w(k,j,i)
           flux_t(k) = w_comp * (                                              &
-                      37.  * ( v(k+1,j,i) + v(k,j,i)   )                      &           
-                      - 8. * ( v(k+2,j,i) + v(k-1,j,i) )                      &
-                      +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
-          diss_t(k) = - abs(w_comp) * (                                       &
-                      10.  * ( v(k+1,j,i) - v(k,j,i)   )                      & 
-                      - 5. * ( v(k+2,j,i) - v(k-1,j,i) )                      &
-                      +      ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
+                      37.0 * ( v(k+1,j,i) + v(k,j,i)   )                      &           
+                    -  8.0 * ( v(k+2,j,i) + v(k-1,j,i) )                      &
+                    +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
+          diss_t(k) = - ABS( w_comp ) * (                                     &
+                      10.0 * ( v(k+1,j,i) - v(k,j,i)   )                      & 
+                    -  5.0 * ( v(k+2,j,i) - v(k-1,j,i) )                      &
+                    +        ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
 
           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                 -   flux_d    - diss_d       ) * ddzw(k)
        ENDDO
 !
-!--    WS3 as an intermediate step
+!--    WS3 as an intermediate step (top)
        k         = nzt - 1
        flux_d    = flux_t(k-1)
@@ -1290,13 +1298,13 @@
        w_comp    = w(k,j-1,i) + w(k,j,i)
        flux_t(k) = w_comp * (                                                 &
-                   7. * ( v(k+1,j,i) + v(k,j,i)   )                           &    
-                   -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
-       diss_t(k) = - abs(w_comp) * (                                          &
-                   3. * ( v(k+1,j,i) - v(k,j,i) )                             & 
-                   -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
+                   7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &    
+                 -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
+       diss_t(k) = - ABS( w_comp ) * (                                        &
+                   3.0 * ( v(k+1,j,i) - v(k,j,i) )                            & 
+                 -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!--    2nd order scheme 
+                                 -   flux_d    - diss_d       ) * ddzw(k)
+!
+!--    2nd order scheme (top)
        k         = nzt
        flux_d    = flux_t(k-1)
@@ -1308,10 +1316,10 @@
  
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzw(k)
+                                 -   flux_d    - diss_d       ) * ddzw(k)
               
        DO  k = nzb_v_inner(j,i), nzt
           sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:)                           &
-            + ( flux_t(k) + diss_t(k) )                                       &
-            * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+                 + ( flux_t(k) + diss_t(k) )                                  &
+                 *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
        ENDDO
 
@@ -1352,16 +1360,16 @@
          
             CASE ( 1 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                   flux_r(k) = u_comp * (                                      &
-                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
-                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
-                  diss_r(k) = -abs(u_comp) * (                                &
-                              3. * ( w(k,j,i+1) - w(k,j,i)   )                & 
-                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3  
+                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
+                            -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
+                  diss_r(k) = -ABS( u_comp ) * (                              &
+                              3.0 * ( w(k,j,i+1) - w(k,j,i)   )               & 
+                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3  
                ENDDO
                
             CASE ( 2 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                   flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 
@@ -1372,16 +1380,16 @@
                
             CASE ( 3 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                   flux_n(k) = v_comp * (                                      &
-                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
-                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
-                  diss_n(k) = -abs(v_comp) * (                                &
-                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
-                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
+                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
+                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
+                  diss_n(k) = -ABS( v_comp ) * (                              &
+                              3.0 * ( w(k,j+1,i) - w(k,j,i)   )               & 
+                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
                ENDDO
                
             CASE ( 4 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                   flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 
@@ -1392,25 +1400,25 @@
                
             CASE ( 5 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu 
                   flux_r(k) = u_comp * (                                      &
-                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
-                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
-                  diss_r(k) = - abs(u_comp) * (                               &
-                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
-                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
+                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
+                            -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
+                  diss_r(k) = - ABS( u_comp ) * (                             &
+                              3.0 * ( w(k,j,i+1) - w(k,j,i) )                 & 
+                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
                ENDDO
                
             CASE ( 6 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
 !               
 !--               Compute leftside fluxes for the left boundary of PE domain
                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                   flux_r(k) = u_comp *(                                       &
-                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
-                              -    ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
-                  diss_r(k) = - abs(u_comp) * (                               &
-                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
-                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
+                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
+                            -       ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
+                  diss_r(k) = - ABS( u_comp ) * (                             &
+                              3.0 * ( w(k,j,i+1) - w(k,j,i) )                 & 
+                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
                   
                   u_comp        = u(k+1,j,i) + u(k,j,i) - gu
@@ -1418,28 +1426,28 @@
                   diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
                                             w(k,j,i-1), w(k,j,i-1), u_comp,   &
-                                            0.25,ddx)
+                                            0.25, ddx )
                ENDDO
                degraded_l = .TRUE.
                
             CASE ( 7 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                   flux_n(k) = v_comp *(                                       &
-                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
-                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
-                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
+                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
+                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
+                  diss_n(k) = - ABS( v_comp ) * (                             &
+                              3.0 * ( w(k,j+1,i) - w(k,j,i)   )               & 
+                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
                 ENDDO
                 
             CASE ( 8 )
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                   flux_n(k) = v_comp * (                                      &
-                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
-                             -     ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              3. * ( w(k,j+1,i) - w(k,j,i) )                  & 
-                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
+                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
+                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
+                  diss_n(k) = - ABS( v_comp ) * (                             &
+                              3.0 * ( w(k,j+1,i) - w(k,j,i) )                 & 
+                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
 !                    
 !--              Compute southside fluxes for the south boundary of PE domain           
@@ -1457,31 +1465,31 @@
 !         
 !--      Compute the crosswise 5th order fluxes at the outflow
-         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
-         .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
+         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
+              boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
          
-            DO  k = nzb_w_inner(j,i) + 1, nzt
+            DO  k = nzb_w_inner(j,i)+1, nzt
                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                flux_n(k) = v_comp * (                                         &
-                           37.  * ( w(k,j+1,i) + w(k,j,i)   )                 &
-                           - 8. * ( 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) * (                                  &
-                           10.  * ( w(k,j+1,i) - w(k,j,i)   )                 &
-                           - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                 &
-                           +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
+                           37.0 * ( w(k,j+1,i) + w(k,j,i)   )                 &
+                         -  8.0 * ( 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 ) * (                                &
+                           10.0 * ( w(k,j+1,i) - w(k,j,i)   )                 &
+                         -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                 &
+                         +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
             ENDDO
             
          ELSE
          
-            DO  k = nzb_w_inner(j,i) + 1, nzt         
+            DO  k = nzb_w_inner(j,i)+1, nzt         
                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                flux_r(k) = u_comp * (                                         &
-                           37.  * ( w(k,j,i+1) + w(k,j,i) )                   &
-                           - 8. * ( 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) * (                                  &
-                           10.  * ( w(k,j,i+1) - w(k,j,i)   )                 &
-                           - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                 &
-                           +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
+                           37.0 * ( w(k,j,i+1) + w(k,j,i) )                   &
+                         -  8.0 * ( 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 ) * (                                &
+                           10.0 * ( w(k,j,i+1) - w(k,j,i)   )                 &
+                         -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                 &
+                         +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
             ENDDO
             
@@ -1491,24 +1499,24 @@
 !       
 !--      Compute the fifth order fluxes for the interior of PE domain.               
-         DO  k = nzb_w_inner(j,i) + 1, nzt
+         DO  k = nzb_w_inner(j,i)+1, nzt
             u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
             flux_r(k) = u_comp * (                                            &
-                        37.  * ( w(k,j,i+1) + w(k,j,i)   )                    &
-                        - 8. * ( 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) * (                                     &
-                        10.  * ( w(k,j,i+1) - w(k,j,i)   )                    &
-                        - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                    &
-                        +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
+                        37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
+                      -  8.0 * ( 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 ) * (                                   &
+                        10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
+                      -  5.0 * ( 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 * (                                            &
-                        37.  * ( w(k,j+1,i) + w(k,j,i)   )                    &
-                        - 8. * ( 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) * (                                     &
-                        10.  * ( w(k,j+1,i) - w(k,j,i)   )                    &
-                        - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                    &
-                        +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
+                        37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
+                      -  8.0 * ( 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 ) * (                                   &
+                        10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
+                      -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
+                      +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
          ENDDO 
          
@@ -1518,14 +1526,14 @@
        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
        
-          DO  k = nzb_w_inner(j,i) + 1, nzt
+          DO  k = nzb_w_inner(j,i)+1, nzt
              v_comp      = v(k+1,j,i) + v(k,j,i) - gv
              flux_s_w(k) = v_comp * (                                         &
-                           37.  * ( w(k,j,i) + w(k,j-1,i)   )                 &
-                           - 8. * ( 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) = - abs(v_comp) * (                                  &
-                           10.  * ( w(k,j,i) - w(k,j-1,i)   )                 &
-                           - 5. * ( w(k,j+1,i) - w(k,j-2,i) )                 &
-                           +      ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
+                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
+                         -  8.0 * ( 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) = - ABS( v_comp ) * (                                &
+                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
+                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
+                         +        ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
           ENDDO
           
@@ -1534,14 +1542,14 @@
        IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN
        
-          DO  k = nzb_w_inner(j,i) + 1, nzt
+          DO  k = nzb_w_inner(j,i)+1, nzt
             u_comp        = u(k+1,j,i) + u(k,j,i) - gu
             flux_l_w(k,j) = u_comp * (                                        &
-                            37.  * ( w(k,j,i) + w(k,j,i-1)   )                &
-                            - 8. * ( 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) = - abs(u_comp) * (                                 &
-                            10.  * ( w(k,j,i) - w(k,j,i-1)   )                &
-                            - 5. * ( w(k,j,i+1) - w(k,j,i-2) )                &
-                            +      ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 
+                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
+                          -  8.0 * ( 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) = - ABS( u_comp ) * (                               &
+                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
+                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
+                          +        ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 
           ENDDO
           
@@ -1549,10 +1557,10 @@
 !       
 !--    Now compute the tendency terms for the horizontal parts.         
-       DO  k = nzb_w_inner(j,i) + 1, nzt
+       DO  k = nzb_w_inner(j,i)+1, nzt
           tend(k,j,i) = tend(k,j,i) - (                                       &
                       ( flux_r(k) + diss_r(k)                                 &
-                    - ( flux_l_w(k,j) + diss_l_w(k,j) ) ) * ddx               &
+                    -   flux_l_w(k,j) - diss_l_w(k,j)   ) * ddx               &
                     + ( flux_n(k) + diss_n(k)                                 &
-                    - ( flux_s_w(k) + diss_s_w(k)     ) ) * ddy  )
+                    -   flux_s_w(k) - diss_s_w(k)       ) * ddy  )
 
           flux_l_w(k,j) = flux_r(k)
@@ -1566,52 +1574,52 @@
 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
 !--    statistical evaluation the top flux at the surface should be 0
-       flux_t(nzb_w_inner(j,i)) = 0.  !statistical reasons
-       diss_t(nzb_w_inner(j,i)) = 0.
-!
-!--    2nd order scheme        
-       k         = nzb_w_inner(j,i) + 1
+       flux_t(nzb_w_inner(j,i)) = 0.0  !statistical reasons
+       diss_t(nzb_w_inner(j,i)) = 0.0
+!
+!--    2nd order scheme (bottom)       
+       k         = nzb_w_inner(j,i)+1
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
        w_comp    = w(k+1,j,i) + w(k,j,i)
        flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
-       diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0.,        &
+       diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0,        &
                              w_comp, 0.25, ddzu(k+1) )
                       
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)               
-!
-!--    WS3 as an intermediate step
-       k         = nzb_w_inner(j,i) + 2
+                                 -   flux_d    - diss_d      ) * ddzu(k+1)               
+!
+!--    WS3 as an intermediate step (bottom)
+       k         = nzb_w_inner(j,i)+2
        flux_d    = flux_t(k-1)
        diss_d    = diss_t(k-1)
        w_comp    = w(k+1,j,i) + w(k,j,i)
        flux_t(k) = w_comp * (                                                 &
-                   7. * ( w(k+1,j,i) + w(k,j,i) )                             &
-                   -    ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
-       diss_t(k) = - abs(w_comp) * (                                          &
-                   3. * ( w(k+1,j,i) - w(k,j,i) )                             & 
-                   -    ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
+                   7.0 * ( w(k+1,j,i) + w(k,j,i) )                            &
+                 -       ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
+       diss_t(k) = - ABS( w_comp ) * (                                        &
+                   3.0 * ( w(k+1,j,i) - w(k,j,i) )                            & 
+                 -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
 
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
+                                 -   flux_d    - diss_d      ) * ddzu(k+1)
 !                                 
 !--    WS5
-       DO  k = nzb_w_inner(j,i) + 3, nzt -2 
+       DO  k = nzb_w_inner(j,i)+3, nzt-2 
           flux_d    = flux_t(k-1)
           diss_d    = diss_t(k-1)
           w_comp    = w(k+1,j,i) + w(k,j,i)
           flux_t(k) = w_comp * (                                              &
-                      37.  * ( w(k+1,j,i) + w(k,j,i)   )                      &
-                      - 8. * ( w(k+2,j,i) + w(k-1,j,i) )                      &
-                      +      ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
-          diss_t(k) = - abs(w_comp) * (                                       &
-                      10.  * ( w(k+1,j,i) - w(k,j,i)   )                      & 
-                      - 5. * ( w(k+2,j,i) - w(k-1,j,i) )                      & 
-                      +      ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
+                      37.0 * ( w(k+1,j,i) + w(k,j,i)   )                      &
+                    -  8.0 * ( w(k+2,j,i) + w(k-1,j,i) )                      &
+                    +        ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
+          diss_t(k) = - ABS( w_comp ) * (                                     &
+                      10.0 * ( w(k+1,j,i) - w(k,j,i)   )                      & 
+                    -  5.0 * ( w(k+2,j,i) - w(k-1,j,i) )                      & 
+                    +        ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
 
           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
+                                    -   flux_d    - diss_d      ) * ddzu(k+1)
        ENDDO
-!--    WS3 as an intermediate step
+!--    WS3 as an intermediate step (top)
        k         = nzt - 1
        flux_d    = flux_t(k-1)
@@ -1619,14 +1627,14 @@
        w_comp    = w(k+1,j,i) + w(k,j,i)
        flux_t(k) = w_comp * (                                                 &
-                   7. * ( w(k+1,j,i) + w(k,j,i)   )                           &
-                   -    ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
-       diss_t(k) = - abs(w_comp) * (                                          &
-                   3. * ( w(k+1,j,i) - w(k,j,i)     )                         & 
-                   -    ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
+                   7.0 * ( w(k+1,j,i) + w(k,j,i)   )                          &
+                 -       ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
+       diss_t(k) = - ABS( w_comp ) * (                                        &
+                   3.0 * ( w(k+1,j,i) - w(k,j,i)     )                        & 
+                 -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
                    
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
-!
-!--    2nd order scheme 
+                                 -   flux_d    - diss_d      ) * ddzu(k+1)
+!
+!--    2nd order scheme (top)
        k         = nzt
        flux_d    = flux_t(k-1)
@@ -1638,10 +1646,10 @@
 
        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
+                                 -   flux_d    - diss_d      ) * ddzu(k+1)
        
        DO  k = nzb_w_inner(j,i), nzt
            sums_ws2_ws_l(k,:)  = sums_ws2_ws_l(k,:)                           &
-               + ( flux_t(k) + diss_t(k) )                                    &
-               * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+                 + ( flux_t(k) + diss_t(k) )                                  &
+                 *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
        ENDDO
 
@@ -1649,5 +1657,5 @@
     
 
-!
+!------------------------------------------------------------------------------!
 ! Scalar advection - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -1680,5 +1688,5 @@
           IF ( boundary_flags(j,i) == 6 )  THEN
           
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 u_comp                 = u(k,j,i) - u_gtrans
                 swap_flux_x_local(k,j) = u_comp * (                           &
@@ -1692,16 +1700,16 @@
           ELSE
           
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 u_comp                 = u(k,j,i) - u_gtrans
-                swap_flux_x_local(k,j) = u_comp*(                             &
-                                         37.  * (sk(k,j,i)+sk(k,j,i-1)    )   &
-                                         -8. * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
-                                         +     ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
-                                         * adv_sca_5
-                swap_diss_x_local(k,j) = - abs(u_comp) * (                    &
-                                         10. * (sk(k,j,i) - sk(k,j,i-1)    )  &
-                                         -5. * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
-                                         +     ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
-                                         * adv_sca_5 
+                swap_flux_x_local(k,j) = u_comp*(                              &
+                                         37.0 * ( sk(k,j,i)+sk(k,j,i-1)     )  &
+                                       -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
+                                       +        ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
+                                       * adv_sca_5
+                swap_diss_x_local(k,j) = - ABS( u_comp ) * (                   &
+                                         10.0 * (sk(k,j,i) - sk(k,j,i-1)    )  &
+                                       -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
+                                       +        ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
+                                       * adv_sca_5 
              ENDDO
           ENDIF
@@ -1715,27 +1723,27 @@
           IF ( boundary_flags(j,i) == 8 )  THEN
           
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 v_comp               = v(k,j,i) - v_gtrans
-                swap_flux_y_local(k) = v_comp *                               &
+                swap_flux_y_local(k) = v_comp *                                &
                                        ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 
-                swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),    &
-                                                 sk(k,j,i), sk(k,j-1,i),      &
-                                                 sk(k,j-1,i), v_comp, 0.5, ddy)
+                swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),     &
+                                                 sk(k,j,i), sk(k,j-1,i),       &
+                                                 sk(k,j-1,i), v_comp, 0.5, ddy )
              ENDDO
              
           ELSE
           
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             DO  k = nzb_s_inner(j,i)+1, nzt
                 v_comp               = v(k,j,i) - v_gtrans
-                swap_flux_y_local(k) = v_comp * (                             &
-                                       37.  * ( sk(k,j,i) + sk(k,j-1,i)   )   &
-                                       - 8. * ( sk(k,j+1,i) + sk(k,j-2,i) )   &
-                                       +      ( sk(k,j+2,i) + sk(k,j-3,i) ) ) &  
-                                       * adv_sca_5
-                swap_diss_y_local(k)= - abs(v_comp) * (                       &
-                                      10.  * ( sk(k,j,i) - sk(k,j-1,i)   )    &
-                                      - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
-                                      +      ( sk(k,j+2,i)-sk(k,j-3,i) ) )    &
-                                      * adv_sca_5
+                swap_flux_y_local(k) = v_comp * (                              &
+                                       37.0 * ( sk(k,j,i) + sk(k,j-1,i)   )    &
+                                     -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
+                                     +        ( sk(k,j+2,i) + sk(k,j-3,i) ) )  &  
+                                     * adv_sca_5
+                swap_diss_y_local(k)= - ABS( v_comp ) * (                      &
+                                      10.0 * ( sk(k,j,i) - sk(k,j-1,i)   )     &
+                                    -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
+                                    +        ( sk(k,j+2,i)-sk(k,j-3,i) ) )     &
+                                    * adv_sca_5
              ENDDO
              
@@ -1749,85 +1757,97 @@
                
                   CASE ( 1 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         u_comp    = u(k,j,i+1) - u_gtrans
-                        flux_r(k) = u_comp * (                                &
-                                  7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
-                                  - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                                  3. * ( sk(k,j,i+1) - sk(k,j,i) )            & 
-                                  - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+                        flux_r(k) = u_comp * (                                 &
+                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
+                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
+                                  * adv_sca_3
+                        diss_r(k) = - ABS( u_comp ) * (                        &
+                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        & 
+                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
+                                  * adv_sca_3
                      ENDDO
                      
                   CASE ( 2 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         u_comp    = u(k,j,i+1) - u_gtrans
                         flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
-                        diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),       &
-                                              sk(k,j,i), sk(k,j,i-1),         &
+                        diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),        &
+                                              sk(k,j,i), sk(k,j,i-1),          &
                                               sk(k,j,i-2), u_comp, 0.5, ddx )
                      ENDDO
                      
                   CASE ( 3 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         v_comp    = v(k,j+1,i) - v_gtrans
-                        flux_n(k) = v_comp * (                                &
-                                  7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
-                                  - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                                  3. * ( sk(k,j+1,i) - sk(k,j,i) )            & 
-                                  - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
+                        flux_n(k) = v_comp * (                                 &
+                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
+                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
+                                  * adv_sca_3
+                        diss_n(k) = - ABS( v_comp ) * (                        &
+                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        & 
+                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
+                                  * adv_sca_3
                      ENDDO
                      
                   CASE ( 4 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         v_comp    = v(k,j+1,i) - v_gtrans
                         flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
-                        diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),       &
-                                              sk(k,j,i), sk(k,j-1,i),         &
+                        diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),        &
+                                              sk(k,j,i), sk(k,j-1,i),          &
                                               sk(k,j-2,i), v_comp, 0.5, ddy )     
                      ENDDO
                      
                   CASE ( 5 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         u_comp    = u(k,j,i+1) - u_gtrans 
-                        flux_r(k) = u_comp * (                                &
-                                  7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
-                                  - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                                  3. * ( sk(k,j,i+1) - sk(k,j,i) )            & 
-                                  - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+                        flux_r(k) = u_comp * (                                 &
+                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
+                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
+                                  * adv_sca_3
+                        diss_r(k) = - ABS( u_comp ) * (                        &
+                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        & 
+                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
+                                  * adv_sca_3
                      ENDDO  
                      
                   CASE ( 6 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         u_comp    = u(k,j,i+1) - u_gtrans
-                        flux_r(k) = u_comp * (                                &
-                                  7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
-                                  - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                                  3. * ( sk(k,j,i+1) - sk(k,j,i) )            & 
-                                  - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
+                        flux_r(k) = u_comp * (                                 &
+                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
+                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
+                                  * adv_sca_3
+                        diss_r(k) = - ABS( u_comp ) * (                        &
+                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        & 
+                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
+                                  * adv_sca_3
                      ENDDO
                      
                   CASE ( 7 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         v_comp    = v(k,j+1,i) - v_gtrans
-                        flux_n(k) = v_comp * (                                &
-                                  7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
-                                  - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                                  3. * ( sk(k,j+1,i) - sk(k,j,i) )            & 
-                                  - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
+                        flux_n(k) = v_comp * (                                 &
+                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
+                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
+                                  * adv_sca_3
+                        diss_n(k) = - ABS( v_comp ) * (                        &
+                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        & 
+                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
+                                  * adv_sca_3
                      ENDDO
                      
                   CASE ( 8 )
-                     DO  k = nzb_s_inner(j,i) + 1, nzt
+                     DO  k = nzb_s_inner(j,i)+1, nzt
                         v_comp    = v(k,j+1,i) - v_gtrans
-                        flux_n(k) = v_comp * (                                &
-                                  7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
-                                  - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
-                        diss_n(k) = -  abs(v_comp) * (                        &
-                                  3. * ( sk(k,j+1,i) - sk(k,j,i) )            & 
-                                  - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
+                        flux_n(k) = v_comp * (                                 &
+                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
+                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
+                                  * adv_sca_3
+                        diss_n(k) = -  ABS( v_comp ) * (                       &
+                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        & 
+                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
+                                  * adv_sca_3
                      ENDDO  
                      
@@ -1836,35 +1856,36 @@
                END SELECT
                
-               IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2     &
-             .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6)  THEN
+               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
+                   boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )    &
+               THEN
               
-                  DO  k = nzb_s_inner(j,i) + 1, nzt
+                  DO  k = nzb_s_inner(j,i)+1, nzt
                      v_comp    = v(k,j+1,i) - v_gtrans
-                     flux_n(k) = v_comp * (                                   &
-                                 37.  * ( sk(k,j+1,i) + sk(k,j,i)  )          &
-                                 - 8. * (sk(k,j+2,i) + sk(k,j-1,i) )          &
-                                 +      ( sk(k,j+3,i) + sk(k,j-2,i) ) )       &
-                                 * adv_sca_5
-                     diss_n(k) = - abs(v_comp) * (                            &
-                                 10.  * ( sk(k,j+1,i) - sk(k,j,i)   )         &
-                                 - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )         &
-                                 +      ( sk(k,j+3,i) - sk(k,j-2,i) ) )       &
-                                 * adv_sca_5
+                     flux_n(k) = v_comp * (                                    &
+                                 37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )          &
+                               -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )          &
+                               +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )        &
+                               * adv_sca_5
+                     diss_n(k) = - ABS( v_comp ) * (                           &
+                                 10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )          &
+                               -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )          &
+                               +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )        &
+                               * adv_sca_5
                   ENDDO
                   
                ELSE
                
-                  DO  k = nzb_s_inner(j,i) + 1, nzt
+                  DO  k = nzb_s_inner(j,i)+1, nzt
                      u_comp    = u(k,j,i+1) - u_gtrans  
-                     flux_r(k) = u_comp * (                                   & 
-                                 37.  * ( sk(k,j,i+1) + sk(k,j,i)   )         &
-                                 - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )         &
-                                 +      ( sk(k,j,i+3) + sk(k,j,i-2) ) )       &
-                                 * adv_sca_5
-                     diss_r(k) = - abs(u_comp) * (                            &
-                                 10.  * ( sk(k,j,i+1) - sk(k,j,i)   )         &
-                                 - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )         &
-                                 +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )       &
-                                 * adv_sca_5
+                     flux_r(k) = u_comp * (                                    & 
+                                 37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )          &
+                               -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )          &
+                               +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )        &
+                               * adv_sca_5
+                     diss_r(k) = - ABS( u_comp ) * (                           &
+                                 10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )          &
+                               -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )          &
+                               +        ( sk(k,j,i+3) - sk(k,j,i-2) ) )        &
+                               * adv_sca_5
                   ENDDO
                   
@@ -1873,38 +1894,38 @@
             ELSE
             
-               DO  k = nzb_s_inner(j,i) + 1, nzt
+               DO  k = nzb_s_inner(j,i)+1, nzt
                   u_comp    = u(k,j,i+1) - u_gtrans
-                  flux_r(k) = u_comp * (                                      &
-                              37.  * ( sk(k,j,i+1) + sk(k,j,i)   )            &
-                              - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )            &
-                              +      ( sk(k,j,i+3) + sk(k,j,i-2) ) )          &
-                              * adv_sca_5
-                  diss_r(k) = - abs(u_comp) * (                               &
-                              10.  * ( sk(k,j,i+1) - sk(k,j,i)   )            &
-                              - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )            &
-                              +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )          &
-                              * adv_sca_5
+                  flux_r(k) = u_comp * (                                       &
+                              37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
+                            -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )             &
+                            +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )           &
+                            * adv_sca_5
+                  diss_r(k) = - ABS( u_comp ) * (                              &
+                              10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
+                            -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )             &
+                            +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )             &
+                            * adv_sca_5
                   
                   v_comp    = v(k,j+1,i) - v_gtrans
-                  flux_n(k) = v_comp * (                                      &
-                              37.  * ( sk(k,j+1,i) + sk(k,j,i)   )            &
-                              - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )            &
-                              +      ( sk(k,j+3,i) + sk(k,j-2,i) ) )          &
-                              * adv_sca_5
-                  diss_n(k) = - abs(v_comp) * (                               &
-                              10.  * ( sk(k,j+1,i) - sk(k,j,i)   )            &
-                              - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )            &
-                              +      ( sk(k,j+3,i) - sk(k,j-2,i) ) )          &
-                              * adv_sca_5                  
+                  flux_n(k) = v_comp * (                                       &
+                              37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
+                            -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )             &
+                            +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )           &
+                            * adv_sca_5
+                  diss_n(k) = - ABS( v_comp ) * (                              &
+                              10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
+                            -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )             &
+                            +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )           &
+                            * adv_sca_5                  
                ENDDO
                
             ENDIF
             
-            DO  k = nzb_s_inner(j,i) + 1, nzt
+            DO  k = nzb_s_inner(j,i)+1, nzt
                tend(k,j,i) = tend(k,j,i) - (                                  &
                   ( flux_r(k) + diss_r(k)                                     &
-                - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx &
+                -   swap_flux_x_local(k,j) - swap_diss_x_local(k,j)   ) * ddx &
                 + ( flux_n(k) + diss_n(k)                                     &
-                - ( swap_flux_y_local(k) + swap_diss_y_local(k) )     ) * ddy)
+                -   swap_flux_y_local(k) - swap_diss_y_local(k)       ) * ddy)
                    
                 swap_flux_x_local(k,j) = flux_r(k)
@@ -1923,65 +1944,65 @@
          DO  j = nys, nyn
 !
-!--         2nd order scheme
-            k=nzb_s_inner(j,i) + 1
+!--         2nd order scheme (bottom)
+            k=nzb_s_inner(j,i)+1
 !            
 !--         The fluxes flux_d and diss_d at the surface are 0. Due to static
 !--         reasons the top flux at the surface should be 0.
-            flux_t(nzb_s_inner(j,i)) = 0.
-            diss_t(nzb_s_inner(j,i)) = 0.
+            flux_t(nzb_s_inner(j,i)) = 0.0
+            diss_t(nzb_s_inner(j,i)) = 0.0
             flux_d = flux_t(k-1)
             diss_d = diss_t(k-1)
-            flux_t(k) = w(k,j,i) * (sk(k+1,j,i) + sk(k,j,i) ) *0.5
+            flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
             diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i),        &
                                   sk(k,j,i), sk(k,j,i), w(k,j,i),             &
                                   0.5, ddzw(k) )
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)    
-!
-!--         WS3 as an intermediate step
-            k         = nzb_s_inner(j,i) + 2
+                                      -   flux_d    - diss_d       ) * ddzw(k)    
+!
+!--         WS3 as an intermediate step (bottom)
+            k         = nzb_s_inner(j,i)+2
             flux_d    = flux_t(k-1)
             diss_d    = diss_t(k-1)
             flux_t(k) = w(k,j,i) * (                                          &
-                        7. * ( sk(k+1,j,i) + sk(k,j,i) )                      &
+                        7.0 * ( sk(k+1,j,i) + sk(k,j,i) )                     &
                        - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 
-            diss_t(k) = - abs(w(k,j,i)) * (                                   &
-                        3. * ( sk(k+1,j,i) - sk(k,j,i) )                      & 
-                        - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
+            diss_t(k) = - ABS( w(k,j,i) ) * (                                 &
+                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   & 
+                      -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                      -   flux_d    - diss_d       ) * ddzw(k)
 !
 !--         WS5
-            DO  k = nzb_s_inner(j,i) + 3, nzt - 2
+            DO  k = nzb_s_inner(j,i)+3, nzt-2
                flux_d    = flux_t(k-1)
                diss_d    = diss_t(k-1)
                flux_t(k) = w(k,j,i) * (                                       &
-                           37.  * ( sk(k+1,j,i) + sk(k,j,i)   )               &
-                           - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
-                           +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
-               diss_t(k) = - abs(w(k,j,i)) * (                                &
-                           10.  * ( sk(k+1,j,i) -sk(k,j,i) )                  &
-                           - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
-                           +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
+                           37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )               &
+                         -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
+                         +        ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
+               diss_t(k) = - ABS(w(k,j,i)) * (                                &
+                           10.0 * ( sk(k+1,j,i) -sk(k,j,i)    )               &
+                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
+                         +        ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
 
                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                         -   flux_d    - diss_d     ) * ddzw(k)
             ENDDO
 !
-!--         WS3 as an intermediate step
+!--         WS3 as an intermediate step (top)
             k         = nzt - 1
             flux_d    = flux_t(k-1)
             diss_d    = diss_t(k-1)
             flux_t(k) = w(k,j,i) * (                                          &
-                        7. * ( sk(k+1,j,i) + sk(k,j,i)  )                     &
-                        - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
-            diss_t(k) = - abs(w(k,j,i)) * (                                   &
-                        3. * ( sk(k+1,j,i) - sk(k,j,i)   )                    & 
-                        -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
+                        7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )                   &
+                      -       ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
+            diss_t(k) = - ABS(w(k,j,i)) * (                                   &
+                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   & 
+                      -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
                         
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!--         2nd order scheme
+                                      -   flux_d    - diss_d       ) * ddzw(k)
+!
+!--         2nd order scheme (top)
             k         = nzt
             flux_d    = flux_t(k-1)
@@ -1993,5 +2014,5 @@
 
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                      -   flux_d    - diss_d       ) * ddzw(k)
 !
 !--         evaluation of statistics
@@ -2000,22 +2021,22 @@
                CASE ( 'pt' )
                  DO  k = nzb_s_inner(j,i), nzt
-                   sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                 &
-                      + ( flux_t(k) + diss_t(k) )                              &
-                      * weight_substep(intermediate_timestep_count)            &
-                      * rmask(j,i,:)
+                   sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                &
+                      + ( flux_t(k) + diss_t(k) )                             &
+                      *   weight_substep(intermediate_timestep_count)         &
+                      *   rmask(j,i,:)
                  ENDDO
                CASE ( 'sa' )
                  DO  k = nzb_s_inner(j,i), nzt
-                   sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                 &
-                      + ( flux_t(k) + diss_t(k) )                              &
-                      * weight_substep(intermediate_timestep_count)            &
-                      * rmask(j,i,:)
+                   sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                &
+                      + ( flux_t(k) + diss_t(k) )                             &
+                      *   weight_substep(intermediate_timestep_count)         &
+                      *   rmask(j,i,:)
                  ENDDO
                CASE ( 'q' )
                  DO  k = nzb_s_inner(j,i), nzt
-                   sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                   &
-                      + ( flux_t(k) + diss_t(k) )                              &
-                      * weight_substep(intermediate_timestep_count)            &
-                      * rmask(j,i,:)
+                   sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                  &
+                      + ( flux_t(k) + diss_t(k) )                             &
+                      *   weight_substep(intermediate_timestep_count)         &
+                      *   rmask(j,i,:)
                  ENDDO
 
@@ -2028,5 +2049,5 @@
 
 
-!
+!------------------------------------------------------------------------------!
 ! Advection of u - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -2059,11 +2080,11 @@
           IF( boundary_flags(j,i) == 5 )  THEN
           
-             DO  k = nzb_u_inner(j,i) + 1, nzt
+             DO  k = nzb_u_inner(j,i)+1, nzt
                 u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
-                swap_flux_x_local_u(k,j) = u_comp(k) *                        &
+                swap_flux_x_local_u(k,j) = u_comp(k) *                         &
                                            ( u(k,j,i) + u(k,j,i-1) ) * 0.25
-                swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),  &
-                                                     u(k,j,i), u(k,j,i-1),    &
-                                                     u(k,j,i-1), u_comp(k),   &
+                swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),   &
+                                                     u(k,j,i), u(k,j,i-1),     &
+                                                     u(k,j,i-1), u_comp(k),    &
                                                      0.25, ddx )
              ENDDO
@@ -2071,15 +2092,15 @@
           ELSE
           
-             DO  k = nzb_u_inner(j,i) + 1, nzt
+             DO  k = nzb_u_inner(j,i)+1, nzt
                 u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
-                swap_flux_x_local_u(k,j) = u_comp(k) * (                      &
-                                         37.  * ( u(k,j,i) + u(k,j,i-1)   )   &
-                                         - 8. * ( u(k,j,i+1) + u(k,j,i-2) )   &
-                                         +      (u(k,j,i+2)+u(k,j,i-3) )  )   &
+                swap_flux_x_local_u(k,j) = u_comp(k) * (                       &
+                                           37.0 * ( u(k,j,i) + u(k,j,i-1)   )  &
+                                         -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )  &
+                                         +        (u(k,j,i+2)+u(k,j,i-3) )  )  &
                                          * adv_mom_5
-                swap_diss_x_local_u(k,j) = - abs(u_comp(k)) * (               &
-                                         10.  * ( u(k,j,i) - u(k,j,i-1)   )   &
-                                         - 5. * ( u(k,j,i+1) - u(k,j,i-2) )   &
-                                         +      ( u(k,j,i+2) - u(k,j,i-3) ) ) &
+                swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                &
+                                           10.0 * ( u(k,j,i) - u(k,j,i-1)   )  &
+                                         -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )  &
+                                         +        ( u(k,j,i+2) - u(k,j,i-3) ) )&
                                          * adv_mom_5
              ENDDO
@@ -2096,5 +2117,5 @@
 !         
 !--         Compute southside fluxes for the south boundary of PE domain 
-            DO  k = nzb_u_inner(j,i) + 1, nzt
+            DO  k = nzb_u_inner(j,i)+1, nzt
                v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
                swap_flux_y_local_u(k) = v_comp *                              &
@@ -2108,15 +2129,15 @@
          ELSE
          
-            DO  k = nzb_u_inner(j,i) + 1, nzt
+            DO  k = nzb_u_inner(j,i)+1, nzt
                v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
-               swap_flux_y_local_u(k) = v_comp * (                            &
-                                      37.  * ( u(k,j,i) + u(k,j-1,i)   )      &
-                                      - 8. * ( u(k,j+1,i) + u(k,j-2,i) )      &
-                                      +      ( u(k,j+2,i) + u(k,j-3,i) ) )    &
+               swap_flux_y_local_u(k) = v_comp * (                             &
+                                        37.0 * ( u(k,j,i) + u(k,j-1,i)   )     &
+                                      -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )     &
+                                      +        ( u(k,j+2,i) + u(k,j-3,i) ) )   &
                                       * adv_mom_5
-               swap_diss_y_local_u(k) = - abs(v_comp) * (                     &
-                                      10.  * ( u(k,j,i) - u(k,j-1,i)   )      &
-                                      - 5. * ( u(k,j+1,i) - u(k,j-2,i) )      &
-                                      +      ( u(k,j+2,i) - u(k,j-3,i) ) )    &
+               swap_diss_y_local_u(k) = - ABS( v_comp ) * (                    &
+                                        10.0 * ( u(k,j,i) - u(k,j-1,i)   )     &
+                                      -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )     &
+                                      +        ( u(k,j+2,i) - u(k,j-3,i) ) )   &
                                       * adv_mom_5
             ENDDO
@@ -2132,16 +2153,18 @@
                
                   CASE ( 1 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt
+                     DO  k = nzb_u_inner(j,i)+1, nzt
                         u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                        flux_r(k) = ( u_comp(k) - gu ) * (                    &
-                                    7. * ( u(k,j,i+1) + u(k,j,i) )            &
-                                    - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = - abs(u_comp(k) - gu) * (                 &
-                                    3. * ( u(k,j,i+1) - u(k,j,i) )            & 
-                                    - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
+                        flux_r(k) = ( u_comp(k) - gu ) * (                     &
+                                    7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
+                                  -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
+                                    3.0 * ( u(k,j,i+1) - u(k,j,i) )            & 
+                                  -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
+                                  * adv_mom_3
                      ENDDO
                      
                   CASE ( 2 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt
+                     DO  k = nzb_u_inner(j,i)+1, nzt
                         u_comp(k) = u(k,j,i+1) + u(k,j,i)
                         flux_r(k) = ( u_comp(k) - gu ) *                      &
@@ -2153,16 +2176,18 @@
                      
                   CASE ( 3 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt
+                     DO  k = nzb_u_inner(j,i)+1, nzt
                         v_comp   = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                        flux_n(k) = v_comp * (                                &
-                                    7. * ( u(k,j+1,i) + u(k,j,i) )            &
-                                    - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                                    3. * ( u(k,j+1,i) - u(k,j,i) )            & 
-                                    - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
+                        flux_n(k) = v_comp * (                                 &
+                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
+                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS( v_comp ) * (                        &
+                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          & 
+                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
+                                  * adv_mom_3
                      ENDDO
                      
                   CASE ( 4 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt
+                     DO  k = nzb_u_inner(j,i)+1, nzt
                         v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                         flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
@@ -2173,34 +2198,40 @@
                      
                   CASE ( 5 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt        
+                     DO  k = nzb_u_inner(j,i)+1, nzt        
                         u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                        flux_r(k) = ( u_comp(k) - gu ) * (                    &
-                                    7. * ( u(k,j,i+1) + u(k,j,i) )            &
-                                    - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = - abs(u_comp(k) - gu) * (                 &
-                                    3. * ( u(k,j,i+1) - u(k,j,i) )            & 
-                                    - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
+                        flux_r(k) = ( u_comp(k) - gu ) * (                     &
+                                    7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
+                                  -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
+                                    3.0 * ( u(k,j,i+1) - u(k,j,i)   )          & 
+                                  -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
+                                  * adv_mom_3
                      ENDDO
                      
                   CASE ( 7 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt                           
+                     DO  k = nzb_u_inner(j,i)+1, nzt                           
                         v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                        flux_n(k) = v_comp * (                                &
-                                    7. * ( u(k,j+1,i) + u(k,j,i) )            &
-                                    - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                                    3. * ( u(k,j+1,i) - u(k,j,i) )            & 
-                                    - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
+                        flux_n(k) = v_comp * (                                 &
+                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
+                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS( v_comp ) * (                        &
+                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          & 
+                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
+                                  * adv_mom_3
                      ENDDO 
                      
                   CASE ( 8 )
-                     DO  k = nzb_u_inner(j,i) + 1, nzt
+                     DO  k = nzb_u_inner(j,i)+1, nzt
                         v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                        flux_n(k) = v_comp * (                                &
-                                    7. * ( u(k,j+1,i) + u(k,j,i) )            &
-                                    - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                                    3. * ( u(k,j+1,i) - u(k,j,i) )            & 
-                                    - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
+                        flux_n(k) = v_comp * (                                 &
+                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
+                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS( v_comp ) * (                        &
+                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          & 
+                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
+                                  * adv_mom_3
                      ENDDO  
                      
@@ -2210,31 +2241,35 @@
 !               
 !--            Compute the crosswise 5th order fluxes at the outflow
-               IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2     &
-              .OR. boundary_flags(j,i) == 5)  THEN
+               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
+                    boundary_flags(j,i) == 5 )  THEN
               
-                  DO  k = nzb_u_inner(j,i) + 1, nzt
+                  DO  k = nzb_u_inner(j,i)+1, nzt
                      v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                     flux_n(k) = v_comp * (                                   &
-                               37.  * ( u(k,j+1,i) + u(k,j,i)   )             &
-                               - 8. * ( 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) * (                            &
-                               10.  * ( u(k,j+1,i) - u(k,j,i) )               &
-                               - 5. * ( u(k,j+2,i) - u(k,j-1,i) )             &
-                               +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
+                     flux_n(k) = v_comp * (                                    &
+                                 37.0 * ( u(k,j+1,i) + u(k,j,i)   )            &
+                               -  8.0 * ( 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 ) * (                           &
+                                 10.0 * ( u(k,j+1,i) - u(k,j,i)   )            &
+                               -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )            &
+                               +        ( u(k,j+3,i) - u(k,j-2,i) ) )          &
+                               * adv_mom_5
                   ENDDO
                   
                ELSE
                
-                  DO  k = nzb_u_inner(j,i) + 1, nzt
+                  DO  k = nzb_u_inner(j,i)+1, nzt
                      u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                     flux_r(k) = ( u_comp(k) - gu ) * (                       &
-                               37.  * ( u(k,j,i+1) + u(k,j,i)   )             &
-                               - 8. * ( u(k,j,i+2) + u(k,j,i-1) )             &
-                               +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-                     diss_r(k) = - abs(u_comp(k) - gu) * (                    &
-                               10.  * ( u(k,j,i+1) - u(k,j,i)   )             &
-                               - 5. * ( u(k,j,i+2) - u(k,j,i-1) )             &
-                               +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
+                     flux_r(k) = ( u_comp(k) - gu ) * (                        &
+                                 37.0 * ( u(k,j,i+1) + u(k,j,i)   )            &
+                               -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )            &
+                               +        ( u(k,j,i+3) + u(k,j,i-2) ) )          &
+                               * adv_mom_5
+                     diss_r(k) = - ABS( u_comp(k) - gu ) * (                   &
+                                 10.0 * ( u(k,j,i+1) - u(k,j,i)   )            &
+                               -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )            &
+                               +        ( u(k,j,i+3) - u(k,j,i-2) ) )          &
+                               * adv_mom_5
                   ENDDO
                   
@@ -2243,24 +2278,25 @@
             ELSE
             
-               DO  k = nzb_u_inner(j,i) + 1, nzt
+               DO  k = nzb_u_inner(j,i)+1, nzt
                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
-                            37.  * ( u(k,j,i+1) + u(k,j,i)   )                &
-                            - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                &
-                            +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-                  diss_r(k) = - abs(u_comp(k) - gu) * (                       &
-                            10.  * ( u(k,j,i+1) - u(k,j,i)   )                &
-                            - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                &
-                            +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
+                  flux_r(k) = ( u_comp(k) - gu ) * (                           &
+                              37.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
+                            -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )               &
+                            +        ( u(k,j,i+3) + u(k,j,i-2) ) )             &
+                            * adv_mom_5
+                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
+                              10.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
+                            -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )               &
+                            +        ( 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 * (                                      &
-                            37.  * ( u(k,j+1,i) + u(k,j,i)   )                &
-                            - 8. * ( 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) * (                               &
-                            10.  * ( u(k,j+1,i) - u(k,j,i)   )                &
-                            - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                &
-                            +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
+                  flux_n(k) = v_comp * (                                       &
+                              37.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
+                            -  8.0 * ( 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 ) * (                              &
+                              10.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
+                            -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )               &
+                            +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
                                 
                ENDDO
@@ -2268,11 +2304,11 @@
             ENDIF
             
-            DO  k = nzb_u_inner(j,i) + 1, nzt
+            DO  k = nzb_u_inner(j,i)+1, nzt
             
-               tend(k,j,i) = tend(k,j,i) - (                                  &
-            ( flux_r(k) + diss_r(k)               &
-            - ( swap_flux_x_local_u(k,j) + swap_diss_x_local_u(k,j) ) ) * ddx &
-            + ( flux_n(k) + diss_n(k)               &
-            - ( swap_flux_y_local_u(k) + swap_diss_y_local_u(k)     ) ) * ddy )
+               tend(k,j,i) = tend(k,j,i) - (                                   &
+              ( flux_r(k) + diss_r(k)                                          &
+            -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j)   ) * ddx  &
+            + ( flux_n(k) + diss_n(k)                                          &
+            -   swap_flux_y_local_u(k) - swap_diss_y_local_u(k)       ) * ddy )
                   
                swap_flux_x_local_u(k,j) = flux_r(k)   
@@ -2281,14 +2317,14 @@
                swap_diss_y_local_u(k)   = diss_n(k)
                       
-               sums_us2_ws_l(k,:)  = sums_us2_ws_l(k,:)                       &
-                 + ( flux_r(k)                                                &
-                 * ( u_comp(k) - 2. * hom(k,1,1,:) )                          &
-                 / ( u_comp(k) - gu + 1.0E-20 )                               & 
-                 + diss_r(k)                                                  &
-                 * abs(u_comp(k) - 2. * hom(k,1,1,:) )                        &
-                 / (abs(u_comp(k) - gu) + 1.0E-20) )                          &
-                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+               sums_us2_ws_l(k,:)  = sums_us2_ws_l(k,:)                        &
+                 + ( flux_r(k)                                                 &
+                 * ( u_comp(k) - 2.0 * hom(k,1,1,:) )                          &
+                 / ( u_comp(k) - gu + 1.0E-20 )                                & 
+                 +   diss_r(k)                                                 &
+                 *   ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                     &
+                 / ( ABS( u_comp(k) - gu) + 1.0E-20) )                         &
+                 *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
             ENDDO
-            sums_us2_ws_l(nzb_u_inner(j,i),:) =                               &
+            sums_us2_ws_l(nzb_u_inner(j,i),:) =                                &
                                            sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
          ENDDO
@@ -2301,39 +2337,39 @@
        DO i = nxlu, nxr
           DO  j = nys, nyn
-             k = nzb_u_inner(j,i) + 1
+             k = nzb_u_inner(j,i)+1
 !             
 !--         The fluxes flux_d and diss_d at the surface are 0. Due to static
 !--         reasons the top flux at the surface should be 0.
-            flux_t(nzb_u_inner(j,i)) = 0.
-            diss_t(nzb_u_inner(j,i)) = 0.
+            flux_t(nzb_u_inner(j,i)) = 0.0
+            diss_t(nzb_u_inner(j,i)) = 0.0
             flux_d = flux_t(k-1)
             diss_d = diss_t(k-1)             
 !
-!--         2nd order scheme
+!--         2nd order scheme (bottom)
              w_comp    = w(k,j,i) + w(k,j,i-1)
              flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
-             diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),          &
-                                   0., 0., w_comp, 0.25, ddzw(k) )
+             diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),           &
+                                   0.0, 0.0, w_comp, 0.25, ddzw(k) )
                                    
              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!--         WS3 as an intermediate step
-            k         = nzb_u_inner(j,i) + 2
+                                       -   flux_d    -  diss_d       ) * ddzw(k)
+!
+!--         WS3 as an intermediate step (bottom)
+            k         = nzb_u_inner(j,i)+2
             flux_d    = flux_t(k-1)
             diss_d    = diss_t(k-1)
             w_comp    = w(k,j,i) + w(k,j,i-1)
             flux_t(k) = w_comp * (                                            &
-                        7. * (u(k+1,j,i) + u(k,j,i)    )                      &
-                        -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
-            diss_t(k) = - abs(w_comp) * (                                     &
-                        3. * ( u(k+1,j,i) - u(k,j,i) )                        & 
-                        -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
+                        7.0 * ( u(k+1,j,i) + u(k,j,i)   )                     &
+                      -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
+            diss_t(k) = - ABS( w_comp ) * (                                   &
+                        3.0 * ( u(k+1,j,i) - u(k,j,i)   )                     & 
+                      -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
 
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                      -   flux_d    -  diss_d      ) * ddzw(k)
 !
 !WS5
-             DO  k = nzb_u_inner(j,i) + 3, nzt - 3
+             DO  k = nzb_u_inner(j,i)+3, nzt-2
 
                 flux_d    = flux_t(k-1)
@@ -2341,34 +2377,33 @@
                 w_comp    = w(k,j,i) + w(k,j,i-1)
                 flux_t(k) = w_comp * (                                        &
-                            37.  * ( u(k+1,j,i) + u(k,j,i)   )                &
-                            - 8. * ( u(k+2,j,i) + u(k-1,j,i) )                &
-                            +      ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
-                diss_t(k) = - abs(w_comp) * (                                 &
-                            10.  * ( u(k+1,j,i) - u(k,j,i)   )                & 
-                            - 5. * ( u(k+2,j,i) - u(k-1,j,i) )                & 
-                            +      ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
+                            37.0 * ( u(k+1,j,i) + u(k,j,i)   )                &
+                          -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                &
+                          +        ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
+                diss_t(k) = - ABS( w_comp ) * (                               &
+                            10.0 * ( u(k+1,j,i) - u(k,j,i)   )                & 
+                          -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                & 
+                          +        ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
 
                 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                          -   flux_d    -  diss_d   ) * ddzw(k)
 
              ENDDO
 !
-!--         WS3 as an intermediate step
-             DO k = nzt - 2, nzt - 1
-                flux_d = flux_t(k-1)
-                diss_d = diss_t(k-1)
-                w_comp = w(k,j,i) + w(k,j,i-1)
-                flux_t(k) = w_comp * (                                        &
-                            7. * ( u(k+1,j,i) + u(k,j,i) )                    &
-                            -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
-                diss_t(k) = - abs(w_comp) * (                                 &
-                            3. * ( u(k+1,j,i) - u(k,j,i) )                    & 
-                            -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
+!--          WS3 as an intermediate step (top)
+             k         = nzt-1
+             flux_d    = flux_t(k-1)
+             diss_d    = diss_t(k-1)
+             w_comp    = w(k,j,i) + w(k,j,i-1)
+             flux_t(k) = w_comp * (                                           &
+                         7.0 * ( u(k+1,j,i) + u(k,j,i) )                      &
+                       -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
+             diss_t(k) = - ABS( w_comp ) * (                                  &
+                         3.0 * ( u(k+1,j,i) - u(k,j,i) )                      & 
+                       -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
                             
-                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
-             ENDDO
-!
-!--         2nd order scheme
+             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
+                                       -   flux_d    -  diss_d      ) * ddzw(k)
+!
+!--         2nd order scheme (top)
              k         = nzt
              flux_d    = flux_t(k-1)
@@ -2381,12 +2416,12 @@
 
              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                       -   flux_d    -  diss_d      ) * ddzw(k)
 !
 !-- at last vertical momentum flux is accumulated
             DO  k = nzb_u_inner(j,i), nzt
                sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:)                       &
-                   + ( flux_t(k) + diss_t(k) )                                 &
-                   * weight_substep(intermediate_timestep_count)               &
-                   * rmask(j,i,:)
+                              + ( flux_t(k) + diss_t(k) )                      &
+                              *   weight_substep(intermediate_timestep_count)  &
+                              *   rmask(j,i,:)
             ENDDO
           ENDDO
@@ -2397,5 +2432,5 @@
     
     
-!
+!------------------------------------------------------------------------------!
 ! Advection of v - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -2429,5 +2464,5 @@
        
           IF ( boundary_flags(j,i) == 6 )  THEN
-             DO  k = nzb_v_inner(j,i) + 1, nzt
+             DO  k = nzb_v_inner(j,i)+1, nzt
                 u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
                 swap_flux_x_local_v(k,j) = u_comp *                           &
@@ -2443,13 +2478,13 @@
              DO  k = nzb_v_inner(j,i)+1, nzt
                 u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
-                swap_flux_x_local_v(k,j) = u_comp * (                         &
-                                         37.  * ( v(k,j,i) + v(k,j,i-1)   )   &
-                                         - 8. * ( v(k,j,i+1) + v(k,j,i-2) )   &
-                                         +      ( v(k,j,i+2) + v(k,j,i-3) ) ) &
+                swap_flux_x_local_v(k,j) = u_comp * (                          &
+                                           37.0 * ( v(k,j,i) + v(k,j,i-1)   )  &
+                                         -  8.0 * ( 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. * ( v(k,j,i) - v(k,j,i-1)   )    &
-                                         -5. * ( v(k,j,i+1) - v(k,j,i-2) )    &
-                                         +     ( v(k,j,i+2) - v(k,j,i-3) ) )  &
+                swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                 &
+                                           10.0 * ( v(k,j,i) - v(k,j,i-1)   )  &
+                                         -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )  &
+                                         +        ( v(k,j,i+2) - v(k,j,i-3) ) )&
                                          * adv_mom_5
              ENDDO
@@ -2463,5 +2498,5 @@
          IF ( boundary_flags(j,i) == 7 )  THEN
          
-            DO  k = nzb_v_inner(j,i) + 1, nzt
+            DO  k = nzb_v_inner(j,i)+1, nzt
                v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
                swap_flux_y_local_v(k) = v_comp(k) *                           &
@@ -2475,15 +2510,15 @@
          ELSE
          
-            DO  k = nzb_v_inner(j,i) + 1, nzt
+            DO  k = nzb_v_inner(j,i)+1, nzt
                v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
-               swap_flux_y_local_v(k) = v_comp(k) * (                         &
-                                      37.  * ( v(k,j,i) + v(k,j-1,i)   )      &
-                                      - 8. * ( v(k,j+1,i) + v(k,j-2,i) )      &
-                                      +      ( v(k,j+2,i) + v(k,j-3,i) ) )    &
+               swap_flux_y_local_v(k) = v_comp(k) * (                          &
+                                        37.0 * ( v(k,j,i) + v(k,j-1,i)   )     &
+                                      -  8.0 * ( 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(k)) * (                  &
-                                      10.  * ( v(k,j,i) - v(k,j-1,i)   )      &
-                                      - 5. * ( v(k,j+1,i) - v(k,j-2,i) )      &
-                                      +      ( v(k,j+2,i) - v(k,j-3,i) ) )    &
+               swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                 &
+                                        10.0 * ( v(k,j,i) - v(k,j-1,i)   )     &
+                                      -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )     &
+                                      +        ( v(k,j+2,i) - v(k,j-3,i) ) )   &
                                       * adv_mom_5
             ENDDO
@@ -2498,16 +2533,18 @@
           
                   CASE ( 1 )
-                     DO  k = nzb_v_inner(j,i) + 1, nzt
+                     DO  k = nzb_v_inner(j,i)+1, nzt
                         u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
-                        flux_r(k) = u_comp * (                                &
-                                 7. * (v(k,j,i+1) + v(k,j,i)    )             &
-                                 -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                               3. * ( v(k,j,i+1) - v(k,j,i)   )               & 
-                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
+                        flux_r(k) = u_comp * (                                 &
+                                    7.0 * (v(k,j,i+1) + v(k,j,i)    )          &
+                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )        &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS( u_comp ) * (                        &
+                                    3.0 * ( v(k,j,i+1) - v(k,j,i)   )          & 
+                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )        &
+                                  * adv_mom_3
                      ENDDO
                 
                   CASE ( 2 )
-                     DO  k = nzb_v_inner(j,i) + 1, nzt
+                     DO  k = nzb_v_inner(j,i)+1, nzt
                         u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                         flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 
@@ -2518,16 +2555,18 @@
                 
                   CASE ( 3 )
-                     DO  k = nzb_v_inner(j,i) + 1, nzt
+                     DO  k = nzb_v_inner(j,i)+1, nzt
                         v_comp(k) = v(k,j+1,i) + v(k,j,i)
                         flux_n(k) = ( v_comp(k)- gv ) * (                     &
-                                 7. * ( v(k,j+1,i) + v(k,j,i)   )             &
-                                 -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp(k) - gv) * (                 &
-                                 3. * ( v(k,j+1,i) - v(k,j,i)   )             & 
-                                 -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
+                                    7.0 * ( v(k,j+1,i) + v(k,j,i)   )         &
+                                  -       ( v(k,j+2,i) + v(k,j-1,i) ) )       &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS(v_comp(k) - gv) * (                 &
+                                    3.0 * ( v(k,j+1,i) - v(k,j,i)   )         & 
+                                  -       ( v(k,j+2,i) - v(k,j-1,i) ) )       &
+                                  * adv_mom_3
                      ENDDO
                 
                   CASE ( 4 )
-                     DO  k = nzb_v_inner(j,i) + 1, nzt
+                     DO  k = nzb_v_inner(j,i)+1, nzt
                         v_comp(k) = v(k,j+1,i) + v(k,j,i)
                         flux_n(k) = ( v_comp(k) - gv ) *                      &
@@ -2539,35 +2578,41 @@
                 
                   CASE ( 5 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         u_comp    = u(k,j-1,i) + u(k,j,i) - gu 
                         flux_r(k) = u_comp * (                                &
-                               7. * ( v(k,j,i+1) + v(k,j,i)   )               &
-                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                               3. * ( w(k,j,i+1) - w(k,j,i)   )               & 
-                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
+                                    7.0 * ( v(k,j,i+1) + v(k,j,i)   )         &
+                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )       &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS (u_comp ) * (                       &
+                                    3.0 * ( w(k,j,i+1) - w(k,j,i)   )         & 
+                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )       &
+                                  * adv_mom_3
                      ENDDO
                  
                   CASE ( 6 )
-                     DO  k = nzb_v_inner(j,i) + 1, nzt
+                     DO  k = nzb_v_inner(j,i)+1, nzt
                                   
                         u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
                         flux_r(k) = u_comp * (                                &
-                                 7. * ( v(k,j,i+1) + v(k,j,i)   )             &
-                                 -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                                 3. * ( v(k,j,i+1) - v(k,j,i)   )             & 
-                                 -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 
+                                    7.0 * ( v(k,j,i+1) + v(k,j,i)   )         &
+                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )       &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS( u_comp ) * (                       &
+                                    3.0 * ( v(k,j,i+1) - v(k,j,i)   )         & 
+                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )       &
+                                  * adv_mom_3 
                      ENDDO
                 
                   CASE ( 7 )
-                     DO  k = nzb_v_inner(j,i) + 1, nzt                                 
+                     DO  k = nzb_v_inner(j,i)+1, nzt                                 
                         v_comp(k) = v(k,j+1,i) + v(k,j,i)
                         flux_n(k) = ( v_comp(k) - gv ) * (                    &
-                                 7. * ( v(k,j+1,i) + v(k,j,i)   )             &
-                                 -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp(k) - gv) * (                 &
-                                 3. * ( v(k,j+1,i) - v(k,j,i)   )             & 
-                                 -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
+                                    7.0 * ( v(k,j+1,i) + v(k,j,i)   )         &
+                                  -       ( v(k,j+2,i) + v(k,j-1,i) ) )       &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS( v_comp(k) - gv ) * (               &
+                                    3.0 * ( v(k,j+1,i) - v(k,j,i)   )         & 
+                                  -       ( v(k,j+2,i) - v(k,j-1,i) ) )       &
+                                  * adv_mom_3
                      ENDDO
                 
@@ -2576,32 +2621,37 @@
                END SELECT
 !          
-!--         Compute the crosswise 5th order fluxes at the outflow
-               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2    &
-            .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
+!--            Compute the crosswise 5th order fluxes at the outflow
+               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
+                    boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )   &
+               THEN
           
-                  DO  k = nzb_v_inner(j,i) + 1, nzt
+                  DO  k = nzb_v_inner(j,i)+1, nzt
                      v_comp(k) = v(k,j+1,i) + v(k,j,i)
                      flux_n(k) = ( v_comp(k) - gv ) * (                       &
-                               37.  * ( v(k,j+1,i) + v(k,j,i)   )             &
-                               - 8. * ( 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(k) = - abs(v_comp(k) - gv ) * (                   &
-                               10.  * ( v(k,j+1,i) - v(k,j,i)   )             &
-                               - 5. * ( v(k,j+2,i) - v(k,j-1,i) )             &
-                               +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
+                                 37.0 * ( v(k,j+1,i) + v(k,j,i)   )           &
+                               -  8.0 * ( 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(k) = - ABS( v_comp(k) - gv ) * (                  &
+                                 10.0 * ( v(k,j+1,i) - v(k,j,i)   )           &
+                               -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )           &
+                               +        ( v(k,j+3,i) - v(k,j-2,i) ) )         &
+                               * adv_mom_5
                    ENDDO
               
                ELSE
            
-                  DO  k = nzb_v_inner(j,i) + 1, nzt
+                  DO  k = nzb_v_inner(j,i)+1, nzt
                      u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                      flux_r(k) = u_comp * (                                   &
-                               37.  * ( v(k,j,i+1) + v(k,j,i)   )             &
-                              - 8. * ( 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) * (                            &
-                               10. * ( v(k,j,i+1) - v(k,j,i)   )              &
-                               -5. * ( v(k,j,i+2) - v(k,j,i-1) )              &
-                               +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
+                                 37.0 * ( v(k,j,i+1) + v(k,j,i)   )           &
+                              -   8.0 * ( 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 ) * (                          &
+                                 10.0 * ( v(k,j,i+1) - v(k,j,i)   )           &
+                               -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )           &
+                               +        ( v(k,j,i+3) - v(k,j,i-2) ) )         &
+                               * adv_mom_5
                   ENDDO
                
@@ -2611,34 +2661,35 @@
             ELSE
          
-               DO  k = nzb_v_inner(j,i) + 1, nzt
+               DO  k = nzb_v_inner(j,i)+1, nzt
                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                   flux_r(k) = u_comp * (                                      &
-                            37.  * ( v(k,j,i+1) + v(k,j,i)   )                &
-                            - 8. * ( 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) * (                               &
-                            10. * ( v(k,j,i+1) - v(k,j,i)   )                 &
-                            -5. * ( v(k,j,i+2) - v(k,j,i-1) )                 &
-                            +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
+                              37.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
+                            -  8.0 * ( 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 ) * (                            &
+                              10.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
+                            -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )              &
+                            +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
 
                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
                   flux_n(k) = ( v_comp(k) - gv ) * (                          &
-                            37.  * ( v(k,j+1,i) + v(k,j,i)   )                &
-                            - 8. * ( 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(k) = - abs(v_comp(k) - gv ) * (                      &
-                            10.  * ( v(k,j+1,i) - v(k,j,i)   )                &
-                            - 5. * ( v(k,j+2,i) - v(k,j-1,i) )                &
-                            +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
+                              37.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
+                            -  8.0 * ( 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(k) = - ABS( v_comp(k) - gv ) * (                     &
+                              10.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
+                            -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )              &
+                            +        ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
 
                ENDDO
             ENDIF
             
-            DO  k = nzb_v_inner(j,i) + 1, nzt
+            DO  k = nzb_v_inner(j,i)+1, nzt
                tend(k,j,i) = tend(k,j,i) - (                                  &
               ( flux_r(k) + diss_r(k)                                         &
-            - ( swap_flux_x_local_v(k,j) + swap_diss_x_local_v(k,j) ) ) * ddx &
+            -   swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j)   ) * ddx &
             + ( flux_n(k) + diss_n(k)                                         &
-            - ( swap_flux_y_local_v(k) + swap_diss_y_local_v(k) )     ) * ddy )
+            -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)       ) * ddy )
                   
                swap_flux_x_local_v(k,j) = flux_r(k)   
@@ -2647,12 +2698,12 @@
                swap_diss_y_local_v(k)   = diss_n(k)   
 
-               sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                        &
-                  + (flux_n(k) *( v_comp(k) - 2. * hom(k,1,2,:) )             &
-                  / ( v_comp(k) - gv + 1.0E-20)                               &
-                  + diss_n(k) * abs(v_comp(k) - 2. * hom(k,1,2,:) )           &
-                  /(abs(v_comp(k) - gv) + 1.0E-20 ) )                         &
-                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
+               sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                         &
+                  + ( flux_n(k) * ( v_comp(k) - 2.0 * hom(k,1,2,:) )           &
+                  / ( v_comp(k) - gv + 1.0E-20 )                               &
+                  +   diss_n(k) * ABS( v_comp(k) - 2.0 * hom(k,1,2,:) )        &
+                  / ( ABS( v_comp(k) - gv ) + 1.0E-20 ) )                      &
+                  *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
             ENDDO
-            sums_vs2_ws_l(nzb_v_inner(j,i),:) =                              &  
+            sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                &  
                                            sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
          ENDDO
@@ -2667,20 +2718,20 @@
 !--         The fluxes flux_d and diss_d at the surface are 0. Due to static
 !--         reasons the top flux at the surface should be 0.
-             flux_t(nzb_v_inner(j,i)) = 0.
-             diss_t(nzb_v_inner(j,i)) = 0.
-             k      = nzb_v_inner(j,i) + 1             
+             flux_t(nzb_v_inner(j,i)) = 0.0
+             diss_t(nzb_v_inner(j,i)) = 0.0
+             k      = nzb_v_inner(j,i)+1             
              flux_d = flux_t(k-1)
              diss_d = diss_t(k-1)
 !
-!--          2nd order scheme
+!--          2nd order scheme (bottom)
              w_comp    = w(k,j-1,i) + w(k,j,i)
              flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
              diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i),          &
-                                   0., 0., w_comp, 0.25, ddzw(k) )
+                                   0.0, 0.0, w_comp, 0.25, ddzw(k) )
 
              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!--         WS3 as an intermediate step
+                                       -   flux_d    - diss_d       ) * ddzw(k)
+!
+!--          WS3 as an intermediate step (bottom)
              k          = nzb_v_inner(j,i) + 2
              flux_d     = flux_t(k-1)
@@ -2688,49 +2739,48 @@
              w_comp     = w(k,j-1,i) + w(k,j,i)
              flux_t(k) = w_comp * (                                           &
-                         7. * ( v(k+1,j,i) + v(k,j,i) )                       &
-                         -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
-             diss_t(k) = - abs(w_comp) * (                                    &
-                         3. * ( v(k+1,j,i) - v(k,j,i) )                       & 
-                         -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
+                         7.0 * ( v(k+1,j,i) + v(k,j,i) )                      &
+                       -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
+             diss_t(k) = - ABS( w_comp ) * (                                  &
+                         3.0 * ( v(k+1,j,i) - v(k,j,i) )                      & 
+                       -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
                          
              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                                       -   flux_d    - diss_d       ) * ddzw(k)
 
 !
 !-- WS5
-             DO  k = nzb_v_inner(j,i) + 3, nzt - 3
+             DO  k = nzb_v_inner(j,i)+3, nzt-2
                 flux_d = flux_t(k-1)
                 diss_d = diss_t(k-1)
                 w_comp = w(k,j-1,i) + w(k,j,i)
                 flux_t(k) = w_comp * (                                        &
-                            37.  * ( v(k+1,j,i) + v(k,j,i)   )                &
-                            - 8. * ( v(k+2,j,i) + v(k-1,j,i) )                &
-                            +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
-                diss_t(k) = - abs(w_comp) * (                                 &
-                            10. * ( v(k+1,j,i) - v(k,j,i)   )                 & 
-                            -5. * ( v(k+2,j,i) - v(k-1,j,i) )                 &
-                            +     ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
+                            37.0 * ( v(k+1,j,i) + v(k,j,i)   )                &
+                          -  8.0 * ( v(k+2,j,i) + v(k-1,j,i) )                &
+                          +        ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
+                diss_t(k) = - ABS( w_comp ) * (                               &
+                            10.0 * ( v(k+1,j,i) - v(k,j,i)   )                & 
+                          -  5.0 * ( v(k+2,j,i) - v(k-1,j,i) )                &
+                          +        ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
                             
-                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
+                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
+                                          -   flux_d    - diss_d    ) * ddzw(k)
             ENDDO
 !
-!--         WS3 as an intermediate step
-            DO k = nzt - 2, nzt - 1 
-               flux_d     = flux_t(k-1)
-               diss_d     = diss_t(k-1)
-               w_comp     = w(k,j-1,i) + w(k,j,i)
-               flux_t(k) = w_comp * (                                         &
-                           7. * ( v(k+1,j,i) + v(k,j,i) )                     &
-                           -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
-               diss_t(k) = - abs(w_comp) * (                                  &
-                           3. * ( v(k+1,j,i) - v(k,j,i) )                     & 
-                           -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
+!--         WS3 as an intermediate step (top)
+            k     = nzt-1 
+            flux_d    = flux_t(k-1)
+            diss_d    = diss_t(k-1)
+            w_comp    = w(k,j-1,i) + w(k,j,i)
+            flux_t(k) = w_comp * (                                            &
+                        7.0 * ( v(k+1,j,i) + v(k,j,i)   )                     &
+                      -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
+            diss_t(k) = - ABS( w_comp ) * (                                   &
+                        3.0 * ( v(k+1,j,i) - v(k,j,i)   )                     & 
+                      -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
                          
-             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
-            ENDDO
-!
-!--          2nd order scheme
+            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
+                                      -   flux_d    - diss_d     ) * ddzw(k)
+!
+!--         2nd order scheme (top)
             k         = nzt
             flux_d    = flux_t(k-1)
@@ -2738,19 +2788,19 @@
             w_comp    = w(k,j-1,i) + w(k,j,i)
             flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
-            diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i),       &
-                                  v(k-1,j,i), v(k-2,j,i), w_comp, 0.25,ddzw(k))
-
-            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                      - ( flux_d + diss_d )     ) * ddzw(k)
-!
-!- at last vertical momentum flux is accumulated
+            diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i),        &
+                                  v(k-1,j,i), v(k-2,j,i), w_comp, 0.25, ddzw(k))
+
+            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                &
+                                      -   flux_d    - diss_d       ) * ddzw(k)
+!
+!-          At last vertical momentum flux is accumulated.
             DO  k = nzb_v_inner(j,i), nzt
                sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:)                       &
-                   + ( flux_t(k) + diss_t(k) )                                 &
-                   * weight_substep(intermediate_timestep_count)               &
-                   * rmask(j,i,:)
+                               + ( flux_t(k) + diss_t(k) )                     &
+                               *   weight_substep(intermediate_timestep_count) &
+                               *   rmask(j,i,:)
             ENDDO
             sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                &
-                   sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
+                                             sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
           ENDDO
        ENDDO
@@ -2759,5 +2809,5 @@
     
     
-!
+!------------------------------------------------------------------------------!
 ! Advection of w - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -2790,5 +2840,5 @@
           IF ( boundary_flags(j,i) == 6 )  THEN
           
-             DO  k = nzb_v_inner(j,i) + 1, nzt
+             DO  k = nzb_v_inner(j,i)+1, nzt
                 u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
                 swap_flux_x_local_w(k,j) = u_comp *                           &
@@ -2802,15 +2852,15 @@
           ELSE
           
-             DO  k = nzb_v_inner(j,i) + 1, nzt
+             DO  k = nzb_v_inner(j,i)+1, nzt
                 u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
-                swap_flux_x_local_w(k,j) = u_comp * (                         &
-                                         37.  * ( w(k,j,i) + w(k,j,i-1)   )   &
-                                         - 8. * ( w(k,j,i+1) + w(k,j,i-2) )   &
-                                         +      ( w(k,j,i+2) + w(k,j,i-3) ) ) &
+                swap_flux_x_local_w(k,j) = u_comp * (                          &
+                                           37.0 * ( w(k,j,i) + w(k,j,i-1)   )  &
+                                         -  8.0 * ( 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.  * ( w(k,j,i) - w(k,j,i-1)   )   &
-                                         - 5. * ( w(k,j,i+1) - w(k,j,i-2) )   &
-                                         +      ( w(k,j,i+2) - w(k,j,i-3) ) ) &
+                swap_diss_x_local_w(k,j) = - ABS( u_comp ) * (                 &
+                                           10.0 * ( w(k,j,i) - w(k,j,i-1)   )  &
+                                         -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )  &
+                                         +        ( w(k,j,i+2) - w(k,j,i-3) ) )&
                                          * adv_mom_5
              ENDDO
@@ -2824,26 +2874,26 @@
          IF ( boundary_flags(j,i) == 8 )  THEN
          
-            DO  k = nzb_v_inner(j,i) + 1, nzt
+            DO  k = nzb_v_inner(j,i)+1, nzt
                v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
-               swap_flux_y_local_w(k) = v_comp *                              &
+               swap_flux_y_local_w(k) = v_comp *                               &
                                         ( w(k,j,i) + w(k,j-1,i) ) * 0.25
-               swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i),     &
-                                                  w(k,j,i), w(k,j-1,i),       &
-                                                  w(k,j-1,i), v_comp, 0.25,ddy)
+               swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i),      &
+                                                  w(k,j,i), w(k,j-1,i),        &
+                                                  w(k,j-1,i), v_comp, 0.25, ddy)
             ENDDO
             
          ELSE
          
-            DO  k = nzb_v_inner(j,i) + 1, nzt            
+            DO  k = nzb_v_inner(j,i)+1, nzt            
                v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
-               swap_flux_y_local_w(k) = v_comp * (                            &
-                                      37. *  ( w(k,j,i) + w(k,j-1,i)  )       &
-                                      - 8. * ( w(k,j+1,i) +w(k,j-2,i) )       &
-                                      +      ( w(k,j+2,i) +w(k,j-3,i) ) )     &
+               swap_flux_y_local_w(k) = v_comp * (                             &
+                                        37.0 * ( w(k,j,i) + w(k,j-1,i)  )      &
+                                      -  8.0 * ( 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.  * ( w(k,j,i) - w(k,j-1,i)   )      &
-                                      - 5. * ( w(k,j+1,i) - w(k,j-2,i) )      &
-                                      +      ( w(k,j+2,i) - w(k,j-3,i) ) )    &
+               swap_diss_y_local_w(k) = - ABS( v_comp ) * (                    &
+                                        10.0 * ( w(k,j,i) - w(k,j-1,i)   )     &
+                                      -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )     &
+                                      +        ( w(k,j+2,i) - w(k,j-3,i) ) )   &
                                       * adv_mom_5
             ENDDO
@@ -2859,16 +2909,18 @@
          
                   CASE ( 1 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                         flux_r(k) = u_comp * (                                &
-                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
-                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = -abs(u_comp) * (                          &
-                              3. * ( w(k,j,i+1) - w(k,j,i)   )                & 
-                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3  
+                                    7.0 * ( w(k,j,i+1) + w(k,j,i)   )         &
+                                  -       ( w(k,j,i+2) + w(k,j,i-1) ) )       &
+                                  * adv_mom_3
+                        diss_r(k) = -ABS( u_comp ) * (                        &
+                                    3.0 * ( w(k,j,i+1) - w(k,j,i)   )         & 
+                                  -       ( w(k,j,i+2) - w(k,j,i-1) ) )       &
+                                  * adv_mom_3  
                      ENDDO
                
                   CASE ( 2 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                         flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 
@@ -2879,16 +2931,18 @@
                
                   CASE ( 3 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                         flux_n(k) = v_comp * (                                &
-                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
-                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = -abs(v_comp) * (                          &
-                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
-                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
+                                     7.0 * ( w(k,j+1,i) + w(k,j,i)   )        &
+                                  -        ( w(k,j+2,i) + w(k,j-1,i) ) )      &
+                                  * adv_mom_3
+                        diss_n(k) = -ABS( v_comp ) * (                        &
+                                    3.0 * ( w(k,j+1,i) - w(k,j,i)   )         & 
+                                  -       ( w(k,j+2,i) - w(k,j-1,i) ) )       &
+                                  * adv_mom_3
                      ENDDO
                
                   CASE ( 4 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                         flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 
@@ -2899,45 +2953,53 @@
                
                   CASE ( 5 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu 
                         flux_r(k) = u_comp * (                                &
-                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
-                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
-                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
+                                    7.0 * ( w(k,j,i+1) + w(k,j,i)   )         &
+                                  -       ( w(k,j,i+2) + w(k,j,i-1) ) )       &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS( u_comp ) * (                       &
+                                    3.0 * ( w(k,j,i+1) - w(k,j,i) )           & 
+                                  -       ( w(k,j,i+2) - w(k,j,i-1) ) )       &
+                                  * adv_mom_3
                      ENDDO
                
                   CASE ( 6 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                         flux_r(k) = u_comp *(                                 &
-                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
-                              -    ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
-                        diss_r(k) = - abs(u_comp) * (                         &
-                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
-                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
+                                    7.0 * ( w(k,j,i+1) + w(k,j,i)   )         &
+                                  -       ( w(k,j,i+2) + w(k,j,i-1) ) )       &
+                                  * adv_mom_3
+                        diss_r(k) = - ABS( u_comp ) * (                       &
+                                    3.0 * ( w(k,j,i+1) - w(k,j,i) )           & 
+                                  -       ( w(k,j,i+2) - w(k,j,i-1) ) )       &
+                                  * adv_mom_3
                      ENDDO
                
                   CASE ( 7 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                         flux_n(k) = v_comp *(                                 &
-                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
-                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
-                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
+                                    7.0 * ( w(k,j+1,i) + w(k,j,i)   )         &
+                                  -       ( w(k,j+2,i) + w(k,j-1,i) ) )       &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS( v_comp ) * (                       &
+                                    3.0 * ( w(k,j+1,i) - w(k,j,i)   )         & 
+                                  -       ( w(k,j+2,i) - w(k,j-1,i) ) )       &
+                                  * adv_mom_3
                       ENDDO
                 
                   CASE ( 8 )
-                     DO  k = nzb_w_inner(j,i) + 1, nzt
+                     DO  k = nzb_w_inner(j,i)+1, nzt
                         v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                         flux_n(k) = v_comp * (                                &
-                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
-                             -     ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
-                        diss_n(k) = - abs(v_comp) * (                         &
-                              3. * ( w(k,j+1,i) - w(k,j,i) )                  & 
-                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
+                                    7.0 * ( w(k,j+1,i) + w(k,j,i)   )         &
+                                  -       ( w(k,j+2,i) + w(k,j-1,i) ) )       &
+                                  * adv_mom_3
+                        diss_n(k) = - ABS( v_comp ) * (                       &
+                                    3.0 * ( w(k,j+1,i) - w(k,j,i) )           & 
+                                  -       ( w(k,j+2,i) - w(k,j-1,i) ) )       &
+                                  * adv_mom_3
                 
                      ENDDO  
@@ -2947,32 +3009,37 @@
                END SELECT
 !         
-!--      Compute the crosswise 5th order fluxes at the outflow
-               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2    &
-            .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
+!--            Compute the crosswise 5th order fluxes at the outflow
+               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
+                    boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )   &
+               THEN
          
-                  DO  k = nzb_w_inner(j,i) + 1, nzt
+                  DO  k = nzb_w_inner(j,i)+1, nzt
                      v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                      flux_n(k) = v_comp * (                                   &
-                           37.  * ( w(k,j+1,i) + w(k,j,i)   )                 &
-                           - 8. * ( 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) * (                            &
-                           10.  * ( w(k,j+1,i) - w(k,j,i)   )                 &
-                           - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                 &
-                           +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
+                                 37.0 * ( w(k,j+1,i) + w(k,j,i)   )           &
+                               -  8.0 * ( 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 ) * (                          &
+                                 10.0 * ( w(k,j+1,i) - w(k,j,i)   )           &
+                               -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )           &
+                               +        ( w(k,j+3,i) - w(k,j-2,i) ) )         &
+                               * adv_mom_5
                   ENDDO
             
                ELSE
          
-                  DO  k = nzb_w_inner(j,i) + 1, nzt         
+                  DO  k = nzb_w_inner(j,i)+1, nzt         
                      u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                      flux_r(k) = u_comp * (                                   &
-                           37.  * ( w(k,j,i+1) + w(k,j,i) )                   &
-                           - 8. * ( 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) * (                            &
-                           10.  * ( w(k,j,i+1) - w(k,j,i)   )                 &
-                           - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                 &
-                            +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
+                                 37.0 * ( w(k,j,i+1) + w(k,j,i) )             &
+                               -  8.0 * ( 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 ) * (                          &
+                                 10.0 * ( w(k,j,i+1) - w(k,j,i)   )           &
+                               -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )           &
+                               +        ( w(k,j,i+3) - w(k,j,i-2) ) )         &
+                               * adv_mom_5
                   ENDDO
             
@@ -2982,24 +3049,25 @@
 !       
 !--            Compute the fifth order fluxes for the interior of PE domain.               
-               DO  k = nzb_w_inner(j,i) + 1, nzt
+               DO  k = nzb_w_inner(j,i)+1, nzt
                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                   flux_r(k) = u_comp * (                                      &
-                        37.  * ( w(k,j,i+1) + w(k,j,i)   )                    &
-                        - 8. * ( 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) * (                               &
-                        10.  * ( w(k,j,i+1) - w(k,j,i)   )                    &
-                        - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                    &
-                        +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
+                              37.0 * ( w(k,j,i+1) + w(k,j,i)   )              &
+                            -  8.0 * ( 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 ) * (                             &
+                              10.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
+                            -  5.0 * ( 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 * (                                      &
-                        37.  * ( w(k,j+1,i) + w(k,j,i)   )                    &
-                        - 8. * ( 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) * (                               &
-                        10.  * ( w(k,j+1,i) - w(k,j,i)   )                    &
-                        - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                    &
-                        +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
+                              37.0 * ( w(k,j+1,i) + w(k,j,i)   )              &
+                            -  8.0 * ( 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 ) * (                             &
+                              10.0 * ( w(k,j+1,i) - w(k,j,i)   )              &
+                            -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )              &
+                            +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
                ENDDO 
          
@@ -3009,7 +3077,7 @@
                tend(k,j,i) = tend(k,j,i) - (                                  &
               ( flux_r(k) + diss_r(k)                                         &
-            - ( swap_flux_x_local_w(k,j) + swap_diss_x_local_w(k,j) ) ) * ddx &
+            -   swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j)  ) * ddx  &
             + ( flux_n(k) + diss_n(k)                                         &
-            - ( swap_flux_y_local_w(k) + swap_diss_y_local_w(k)     ) ) * ddy )
+            -   swap_flux_y_local_w(k) - swap_diss_y_local_w(k)      ) * ddy )
                
                swap_flux_x_local_w(k,j) = flux_r(k)
@@ -3027,5 +3095,5 @@
        DO i = nxl, nxr
           DO  j = nys, nyn
-            k      = nzb_w_inner(j,i) + 1
+            k      = nzb_w_inner(j,i)+1
             flux_d = w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i))
             flux_t(k-1) = flux_d
@@ -3033,64 +3101,63 @@
             diss_t(k-1) = diss_d
 !
-!--         2nd order scheme            
+!--         2nd order scheme (bottom)           
             w_comp = w(k+1,j,i) + w(k,j,i)
             flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
-            diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0.,   &
+            diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0, &
                                   w_comp, 0.25, ddzu(k+1) )
 
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
+                                      -   flux_d    - diss_d      ) * ddzu(k+1)
 !                                 
-!--        WS3 as an intermediate step
+!--        WS3 as an intermediate step (bottom)
             k         = nzb_w_inner(j,i) + 2
+            flux_d    = flux_t(k-1)
+            diss_d    = diss_t(k-1)
+            w_comp    = w(k+1,j,i) + w(k,j,i)
+            flux_t(k) = w_comp * (                                            &
+                        7.0 * ( w(k+1,j,i) + w(k,j,i) )                       &
+                      -       ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
+            diss_t(k) = - ABS( w_comp ) * (                                   &
+                        3.0 * ( w(k+1,j,i) - w(k,j,i) )                       & 
+                      -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
+                        
+            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
+                                      -   flux_d    - diss_d      ) * ddzu(k+1)
+!
+!--         WS5
+            DO  k = nzb_v_inner(j,i)+3, nzt-2
+               flux_d = flux_t(k-1)
+               diss_d = diss_t(k-1)
+
+               w_comp = w(k+1,j,i) + w(k,j,i)
+               flux_t(k) = w_comp * (                                          &
+                           37.0 * ( w(k+1,j,i) + w(k,j,i)   )                  &
+                         -  8.0 * ( w(k+2,j,i) + w(k-1,j,i) )                  &
+                         +        ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
+               diss_t(k) = - ABS( w_comp ) * (                                 &
+                           10.0 * ( w(k+1,j,i) - w(k,j,i)   )                  & 
+                         -  5.0 * ( w(k+2,j,i) - w(k-1,j,i) )                  &
+                         +        ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
+
+               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)             &
+                                         -   flux_d    - diss_d    ) * ddzu(k+1)
+            ENDDO
+!                                 
+!--         WS3 as an intermediate step (top)
+            k         = nzt-1
             flux_d    = flux_t(k-1)
             diss_d    = diss_t(k-1)
             w_comp    = w(k+1,j,i)+w(k,j,i)
             flux_t(k) = w_comp * (                                            &
-                        7. * ( w(k+1,j,i) + w(k,j,i) )                        &
-                        - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
-            diss_t(k) = - abs(w_comp) * (                                     &
-                        3. * ( w(k+1,j,i) - w(k,j,i) )                        & 
-                        - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
+                        7.0 * ( w(k+1,j,i) + w(k,j,i)   )                     &
+                      -       ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
+            diss_t(k) = - ABS( w_comp ) * (                                   &
+                          3.0 * ( w(k+1,j,i) - w(k,j,i)   )                   & 
+                      -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
                         
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
-!
-!--         WS5
-            DO  k = nzb_v_inner(j,i) + 3, nzt
-               flux_d = flux_t(k-1)
-               diss_d = diss_t(k-1)
-
-               w_comp = w(k+1,j,i) + w(k,j,i)
-               flux_t(k) = w_comp * (                                         &
-                           37.  * ( w(k+1,j,i) + w(k,j,i)   )                 &
-                           - 8. * ( w(k+2,j,i) + w(k-1,j,i) )                 &
-                           +      ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
-               diss_t(k) = - abs(w_comp) * (                                  &
-                           10.  * ( w(k+1,j,i) - w(k,j,i)   )                 & 
-                           - 5. * ( w(k+2,j,i) - w(k-1,j,i) )                 &
-                           +      ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
-
-               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
-            ENDDO
-!                                 
-!--        WS3 as an intermediate step
-            DO k = nzt - 2, nzt - 1
-               flux_d    = flux_t(k-1)
-               diss_d    = diss_t(k-1)
-               w_comp    = w(k+1,j,i)+w(k,j,i)
-               flux_t(k) = w_comp * (                                         &
-                           7. * ( w(k+1,j,i) + w(k,j,i) )                     &
-                           - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
-               diss_t(k) = - abs(w_comp) * (                                  &
-                           3. * ( w(k+1,j,i) - w(k,j,i) )                     & 
-                           - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
-                        
-               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
-            ENDDO
-!
-!--         2nd order scheme  
+                                      -   flux_d    - diss_d      ) * ddzu(k+1)
+!
+!--         2nd order scheme (top) 
             k = nzt
             flux_d = flux_t(k-1)
@@ -3103,12 +3170,12 @@
 
             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
-                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
+                                      -   flux_d    - diss_d      ) * ddzu(k+1)
 !
 !--         at last vertical momentum flux is accumulated
             DO  k = nzb_w_inner(j,i), nzt
                sums_ws2_ws_l(k,:)  = sums_ws2_ws_l(k,:)                        &
-                   + ( flux_t(k) + diss_t(k) )                                 &
-                   * weight_substep(intermediate_timestep_count)               &
-                   * rmask(j,i,:)
+                               + ( flux_t(k) + diss_t(k) )                     &
+                               *   weight_substep(intermediate_timestep_count) &
+                               *   rmask(j,i,:)
             ENDDO
 
@@ -3131,9 +3198,9 @@
 !   
 !-- Computation of 2nd order dissipation. This numerical dissipation is 
-!-- necessary to keep a stable numerical solution in region where the
-!-- order of the schemes degraded.
+!-- necessary to keep a stable numerical solution in regions where the
+!-- order of the schemes is degraded.
  
      REAL FUNCTION diss_2nd( v2, v1, v0, vm1, vm2, vel_comp, factor, grid )   &
-                            RESULT(dissip)
+                            RESULT( dissip )
         
        IMPLICIT NONE 
@@ -3144,17 +3211,17 @@
                value_min_p, value_max_p, diss_m, diss_0, diss_p
              
-       value_min_m = MIN(v0, vm1, vm2)
-       value_max_m = MAX(v0, vm1, vm2)
-       value_min = MIN(v1, v0, vm2)
-       value_max = MAX(v1, v0, vm2)
-       value_min_p = MIN(v2, v1, v0)
-       value_max_p = MAX(v2, v1, v0)
+       value_min_m = MIN( v0, vm1, vm2 )
+       value_max_m = MAX( v0, vm1, vm2 )
+       value_min   = MIN( v1, v0,  vm2 )
+       value_max   = MAX( v1, v0,  vm2 )
+       value_min_p = MIN( v2, v1,  v0  )
+       value_max_p = MAX( v2, v1,  v0  )
        
-       diss_m = MAX(0.,vm1 - value_max_m) + MIN(0.,vm1 - value_min_m)
-       diss_0 = MAX(0.,v0 - value_max) + MIN(0.,v0 - value_min)
-       diss_p = MAX(0.,v1 - value_max_p) + MIN(0.,v1 - value_min_p)
+       diss_m = MAX( 0.0, vm1 - value_max_m ) + MIN( 0.0, vm1 - value_min_m )
+       diss_0 = MAX( 0.0, v0  - value_max   ) + MIN( 0.0, v0  - value_min   )
+       diss_p = MAX( 0.0, v1  - value_max_p ) + MIN( 0.0, v1  - value_min_p )
        
-       dissip = abs(vel_comp) * (diss_p - 2.*diss_0 + diss_m)         &
-                         * grid**2 * factor
+       dissip = ABS( vel_comp ) * ( diss_p - 2.0 * diss_0 + diss_m )          &
+                                * grid**2 * factor
               
     END FUNCTION diss_2nd
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 712)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 713)
@@ -7,5 +7,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! ! Reformulate weight_substep and weight_pres as broken numbers.
 !
 ! Former revisions:
@@ -1353,19 +1353,19 @@
     IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
 
-       weight_substep(1) = 0.166666666666666
-       weight_substep(2) = 0.3
-       weight_substep(3) = 0.533333333333333
-
-       weight_pres(1)    = 0.333333333333333
-       weight_pres(2)    = 0.416666666666666
-       weight_pres(3)    = 0.25
+       weight_substep(1) = 1./6.
+       weight_substep(2) = 3./10.
+       weight_substep(3) = 8./15.
+
+       weight_pres(1)    = 1./3.
+       weight_pres(2)    = 5./12.
+       weight_pres(3)    = 1./4.
 
     ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
 
-       weight_substep(1) = 0.5
-       weight_substep(2) = 0.5
+       weight_substep(1) = 1./2.
+       weight_substep(2) = 1./2.
           
-       weight_pres(1)    = 0.5
-       weight_pres(2)    = 0.5          
+       weight_pres(1)    = 1./2.
+       weight_pres(2)    = 1./2.        
 
     ELSE                                     ! for Euler- and leapfrog-method
