Changeset 713


Ignore:
Timestamp:
Mar 30, 2011 2:21:21 PM (10 years ago)
Author:
suehring
Message:

reformatted advec_ws.f90, reformulate constants as broken numbers, bugfix in vertical advection of vertical velocity (concerning vector optimized routine)

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_ws.f90

    r711 r713  
    44! Current revisions:
    55! -----------------
    6 !
     6! File reformatted.
     7! Bugfix in vertical advection of w concerning the optimized version for   
     8! vector architecture.
     9! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
     10! broken numbers.
    711!
    812! Former revisions:
     
    1317! formatting adjustments
    1418!
    15 ! 705 2011-03-25 11:21:43 Z suehring
     19! 705 2011-03-25 11:21:43 Z suehring $
    1620! Bugfix in declaration of logicals concerning outflow boundaries.
    1721!
     
    111115!
    112116!--    Set the appropriate factors for scalar and momentum advection.
    113        adv_sca_5 = 0.016666666666666
    114        adv_sca_3 = 0.083333333333333
    115        adv_mom_5 = 0.0083333333333333
    116        adv_mom_3 = 0.041666666666666
     117       adv_sca_5 = 1./60.
     118       adv_sca_3 = 1./12.
     119       adv_mom_5 = 1./120.
     120       adv_mom_3 = 1./24.
    117121
    118122!         
     
    317321             CASE ( 3 )
    318322
    319                 DO  k = nzb_s_inner(j,i) + 1, nzt
     323                DO  k = nzb_s_inner(j,i)+1, nzt
    320324                   v_comp = v(k,j+1,i) - v_gtrans
    321325                   flux_n(k) = v_comp * (                                      &
     
    473477       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
    474478       
    475            DO  k = nzb_s_inner(j,i) + 1, nzt
     479           DO  k = nzb_s_inner(j,i)+1, nzt
    476480              v_comp               = v(k,j,i) - v_gtrans
    477481              swap_flux_y_local(k) = v_comp * (                               &
     
    534538!       
    535539!--    2nd-order scheme (bottom)
    536        k         = nzb_s_inner(j,i) + 1
     540       k         = nzb_s_inner(j,i)+1
    537541       flux_d    = flux_t(k-1)
    538542       diss_d    = diss_t(k-1)
     
    684688         
    685689            CASE ( 1 )
    686                DO  k = nzb_u_inner(j,i) + 1, nzt
     690               DO  k = nzb_u_inner(j,i)+1, nzt
    687691                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
    688692                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
    689                               7. * (u(k,j,i+1) + u(k,j,i)    )                &
    690                               -    ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    691                   diss_r(k) = - abs(u_comp(k) - gu) * (                       &
    692                               3. * (u(k,j,i+1) - u(k,j,i)    )                &
    693                               -    ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
     693                              7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
     694                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
     695                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                     &
     696                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
     697                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
    694698              ENDDO
    695699             
    696700            CASE ( 2 )
    697                DO  k = nzb_u_inner(j,i) + 1, nzt
     701               DO  k = nzb_u_inner(j,i)+1, nzt
    698702                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
    699                   flux_r(k) = (u_comp(k) - gu) * ( u(k,j,i+1) + u(k,j,i) )*0.25
     703                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
     704                                u(k,j,i+1) + u(k,j,i) ) * 0.25
    700705                  diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i),     &
    701706                                        u(k,j,i-1), u(k,j,i-2), u_comp(k),    &
     
    704709               
    705710            CASE ( 3 )
    706                DO  k = nzb_u_inner(j,i) + 1, nzt
     711               DO  k = nzb_u_inner(j,i)+1, nzt
    707712                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    708713                  flux_n(k) = v_comp * (                                      &
    709                               7. * ( u(k,j+1,i) + u(k,j,i) )                  &
    710                               -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    711                   diss_n(k) = - abs(v_comp) * (                               &
    712                               3. * ( u(k,j+1,i) - u(k,j,i)   )                &
    713                               -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
     714                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
     715                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
     716                  diss_n(k) = - ABS( v_comp ) * (                             &
     717                              3.0 * ( u(k,j+1,i) - u(k,j,i)     )             &
     718                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    714719               ENDDO
    715720               
    716721            CASE ( 4 )
    717                DO  k = nzb_u_inner(j,i) + 1, nzt
     722               DO  k = nzb_u_inner(j,i)+1, nzt
    718723                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    719724                  flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
     
    724729               
    725730            CASE ( 5 )
    726                DO  k = nzb_u_inner(j,i) + 1, nzt
     731               DO  k = nzb_u_inner(j,i)+1, nzt
    727732!               
    728733!--               Compute leftside fluxes for the left boundary of PE domain
    729734                  u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
    730735                  flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
    731                   diss_l_u(k,j) = diss_2nd(u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
    732                                            u(k,j,i-1), u(k,j,i-1), u_comp(k), &
    733                                            0.25, ddx )
     736                  diss_l_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
     737                                            u(k,j,i-1), u(k,j,i-1), u_comp(k),&
     738                                            0.25, ddx )
    734739                 
    735740                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
    736741                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
    737                               7. * (u(k,j,i+1) + u(k,j,i)    )                &
    738                               -    ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    739                   diss_r(k) = - abs( u_comp(k) -gu ) * (                      &
    740                               3. * ( u(k,j,i+1) - u(k,j,i)     )              &
    741                               -    ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
     742                              7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
     743                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
     744                  diss_r(k) = - ABS( u_comp(k) -gu ) * (                      &
     745                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
     746                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
    742747               ENDDO
    743748               degraded_l = .TRUE.
    744749               
    745750            CASE ( 7 )
    746                DO  k = nzb_u_inner(j,i) + 1, nzt                           
     751               DO  k = nzb_u_inner(j,i)+1, nzt                           
    747752                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    748753                  flux_n(k) = v_comp * (                                      &
    749                               7. * ( u(k,j+1,i) + u(k,j,i)   )                &
    750                               -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    751                   diss_n(k) = - abs(v_comp) * (                               &
    752                               3. * ( u(k,j+1,i) - u(k,j,i)   )                &
    753                               -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
     754                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
     755                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
     756                  diss_n(k) = - ABS( v_comp ) * (                             &
     757                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
     758                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    754759               ENDDO
    755760               
    756761            CASE ( 8 )
    757                DO  k = nzb_u_inner(j,i) + 1, nzt
     762               DO  k = nzb_u_inner(j,i)+1, nzt
    758763!               
    759764!--              Compute southside fluxes for the south boundary of PE domain           
     
    766771                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    767772                  flux_n(k) = v_comp * (                                      &
    768                               7. * ( u(k,j+1,i) + u(k,j,i)   )                &
    769                               -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    770                   diss_n(k) = - abs(v_comp) * (                               &
    771                               3. * ( u(k,j+1,i) - u(k,j,i)   )                &
    772                               -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
     773                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
     774                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
     775                  diss_n(k) = - ABS( v_comp ) * (                             &
     776                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
     777                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    773778               ENDDO 
    774779               degraded_s = .TRUE.
     
    779784!         
    780785!--      Compute the crosswise 5th order fluxes at the outflow
    781          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
    782          .OR. boundary_flags(j,i) == 5 )  THEN
     786         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
     787             boundary_flags(j,i) == 5 )  THEN
    783788         
    784789             DO  k = nzb_u_inner(j,i)+1, nzt
    785790               v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    786791               flux_n(k) = v_comp * (                                         &
    787                            37.  * ( u(k,j+1,i) + u(k,j,i)   )                 &
    788                            - 8. * ( u(k,j+2,i) + u(k,j-1,i) )                 &
    789                            +      ( u(k,j+3,i) + u(k,j-2,i) ) ) *adv_mom_5
    790                diss_n(k) = - abs(v_comp) * (                                  &
    791                            10.  * ( u(k,j+1,i) - u(k,j,i)   )                 &
    792                            - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                 &
    793                            +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     792                           37.0 * ( u(k,j+1,i) + u(k,j,i)   )                 &
     793                         -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                 &
     794                         +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     795               diss_n(k) = - ABS ( v_comp ) * (                               &
     796                           10.0 * ( u(k,j+1,i) - u(k,j,i)   )                 &
     797                         -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                 &
     798                               ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    794799             ENDDO
    795800             
    796801         ELSE
    797802         
    798             DO  k = nzb_u_inner(j,i) + 1, nzt
     803            DO  k = nzb_u_inner(j,i)+1, nzt
    799804               u_comp(k) = u(k,j,i+1) + u(k,j,i)
    800805               flux_r(k) = ( u_comp(k) - gu ) * (                             &
    801                            37.  * ( u(k,j,i+1) + u(k,j,i)   )                 &
    802                            - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                 &
    803                            +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    804                diss_r(k) = - abs(u_comp(k) - gu) * (                          &
    805                            10.  * ( u(k,j,i+1) - u(k,j,i) )                   &
    806                            - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                 &
    807                            +      ( u(k,j,i+3) - u(k,j,i-2) ) ) *adv_mom_5
     806                           37.0 * ( u(k,j,i+1) + u(k,j,i)   )                 &
     807                         -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                 &
     808                               ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
     809               diss_r(k) = - ABS( u_comp(k) - gu ) * (                        &
     810                           10.0 * ( u(k,j,i+1) - u(k,j,i) )                   &
     811                         -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                 &
     812                         +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    808813            ENDDO
    809814           
     
    813818!       
    814819!--       Compute the fifth order fluxes for the interior of PE domain.                 
    815           DO  k = nzb_u_inner(j,i) + 1, nzt
     820          DO  k = nzb_u_inner(j,i)+1, nzt
    816821             u_comp(k) = u(k,j,i+1) + u(k,j,i)
    817822             flux_r(k) = ( u_comp(k) - gu ) * (                               &
    818                          37.  * ( u(k,j,i+1) + u(k,j,i)   )                   &
    819                          - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                   &
    820                          +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    821              diss_r(k) = - abs(u_comp(k) - gu) * (                            &
    822                          10.  * ( u(k,j,i+1) - u(k,j,i)   )                   &
    823                          - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                   &
    824                          +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     823                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
     824                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
     825                             ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
     826             diss_r(k) = - ABS( u_comp(k) - gu ) * (                          &
     827                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
     828                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
     829                             ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    825830
    826831             v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    827832             flux_n(k) = v_comp * (                                           &
    828                          37.  * ( u(k,j+1,i) + u(k,j,i)   )                   &
    829                          - 8. * ( u(k,j+2,i) + u(k,j-1,i) )                   &
    830                          +      ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    831              diss_n(k) = - abs(v_comp) * (                                    &
    832                          10.  * ( u(k,j+1,i) - u(k,j,i)   )                   &
    833                          - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                   &
    834                          +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     833                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
     834                       -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
     835                             ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     836             diss_n(k) = - ABS( v_comp ) * (                                  &
     837                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
     838                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
     839                             ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    835840          ENDDO
    836841         
     
    840845       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
    841846       
    842           DO  k = nzb_u_inner(j,i) + 1, nzt
     847          DO  k = nzb_u_inner(j,i)+1, nzt
    843848             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    844849             flux_s_u(k) = v_comp * (                                         &
    845                            37.  * ( u(k,j,i) + u(k,j-1,i)   )                 &
    846                            - 8. * ( u(k,j+1,i) + u(k,j-2,i) )                 &
    847                            +      ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    848              diss_s_u(k) = - abs(v_comp) * (                                  &
    849                            10.  * ( u(k,j,i) - u(k,j-1,i)    )                &
    850                            - 5. * ( u(k,j+1,i) - u(k,j-2,i)  )                &
    851                            +      ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
     850                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
     851                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
     852                               ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
     853             diss_s_u(k) = - ABS(v_comp) * (                                  &
     854                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
     855                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
     856                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
    852857          ENDDO
    853858         
     
    859864             u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
    860865             flux_l_u(k,j) = u_comp_l * (                                     &
    861                              37.  * ( u(k,j,i) + u(k,j,i-1)   )               &
    862                              - 8. * ( u(k,j,i+1) + u(k,j,i-2) )               &
    863                              +      ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
    864              diss_l_u(k,j) = - abs(u_comp_l) * (                              &
    865                              10.  * ( u(k,j,i) - u(k,j,i-1)   )               &
    866                              - 5. * ( u(k,j,i+1) - u(k,j,i-2) )               &
    867                              +      ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
     866                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )               &
     867                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )               &
     868                                 ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
     869             diss_l_u(k,j) = - ABS(u_comp_l) * (                              &
     870                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )               &
     871                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )               &
     872                                 ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
    868873          ENDDO
    869874         
     
    871876!       
    872877!--    Now compute the tendency terms for the horizontal parts.
    873        DO  k = nzb_u_inner(j,i) + 1, nzt                   
     878       DO  k = nzb_u_inner(j,i)+1, nzt                   
    874879          tend(k,j,i) = tend(k,j,i) - (                                       &
    875880                            ( flux_r(k) + diss_r(k)                           &
    876                           - ( flux_l_u(k,j) + diss_l_u(k,j) ) ) * ddx         &
     881                          -   flux_l_u(k,j) - diss_l_u(k,j) ) * ddx         &
    877882                          + ( flux_n(k) + diss_n(k)                           &
    878                           - ( flux_s_u(k) + diss_s_u(k) )     ) * ddy  )
     883                          -   flux_s_u(k) - diss_s_u(k)     ) * ddy  )
    879884
    880885           flux_l_u(k,j) = flux_r(k)
     
    887892           sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:)                            &
    888893             + ( flux_r(k) *                                                  &
    889              ( u_comp(k) - 2. * hom(k,1,1,:) ) / ( u_comp(k) - gu + 1.0E-20 ) &
    890              + diss_r(k)                                                      &
    891              * abs(u_comp(k) - 2. * hom(k,1,1,:) )                            &
    892              / (abs(u_comp(k) - gu) + 1.0E-20)         )                      &
    893              * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     894               ( u_comp(k) - 2.0 * hom(k,1,1,:) )                             &
     895             / ( u_comp(k) - gu + 1.0E-20      )                              &
     896             +   diss_r(k) *                                                  &
     897                 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                        &
     898             / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )                          &
     899             *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    894900       ENDDO
    895        sums_us2_ws_l(nzb_u_inner(j,i),:) =                                   &
    896                                            sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
     901       sums_us2_ws_l(nzb_u_inner(j,i),:) = sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
     902                                           
    897903
    898904!
     
    903909       diss_t(nzb_u_inner(j,i)) = 0.
    904910!
    905 !--    2nd order scheme
    906        k         = nzb_u_inner(j,i) + 1
     911!--    2nd order scheme (bottom)
     912       k         = nzb_u_inner(j,i)+1
    907913       flux_d    = flux_t(k-1)
    908914       diss_d    = diss_t(k-1)
     
    913919             
    914920       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    915                                  - ( flux_d + diss_d )     ) * ddzw(k)
    916 !
    917 !--    WS3 as an intermediate step
    918        k         = nzb_u_inner(j,i) + 2
     921                                 -   flux_d    - diss_d      ) * ddzw(k)
     922!
     923!--    WS3 as an intermediate step (bottom)
     924       k         = nzb_u_inner(j,i)+2
    919925       flux_d    = flux_t(k-1)
    920926       diss_d    = diss_t(k-1)
    921927       w_comp    = w(k,j,i) + w(k,j,i-1)
    922        flux_t(k) = w_comp*(                                                   &
    923                    7. * ( u(k+1,j,i) + u(k,j,i)   )                           &
    924                    -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    925        diss_t(k) = -abs(w_comp)*(                                             &
    926                    3. * ( u(k+1,j,i) - u(k,j,i)   )                           &
    927                    -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
     928       flux_t(k) = w_comp * (                                                 &
     929                   7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
     930                 -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
     931       diss_t(k) = -ABS( w_comp) * (                                          &
     932                   3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          &
     933                 -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    928934
    929935       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    930                                  - ( flux_d + diss_d )     ) * ddzw(k)
    931 
    932        DO  k = nzb_u_inner(j,i) + 3, nzt - 2
     936                                 -   flux_d    - diss_d      ) * ddzw(k)
     937!
     938!--    WS5
     939       DO  k = nzb_u_inner(j,i)+3, nzt-2
    933940          flux_d    = flux_t(k-1)
    934941          diss_d    = diss_t(k-1)
    935942          w_comp    = w(k,j,i) + w(k,j,i-1)
    936           flux_t(k) = w_comp*(                                                &
    937                       37. ( u(k+1,j,i) + u(k,j,i)   )                      &
    938                       - 8. * ( u(k+2,j,i) + u(k-1,j,i) )                      &
    939                       +      ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
    940           diss_t(k) = - abs(w_comp) * (                                       &
    941                       10. *  ( u(k+1,j,i) - u(k,j,i)   )                      &
    942                       - 5. * ( u(k+2,j,i) - u(k-1,j,i) )                      &
    943                       +      ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
     943          flux_t(k) = w_comp * (                                              &
     944                      37.0 * ( u(k+1,j,i) + u(k,j,i)   )                      &
     945                    -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                      &
     946                          ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
     947          diss_t(k) = - ABS( w_comp ) * (                                     &
     948                      10.0 *  ( u(k+1,j,i) - u(k,j,i)   )                     &
     949                    -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                      &
     950                          ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
    944951
    945952          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    946                                  - ( flux_d + diss_d )     ) * ddzw(k)
     953                                    -   flux_d    - diss_d      ) * ddzw(k)
    947954       ENDDO
    948955!
    949 !--    WS3 as an intermediate step
     956!--    WS3 as an intermediate step (top)
    950957       k         = nzt - 1
    951958       flux_d    = flux_t(k-1)
     
    953960       w_comp    = w(k,j,i) + w(k,j,i-1)
    954961       flux_t(k) = w_comp * (                                                 &
    955                    7. * ( u(k+1,j,i) + u(k,j,i)   )                           &
    956                    -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    957        diss_t(k) = - abs(w_comp) * (                                          &
    958                    3. * ( u(k+1,j,i) - u(k,j,i)   )                           &
    959                    -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
     962                   7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
     963                 -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
     964       diss_t(k) = - ABS( w_comp ) * (                                        &
     965                   3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          &
     966                 -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    960967                   
    961968       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    962                                  - ( flux_d + diss_d )     ) * ddzw(k)
     969                                 -   flux_d    - diss_d      ) * ddzw(k)
    963970       
    964971!
    965 !--    2nd order scheme
     972!--    2nd order scheme (top)
    966973       k         = nzt
    967974       flux_d    = flux_t(k-1)
     
    973980
    974981       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    975                                  - ( flux_d + diss_d )     ) * ddzw(k)
     982                                 -   flux_d    - diss_d      ) * ddzw(k)
    976983
    977984!
     
    980987          sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:)                           &
    981988              + ( flux_t(k) + diss_t(k) )                                     &
    982               * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     989              *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    983990       ENDDO
    984991
     
    10211028         
    10221029             CASE ( 1 )
    1023                 DO  k = nzb_v_inner(j,i) + 1, nzt
     1030                DO  k = nzb_v_inner(j,i)+1, nzt
    10241031                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    10251032                   flux_r(k) = u_comp * (                                     &
    1026                                7. * (v(k,j,i+1) + v(k,j,i)    )               &
    1027                                -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    1028                    diss_r(k) = - abs(u_comp) * (                              &
    1029                                3. * ( v(k,j,i+1) - v(k,j,i)   )               &
    1030                                -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
     1033                               7.0 * (v(k,j,i+1) + v(k,j,i)    )              &
     1034                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
     1035                   diss_r(k) = - ABS( u_comp ) * (                            &
     1036                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
     1037                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    10311038                ENDDO
    10321039               
    10331040             CASE ( 2 )
    1034                 DO  k = nzb_v_inner(j,i) + 1, nzt
     1041                DO  k = nzb_v_inner(j,i)+1, nzt
    10351042                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    10361043                   flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25
     
    10411048               
    10421049             CASE ( 3 )
    1043                 DO  k = nzb_v_inner(j,i) + 1, nzt
     1050                DO  k = nzb_v_inner(j,i)+1, nzt
    10441051                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
    10451052                   flux_n(k) = ( v_comp(k)- gv ) * (                          &
    1046                                7. * ( v(k,j+1,i) + v(k,j,i)   )               &
    1047                                -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    1048                    diss_n(k) = - abs(v_comp(k) - gv) * (                      &
    1049                                3. * ( v(k,j+1,i) - v(k,j,i)   )               &
    1050                                -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
     1053                               7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
     1054                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
     1055                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
     1056                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
     1057                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
    10511058                ENDDO
    10521059               
    10531060             CASE ( 4 )
    1054                 DO  k = nzb_v_inner(j,i) + 1, nzt
     1061                DO  k = nzb_v_inner(j,i)+1, nzt
    10551062                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
    10561063                   flux_n(k) = ( v_comp(k) - gv ) *                           &
     
    10621069               
    10631070             CASE ( 5 )
    1064                 DO  k = nzb_w_inner(j,i) + 1, nzt
     1071                DO  k = nzb_w_inner(j,i)+1, nzt
    10651072                   u_comp    = u(k,j-1,i) + u(k,j,i) - gu
    10661073                   flux_r(k) = u_comp * (                                     &
    1067                                7. * ( v(k,j,i+1) + v(k,j,i)   )               &
    1068                                -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    1069                    diss_r(k) = - abs(u_comp) * (                              &
    1070                                3. * ( w(k,j,i+1) - w(k,j,i)   )               &
    1071                                -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
     1074                               7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
     1075                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
     1076                   diss_r(k) = - ABS( u_comp ) * (                            &
     1077                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
     1078                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    10721079                ENDDO
    10731080               
    10741081             CASE ( 6 )
    1075                 DO  k = nzb_v_inner(j,i) + 1, nzt
     1082                DO  k = nzb_v_inner(j,i)+1, nzt
    10761083                   u_comp        = u(k,j-1,i) + u(k,j,i) - gu
    10771084                   flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
     
    10821089                   u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
    10831090                   flux_r(k) = u_comp * (                                     &
    1084                                7. * ( v(k,j,i+1) + v(k,j,i)   )               &
    1085                                -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    1086                    diss_r(k) = - abs(u_comp) * (                              &
    1087                                3. * ( v(k,j,i+1) - v(k,j,i)   )               &
    1088                                -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
     1091                               7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
     1092                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
     1093                   diss_r(k) = - ABS( u_comp ) * (                            &
     1094                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
     1095                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    10891096                ENDDO
    10901097                degraded_l = .TRUE.
    10911098               
    10921099             CASE ( 7 )
    1093                 DO  k = nzb_v_inner(j,i) + 1, nzt
     1100                DO  k = nzb_v_inner(j,i)+1, nzt
    10941101                   v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
    10951102                   flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
     
    11001107                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
    11011108                   flux_n(k) = ( v_comp(k) - gv ) * (                         &
    1102                                7. * ( v(k,j+1,i) + v(k,j,i)   )               &
    1103                                -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    1104                    diss_n(k) = - abs(v_comp(k) - gv) * (                      &
    1105                                3. * ( v(k,j+1,i) - v(k,j,i)   )               &
    1106                                -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
     1109                               7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
     1110                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
     1111                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
     1112                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
     1113                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
    11071114                ENDDO
    11081115                degraded_s = .TRUE.
     
    11131120!         
    11141121!--       Compute the crosswise 5th order fluxes at the outflow
    1115           IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2         &
    1116           .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
     1122          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.   &
     1123              boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
    11171124         
    1118              DO  k = nzb_v_inner(j,i) + 1, nzt
     1125             DO  k = nzb_v_inner(j,i)+1, nzt
    11191126                v_comp(k) = v(k,j+1,i) + v(k,j,i)
    11201127                flux_n(k) = ( v_comp(k) - gv ) * (                            &
    1121                             37.  * ( v(k,j+1,i) + v(k,j,i)   )                &
    1122                             - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                &
    1123                             +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    1124                 diss_n(k) = - abs(v_comp(k) - gv ) * (                        &
    1125                              10.  * ( v(k,j+1,i) - v(k,j,i)   )               &
    1126                              - 5. * ( v(k,j+2,i) - v(k,j-1,i) )               &
    1127                              +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
     1128                            37.0 * ( v(k,j+1,i) + v(k,j,i)   )                &
     1129                          -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                &
     1130                                ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
     1131                diss_n(k) = - ABS( v_comp(k) - gv ) * (                       &
     1132                             10.0 * ( v(k,j+1,i) - v(k,j,i)   )               &
     1133                           -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )               &
     1134                           +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    11281135              ENDDO
    11291136             
    11301137          ELSE
    11311138         
    1132              DO  k = nzb_v_inner(j,i) + 1, nzt
     1139             DO  k = nzb_v_inner(j,i)+1, nzt
    11331140                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    11341141                flux_r(k) = u_comp * (                                        &
    1135                             37.  * ( v(k,j,i+1) + v(k,j,i)   )                &
    1136                             - 8. * ( v(k,j,i+2) + v(k,j,i-1) )                &
    1137                             +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    1138                 diss_r(k) = - abs(u_comp) * (                                 &
    1139                             10. * ( v(k,j,i+1) - v(k,j,i)   )                 &
    1140                             -5. * ( v(k,j,i+2) - v(k,j,i-1) )                 &
    1141                             +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     1142                            37.0 * ( v(k,j,i+1) + v(k,j,i)   )                &
     1143                          -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                &
     1144                                ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
     1145                diss_r(k) = - ABS( u_comp ) * (                               &
     1146                            10.0 * ( v(k,j,i+1) - v(k,j,i)   )                &
     1147                          -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                &
     1148                          +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    11421149             ENDDO
    11431150             
     
    11471154!       
    11481155!--       Compute the fifth order fluxes for the interior of PE domain.                 
    1149           DO  k = nzb_v_inner(j,i) + 1, nzt
     1156          DO  k = nzb_v_inner(j,i)+1, nzt
    11501157             u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    11511158             flux_r(k) = u_comp * (                                           &
    1152                          37.  * ( v(k,j,i+1) + v(k,j,i)   )                   &
    1153                          - 8. * ( v(k,j,i+2) + v(k,j,i-1) )                   &
    1154                          +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    1155              diss_r(k) = - abs(u_comp) * (                                    &
    1156                          10.  * ( v(k,j,i+1) - v(k,j,i) )                     &
    1157                          - 5. * ( v(k,j,i+2) - v(k,j,i-1) )                   &
    1158                          +      ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     1159                         37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
     1160                       -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                   &
     1161                             ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
     1162             diss_r(k) = - ABS( u_comp ) * (                                  &
     1163                         10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
     1164                       -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
     1165                             ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    11591166
    11601167             v_comp(k) = v(k,j+1,i) + v(k,j,i)
    11611168             flux_n(k) = ( v_comp(k) - gv ) * (                               &
    1162                          37.  * ( v(k,j+1,i) + v(k,j,i)   )                   &
    1163                          - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                   &
     1169                         37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
     1170                       -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                   &
    11641171                         +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    1165              diss_n(k) = - abs(v_comp(k) - gv) * (                            &
    1166                          10.  * ( v(k,j+1,i) - v(k,j,i)   )                   &
    1167                          - 5. * ( v(k,j+2,i) - v(k,j-1,i) )                   &
    1168                          +      ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
     1172             diss_n(k) = - ABS( v_comp(k) - gv ) * (                          &
     1173                         10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
     1174                       -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
     1175                             ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    11691176          ENDDO
    11701177         
     
    11731180!--    Compute left- and southside fluxes for the respective boundary       
    11741181       IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
    1175           DO  k = nzb_v_inner(j,i) + 1, nzt
     1182          DO  k = nzb_v_inner(j,i)+1, nzt
    11761183             u_comp        = u(k,j-1,i) + u(k,j,i) - gu
    11771184             flux_l_v(k,j) = u_comp * (                                       &
    1178                              37.  * ( v(k,j,i) + v(k,j,i-1)   )               &
    1179                              - 8. * ( v(k,j,i+1) + v(k,j,i-2) )               &
    1180                              +      ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    1181              diss_l_v(k,j) = - abs(u_comp) * (                                &
    1182                              10.  * ( v(k,j,i) - v(k,j,i-1)   )               &
    1183                              - 5. * ( v(k,j,i+1) - v(k,j,i-2) )               &
    1184                              +      ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
     1185                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
     1186                           -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
     1187                                 ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
     1188             diss_l_v(k,j) = - ABS( u_comp ) * (                              &
     1189                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
     1190                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
     1191                                 ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
    11851192          ENDDO
    11861193         
     
    11891196       IF ( j == nysv .AND. ( .NOT. degraded_s ) )  THEN
    11901197       
    1191           DO  k = nzb_v_inner(j,i) + 1, nzt
     1198          DO  k = nzb_v_inner(j,i)+1, nzt
    11921199             v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
    11931200             flux_s_v(k) = v_comp_l * (                                       &
    1194                            37.  * ( v(k,j,i) + v(k,j-1,i)   )                 &
    1195                            - 8. * ( v(k,j+1,i) + v(k,j-2,i) )                 &
    1196                            +      ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
    1197              diss_s_v(k) = - abs(v_comp_l) * (                                &
    1198                            10.  * ( v(k,j,i) - v(k,j-1,i)   )                 &
    1199                            - 5. * ( v(k,j+1,i) - v(k,j-2,i) )                 &
    1200                            +      ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
     1201                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
     1202                         -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
     1203                               ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
     1204             diss_s_v(k) = - ABS( v_comp_l ) * (                              &
     1205                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
     1206                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
     1207                               ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
    12011208          ENDDO
    12021209         
     
    12041211!       
    12051212!--    Now compute the tendency terms for the horizontal parts.         
    1206        DO  k = nzb_v_inner(j,i) + 1, nzt                 
     1213       DO  k = nzb_v_inner(j,i)+1, nzt                 
    12071214          tend(k,j,i) = tend(k,j,i) - (                                       &
    12081215                         ( flux_r(k) + diss_r(k)                              &
    1209                        - ( flux_l_v(k,j) + diss_l_v(k,j) ) ) * ddx            &
     1216                       -   flux_l_v(k,j) - diss_l_v(k,j)  ) * ddx            &
    12101217                       + ( flux_n(k) + diss_n(k)                              &
    1211                        - ( flux_s_v(k) + diss_s_v(k)     ) ) * ddy     )
     1218                       -   flux_s_v(k) - diss_s_v(k)      ) * ddy     )
    12121219       
    12131220           flux_l_v(k,j) = flux_r(k)
     
    12221229           sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                            &
    12231230             + ( flux_n(k)                                                    &
    1224              *  ( v_comp(k) - 2. * hom(k,1,2,:) )                             &
     1231             * ( v_comp(k) - 2.0 * hom(k,1,2,:) )                             &
    12251232             / ( v_comp(k) - gv + 1.0E-20 )                                   &
    1226              + diss_n(k)                                                      &
    1227              * abs( v_comp(k) - 2. * hom(k,1,2,:) )                           &
    1228              / ( abs( v_comp(k) - gv ) +1.0E-20 ) )                           &
    1229              * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     1233             +   diss_n(k)                                                    &
     1234             *   ABS( v_comp(k) - 2.0 * hom(k,1,2,:) )                        &
     1235             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
     1236             *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    12301237
    12311238       ENDDO
    1232        sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                   & 
    1233                                            sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
     1239       sums_vs2_ws_l(nzb_v_inner(j,i),:) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) 
     1240                                           
    12341241!
    12351242!--    Vertical advection, degradation of order near surface and top.
    12361243!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    12371244!--    statistical evaluation the top flux at the surface should be 0
    1238        flux_t(nzb_v_inner(j,i)) = 0.  !statistical reasons
    1239        diss_t(nzb_v_inner(j,i)) = 0.
    1240 !
    1241 !--    2nd order scheme      
    1242        k         = nzb_v_inner(j,i) + 1
     1245       flux_t(nzb_v_inner(j,i)) = 0.0  !statistical reasons
     1246       diss_t(nzb_v_inner(j,i)) = 0.0
     1247!
     1248!--    2nd order scheme (bottom)     
     1249       k         = nzb_v_inner(j,i)+1
    12431250       flux_d    = flux_t(k-1)
    12441251       diss_d    = diss_t(k-1)
    12451252       w_comp    = w(k,j-1,i) + w(k,j,i)
    12461253       flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
    1247        diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0., 0., w_comp,&
    1248                              0.25, ddzw(k) )
     1254       diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0.0, 0.0,      &
     1255                             w_comp, 0.25, ddzw(k) )
    12491256
    12501257       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1251                                  - ( flux_d + diss_d )     ) * ddzw(k)
     1258                                 -   flux_d    - diss_d       ) * ddzw(k)
    12521259       
    12531260!
    1254 !--    WS3 as an intermediate step
    1255        k         = nzb_v_inner(j,i) + 2
     1261!--    WS3 as an intermediate step (bottom)
     1262       k         = nzb_v_inner(j,i)+2
    12561263       flux_d    = flux_t(k-1)
    12571264       diss_d    = diss_t(k-1)
    12581265       w_comp    = w(k,j-1,i) + w(k,j,i)
    12591266       flux_t(k) = w_comp * (                                                 &
    1260                    7. * ( v(k+1,j,i) + v(k,j,i)   )                           &
    1261                    -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
    1262        diss_t(k) = - abs(w_comp) * (                                          &
    1263                    3. * ( v(k+1,j,i) - v(k,j,i)   )                           &
    1264                    -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
     1267                   7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &
     1268                 -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
     1269       diss_t(k) = - ABS( w_comp ) * (                                        &
     1270                   3.0 * ( v(k+1,j,i) - v(k,j,i)   )                          &
     1271                 -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
    12651272
    12661273       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1267                                  - ( flux_d + diss_d )     ) * ddzw(k)
     1274                                 -   flux_d    - diss_d       ) * ddzw(k)
     1275!                                 
    12681276!--    WS5
    1269        DO  k = nzb_v_inner(j,i) + 3, nzt - 2
     1277       DO  k = nzb_v_inner(j,i)+3, nzt-2
    12701278          flux_d    = flux_t(k-1)
    12711279          diss_d    = diss_t(k-1)
    12721280          w_comp    = w(k,j-1,i) + w(k,j,i)
    12731281          flux_t(k) = w_comp * (                                              &
    1274                       37.  * ( v(k+1,j,i) + v(k,j,i)   )                      &           
    1275                       - 8. * ( v(k+2,j,i) + v(k-1,j,i) )                      &
    1276                       +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
    1277           diss_t(k) = - abs(w_comp) * (                                       &
    1278                       10.  * ( v(k+1,j,i) - v(k,j,i)   )                      &
    1279                       - 5. * ( v(k+2,j,i) - v(k-1,j,i) )                      &
    1280                       +      ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
     1282                      37.0 * ( v(k+1,j,i) + v(k,j,i)   )                      &           
     1283                    -  8.0 * ( v(k+2,j,i) + v(k-1,j,i) )                      &
     1284                    +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
     1285          diss_t(k) = - ABS( w_comp ) * (                                     &
     1286                      10.0 * ( v(k+1,j,i) - v(k,j,i)   )                      &
     1287                    -  5.0 * ( v(k+2,j,i) - v(k-1,j,i) )                      &
     1288                          ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
    12811289
    12821290          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    1283                                  - ( flux_d + diss_d )     ) * ddzw(k)
     1291                                 -   flux_d    - diss_d       ) * ddzw(k)
    12841292       ENDDO
    12851293!
    1286 !--    WS3 as an intermediate step
     1294!--    WS3 as an intermediate step (top)
    12871295       k         = nzt - 1
    12881296       flux_d    = flux_t(k-1)
     
    12901298       w_comp    = w(k,j-1,i) + w(k,j,i)
    12911299       flux_t(k) = w_comp * (                                                 &
    1292                    7. * ( v(k+1,j,i) + v(k,j,i)   )                           &   
    1293                    -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
    1294        diss_t(k) = - abs(w_comp) * (                                          &
    1295                    3. * ( v(k+1,j,i) - v(k,j,i) )                             &
    1296                    -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
     1300                   7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &   
     1301                 -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
     1302       diss_t(k) = - ABS( w_comp ) * (                                        &
     1303                   3.0 * ( v(k+1,j,i) - v(k,j,i) )                            &
     1304                 -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
    12971305       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1298                                  - ( flux_d + diss_d )     ) * ddzw(k)
    1299 !
    1300 !--    2nd order scheme
     1306                                 -   flux_d    - diss_d       ) * ddzw(k)
     1307!
     1308!--    2nd order scheme (top)
    13011309       k         = nzt
    13021310       flux_d    = flux_t(k-1)
     
    13081316 
    13091317       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1310                                  - ( flux_d + diss_d )     ) * ddzw(k)
     1318                                 -   flux_d    - diss_d       ) * ddzw(k)
    13111319             
    13121320       DO  k = nzb_v_inner(j,i), nzt
    13131321          sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:)                           &
    1314             + ( flux_t(k) + diss_t(k) )                                       &
    1315             * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     1322                 + ( flux_t(k) + diss_t(k) )                                  &
     1323                 *  weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    13161324       ENDDO
    13171325
     
    13521360         
    13531361            CASE ( 1 )
    1354                DO  k = nzb_w_inner(j,i) + 1, nzt
     1362               DO  k = nzb_w_inner(j,i)+1, nzt
    13551363                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    13561364                  flux_r(k) = u_comp * (                                      &
    1357                               7. * ( w(k,j,i+1) + w(k,j,i)   )                &
    1358                               -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
    1359                   diss_r(k) = -abs(u_comp) * (                                &
    1360                               3. * ( w(k,j,i+1) - w(k,j,i)   )                &
    1361                               -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
     1365                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
     1366                            -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
     1367                  diss_r(k) = -ABS( u_comp ) * (                              &
     1368                              3.0 * ( w(k,j,i+1) - w(k,j,i)   )               &
     1369                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
    13621370               ENDDO
    13631371               
    13641372            CASE ( 2 )
    1365                DO  k = nzb_w_inner(j,i) + 1, nzt
     1373               DO  k = nzb_w_inner(j,i)+1, nzt
    13661374                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    13671375                  flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25
     
    13721380               
    13731381            CASE ( 3 )
    1374                DO  k = nzb_w_inner(j,i) + 1, nzt
     1382               DO  k = nzb_w_inner(j,i)+1, nzt
    13751383                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    13761384                  flux_n(k) = v_comp * (                                      &
    1377                               7. * ( w(k,j+1,i) + w(k,j,i)   )                &
    1378                               -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    1379                   diss_n(k) = -abs(v_comp) * (                                &
    1380                               3. * ( w(k,j+1,i) - w(k,j,i)   )                &
    1381                               -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
     1385                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
     1386                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
     1387                  diss_n(k) = -ABS( v_comp ) * (                              &
     1388                              3.0 * ( w(k,j+1,i) - w(k,j,i)   )               &
     1389                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
    13821390               ENDDO
    13831391               
    13841392            CASE ( 4 )
    1385                DO  k = nzb_w_inner(j,i) + 1, nzt
     1393               DO  k = nzb_w_inner(j,i)+1, nzt
    13861394                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    13871395                  flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25
     
    13921400               
    13931401            CASE ( 5 )
    1394                DO  k = nzb_w_inner(j,i) + 1, nzt
     1402               DO  k = nzb_w_inner(j,i)+1, nzt
    13951403                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    13961404                  flux_r(k) = u_comp * (                                      &
    1397                               7. * ( w(k,j,i+1) + w(k,j,i)   )                &
    1398                               -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
    1399                   diss_r(k) = - abs(u_comp) * (                               &
    1400                               3. * ( w(k,j,i+1) - w(k,j,i) )                  &
    1401                               -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
     1405                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
     1406                            -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
     1407                  diss_r(k) = - ABS( u_comp ) * (                             &
     1408                              3.0 * ( w(k,j,i+1) - w(k,j,i) )                 &
     1409                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
    14021410               ENDDO
    14031411               
    14041412            CASE ( 6 )
    1405                DO  k = nzb_w_inner(j,i) + 1, nzt
     1413               DO  k = nzb_w_inner(j,i)+1, nzt
    14061414!               
    14071415!--               Compute leftside fluxes for the left boundary of PE domain
    14081416                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    14091417                  flux_r(k) = u_comp *(                                       &
    1410                               7. * ( w(k,j,i+1) + w(k,j,i)   )                &
    1411                               -    ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
    1412                   diss_r(k) = - abs(u_comp) * (                               &
    1413                               3. * ( w(k,j,i+1) - w(k,j,i) )                  &
    1414                               -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
     1418                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
     1419                            -       ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
     1420                  diss_r(k) = - ABS( u_comp ) * (                             &
     1421                              3.0 * ( w(k,j,i+1) - w(k,j,i) )                 &
     1422                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
    14151423                 
    14161424                  u_comp        = u(k+1,j,i) + u(k,j,i) - gu
     
    14181426                  diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
    14191427                                            w(k,j,i-1), w(k,j,i-1), u_comp,   &
    1420                                             0.25,ddx)
     1428                                            0.25, ddx )
    14211429               ENDDO
    14221430               degraded_l = .TRUE.
    14231431               
    14241432            CASE ( 7 )
    1425                DO  k = nzb_w_inner(j,i) + 1, nzt
     1433               DO  k = nzb_w_inner(j,i)+1, nzt
    14261434                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    14271435                  flux_n(k) = v_comp *(                                       &
    1428                               7. * ( w(k,j+1,i) + w(k,j,i)   )                &
    1429                               -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    1430                   diss_n(k) = - abs(v_comp) * (                               &
    1431                               3. * ( w(k,j+1,i) - w(k,j,i)   )                &
    1432                               -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
     1436                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
     1437                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
     1438                  diss_n(k) = - ABS( v_comp ) * (                             &
     1439                              3.0 * ( w(k,j+1,i) - w(k,j,i)   )               &
     1440                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
    14331441                ENDDO
    14341442               
    14351443            CASE ( 8 )
    1436                DO  k = nzb_w_inner(j,i) + 1, nzt
     1444               DO  k = nzb_w_inner(j,i)+1, nzt
    14371445                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    14381446                  flux_n(k) = v_comp * (                                      &
    1439                               7. * ( w(k,j+1,i) + w(k,j,i)   )                &
    1440                              -     ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    1441                   diss_n(k) = - abs(v_comp) * (                               &
    1442                               3. * ( w(k,j+1,i) - w(k,j,i) )                  &
    1443                               -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
     1447                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
     1448                                 ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
     1449                  diss_n(k) = - ABS( v_comp ) * (                             &
     1450                              3.0 * ( w(k,j+1,i) - w(k,j,i) )                 &
     1451                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
    14441452!                   
    14451453!--              Compute southside fluxes for the south boundary of PE domain           
     
    14571465!         
    14581466!--      Compute the crosswise 5th order fluxes at the outflow
    1459          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
    1460          .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
     1467         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
     1468             boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
    14611469         
    1462             DO  k = nzb_w_inner(j,i) + 1, nzt
     1470            DO  k = nzb_w_inner(j,i)+1, nzt
    14631471               v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    14641472               flux_n(k) = v_comp * (                                         &
    1465                            37.  * ( w(k,j+1,i) + w(k,j,i)   )                 &
    1466                            - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                 &
    1467                            +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    1468                diss_n(k) = - abs(v_comp) * (                                  &
    1469                            10.  * ( w(k,j+1,i) - w(k,j,i)   )                 &
    1470                            - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                 &
    1471                            +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     1473                           37.0 * ( w(k,j+1,i) + w(k,j,i)   )                 &
     1474                         -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                 &
     1475                               ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
     1476               diss_n(k) = - ABS( v_comp ) * (                                &
     1477                           10.0 * ( w(k,j+1,i) - w(k,j,i)   )                 &
     1478                         -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                 &
     1479                               ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    14721480            ENDDO
    14731481           
    14741482         ELSE
    14751483         
    1476             DO  k = nzb_w_inner(j,i) + 1, nzt         
     1484            DO  k = nzb_w_inner(j,i)+1, nzt         
    14771485               u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    14781486               flux_r(k) = u_comp * (                                         &
    1479                            37.  * ( w(k,j,i+1) + w(k,j,i) )                   &
    1480                            - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                 &
    1481                            +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    1482                diss_r(k) = - abs(u_comp) * (                                  &
    1483                            10.  * ( w(k,j,i+1) - w(k,j,i)   )                 &
    1484                            - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                 &
    1485                            +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     1487                           37.0 * ( w(k,j,i+1) + w(k,j,i) )                   &
     1488                         -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                 &
     1489                               ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
     1490               diss_r(k) = - ABS( u_comp ) * (                                &
     1491                           10.0 * ( w(k,j,i+1) - w(k,j,i)   )                 &
     1492                         -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                 &
     1493                               ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    14861494            ENDDO
    14871495           
     
    14911499!       
    14921500!--      Compute the fifth order fluxes for the interior of PE domain.               
    1493          DO  k = nzb_w_inner(j,i) + 1, nzt
     1501         DO  k = nzb_w_inner(j,i)+1, nzt
    14941502            u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    14951503            flux_r(k) = u_comp * (                                            &
    1496                         37.  * ( w(k,j,i+1) + w(k,j,i)   )                    &
    1497                         - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                    &
    1498                         +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    1499             diss_r(k) = - abs(u_comp) * (                                     &
    1500                         10.  * ( w(k,j,i+1) - w(k,j,i)   )                    &
    1501                         - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                    &
    1502                         +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     1504                        37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
     1505                      -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
     1506                            ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
     1507            diss_r(k) = - ABS( u_comp ) * (                                   &
     1508                        10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
     1509                      -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
     1510                            ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    15031511                 
    15041512            v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    15051513            flux_n(k) = v_comp * (                                            &
    1506                         37.  * ( w(k,j+1,i) + w(k,j,i)   )                    &
    1507                         - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                    &
    1508                         +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    1509             diss_n(k) = - abs(v_comp) * (                                     &
    1510                         10.  * ( w(k,j+1,i) - w(k,j,i)   )                    &
    1511                         - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                    &
    1512                         +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     1514                        37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
     1515                      -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
     1516                            ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
     1517            diss_n(k) = - ABS( v_comp ) * (                                   &
     1518                        10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
     1519                      -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
     1520                            ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    15131521         ENDDO
    15141522         
     
    15181526       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
    15191527       
    1520           DO  k = nzb_w_inner(j,i) + 1, nzt
     1528          DO  k = nzb_w_inner(j,i)+1, nzt
    15211529             v_comp      = v(k+1,j,i) + v(k,j,i) - gv
    15221530             flux_s_w(k) = v_comp * (                                         &
    1523                            37.  * ( w(k,j,i) + w(k,j-1,i)   )                 &
    1524                            - 8. * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
    1525                            +      ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    1526              diss_s_w(k) = - abs(v_comp) * (                                  &
    1527                            10.  * ( w(k,j,i) - w(k,j-1,i)   )                 &
    1528                            - 5. * ( w(k,j+1,i) - w(k,j-2,i) )                 &
    1529                            +      ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
     1531                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
     1532                         -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
     1533                               ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
     1534             diss_s_w(k) = - ABS( v_comp ) * (                                &
     1535                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
     1536                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
     1537                               ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
    15301538          ENDDO
    15311539         
     
    15341542       IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN
    15351543       
    1536           DO  k = nzb_w_inner(j,i) + 1, nzt
     1544          DO  k = nzb_w_inner(j,i)+1, nzt
    15371545            u_comp        = u(k+1,j,i) + u(k,j,i) - gu
    15381546            flux_l_w(k,j) = u_comp * (                                        &
    1539                             37.  * ( w(k,j,i) + w(k,j,i-1)   )                &
    1540                             - 8. * ( w(k,j,i+1) + w(k,j,i-2) )                &
    1541                             +      ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    1542             diss_l_w(k,j) = - abs(u_comp) * (                                 &
    1543                             10.  * ( w(k,j,i) - w(k,j,i-1)   )                &
    1544                             - 5. * ( w(k,j,i+1) - w(k,j,i-2) )                &
    1545                             +      ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
     1547                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
     1548                          -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
     1549                                ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
     1550            diss_l_w(k,j) = - ABS( u_comp ) * (                               &
     1551                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
     1552                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
     1553                                ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
    15461554          ENDDO
    15471555         
     
    15491557!       
    15501558!--    Now compute the tendency terms for the horizontal parts.         
    1551        DO  k = nzb_w_inner(j,i) + 1, nzt
     1559       DO  k = nzb_w_inner(j,i)+1, nzt
    15521560          tend(k,j,i) = tend(k,j,i) - (                                       &
    15531561                      ( flux_r(k) + diss_r(k)                                 &
    1554                     - ( flux_l_w(k,j) + diss_l_w(k,j) ) ) * ddx               &
     1562                    -   flux_l_w(k,j) - diss_l_w(k,j)  ) * ddx               &
    15551563                    + ( flux_n(k) + diss_n(k)                                 &
    1556                     - ( flux_s_w(k) + diss_s_w(k)     ) ) * ddy  )
     1564                    -   flux_s_w(k) - diss_s_w(k)      ) * ddy  )
    15571565
    15581566          flux_l_w(k,j) = flux_r(k)
     
    15661574!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    15671575!--    statistical evaluation the top flux at the surface should be 0
    1568        flux_t(nzb_w_inner(j,i)) = 0.  !statistical reasons
    1569        diss_t(nzb_w_inner(j,i)) = 0.
    1570 !
    1571 !--    2nd order scheme        
    1572        k         = nzb_w_inner(j,i) + 1
     1576       flux_t(nzb_w_inner(j,i)) = 0.0  !statistical reasons
     1577       diss_t(nzb_w_inner(j,i)) = 0.0
     1578!
     1579!--    2nd order scheme (bottom)       
     1580       k         = nzb_w_inner(j,i)+1
    15731581       flux_d    = flux_t(k-1)
    15741582       diss_d    = diss_t(k-1)
    15751583       w_comp    = w(k+1,j,i) + w(k,j,i)
    15761584       flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
    1577        diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0.,        &
     1585       diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0,        &
    15781586                             w_comp, 0.25, ddzu(k+1) )
    15791587                     
    15801588       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1581                                  - ( flux_d + diss_d )     ) * ddzu(k+1)               
    1582 !
    1583 !--    WS3 as an intermediate step
    1584        k         = nzb_w_inner(j,i) + 2
     1589                                 -   flux_d    - diss_d      ) * ddzu(k+1)               
     1590!
     1591!--    WS3 as an intermediate step (bottom)
     1592       k         = nzb_w_inner(j,i)+2
    15851593       flux_d    = flux_t(k-1)
    15861594       diss_d    = diss_t(k-1)
    15871595       w_comp    = w(k+1,j,i) + w(k,j,i)
    15881596       flux_t(k) = w_comp * (                                                 &
    1589                    7. * ( w(k+1,j,i) + w(k,j,i) )                             &
    1590                    -    ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
    1591        diss_t(k) = - abs(w_comp) * (                                          &
    1592                    3. * ( w(k+1,j,i) - w(k,j,i) )                             &
    1593                    -    ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
     1597                   7.0 * ( w(k+1,j,i) + w(k,j,i) )                            &
     1598                 -       ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
     1599       diss_t(k) = - ABS( w_comp ) * (                                        &
     1600                   3.0 * ( w(k+1,j,i) - w(k,j,i) )                            &
     1601                 -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
    15941602
    15951603       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1596                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
     1604                                 -   flux_d    - diss_d      ) * ddzu(k+1)
    15971605!                                 
    15981606!--    WS5
    1599        DO  k = nzb_w_inner(j,i) + 3, nzt -2
     1607       DO  k = nzb_w_inner(j,i)+3, nzt-2
    16001608          flux_d    = flux_t(k-1)
    16011609          diss_d    = diss_t(k-1)
    16021610          w_comp    = w(k+1,j,i) + w(k,j,i)
    16031611          flux_t(k) = w_comp * (                                              &
    1604                       37.  * ( w(k+1,j,i) + w(k,j,i)   )                      &
    1605                       - 8. * ( w(k+2,j,i) + w(k-1,j,i) )                      &
    1606                       +      ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
    1607           diss_t(k) = - abs(w_comp) * (                                       &
    1608                       10.  * ( w(k+1,j,i) - w(k,j,i)   )                      &
    1609                       - 5. * ( w(k+2,j,i) - w(k-1,j,i) )                      &
    1610                       +      ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
     1612                      37.0 * ( w(k+1,j,i) + w(k,j,i)   )                      &
     1613                    -  8.0 * ( w(k+2,j,i) + w(k-1,j,i) )                      &
     1614                          ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
     1615          diss_t(k) = - ABS( w_comp ) * (                                     &
     1616                      10.0 * ( w(k+1,j,i) - w(k,j,i)   )                      &
     1617                    -  5.0 * ( w(k+2,j,i) - w(k-1,j,i) )                      &
     1618                          ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
    16111619
    16121620          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    1613                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
     1621                                    -   flux_d    - diss_d      ) * ddzu(k+1)
    16141622       ENDDO
    1615 !--    WS3 as an intermediate step
     1623!--    WS3 as an intermediate step (top)
    16161624       k         = nzt - 1
    16171625       flux_d    = flux_t(k-1)
     
    16191627       w_comp    = w(k+1,j,i) + w(k,j,i)
    16201628       flux_t(k) = w_comp * (                                                 &
    1621                    7. * ( w(k+1,j,i) + w(k,j,i)   )                           &
    1622                    -    ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
    1623        diss_t(k) = - abs(w_comp) * (                                          &
    1624                    3. * ( w(k+1,j,i) - w(k,j,i)     )                         &
    1625                    -    ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
     1629                   7.0 * ( w(k+1,j,i) + w(k,j,i)   )                          &
     1630                 -       ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
     1631       diss_t(k) = - ABS( w_comp ) * (                                        &
     1632                   3.0 * ( w(k+1,j,i) - w(k,j,i)     )                        &
     1633                 -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
    16261634                   
    16271635       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1628                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
    1629 !
    1630 !--    2nd order scheme
     1636                                 -   flux_d    - diss_d      ) * ddzu(k+1)
     1637!
     1638!--    2nd order scheme (top)
    16311639       k         = nzt
    16321640       flux_d    = flux_t(k-1)
     
    16381646
    16391647       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1640                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
     1648                                 -   flux_d    - diss_d      ) * ddzu(k+1)
    16411649       
    16421650       DO  k = nzb_w_inner(j,i), nzt
    16431651           sums_ws2_ws_l(k,:)  = sums_ws2_ws_l(k,:)                           &
    1644                + ( flux_t(k) + diss_t(k) )                                    &
    1645                * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     1652                 + ( flux_t(k) + diss_t(k) )                                  &
     1653                 *  weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    16461654       ENDDO
    16471655
     
    16491657   
    16501658
    1651 !
     1659!------------------------------------------------------------------------------!
    16521660! Scalar advection - Call for all grid points
    16531661!------------------------------------------------------------------------------!
     
    16801688          IF ( boundary_flags(j,i) == 6 )  THEN
    16811689         
    1682              DO  k = nzb_s_inner(j,i) + 1, nzt
     1690             DO  k = nzb_s_inner(j,i)+1, nzt
    16831691                u_comp                 = u(k,j,i) - u_gtrans
    16841692                swap_flux_x_local(k,j) = u_comp * (                           &
     
    16921700          ELSE
    16931701         
    1694              DO  k = nzb_s_inner(j,i) + 1, nzt
     1702             DO  k = nzb_s_inner(j,i)+1, nzt
    16951703                u_comp                 = u(k,j,i) - u_gtrans
    1696                 swap_flux_x_local(k,j) = u_comp*(                             &
    1697                                          37.  * (sk(k,j,i)+sk(k,j,i-1)    )   &
    1698                                          -8. * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
    1699                                          +     ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
    1700                                          * adv_sca_5
    1701                 swap_diss_x_local(k,j) = - abs(u_comp) * (                    &
    1702                                          10. * (sk(k,j,i) - sk(k,j,i-1)    )  &
    1703                                          -5. * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
    1704                                          +     ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
    1705                                          * adv_sca_5
     1704                swap_flux_x_local(k,j) = u_comp*(                              &
     1705                                         37.0 * ( sk(k,j,i)+sk(k,j,i-1)     )  &
     1706                                       -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
     1707                                       +        ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
     1708                                       * adv_sca_5
     1709                swap_diss_x_local(k,j) = - ABS( u_comp ) * (                   &
     1710                                         10.0 * (sk(k,j,i) - sk(k,j,i-1)    )  &
     1711                                       -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
     1712                                       +        ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
     1713                                       * adv_sca_5
    17061714             ENDDO
    17071715          ENDIF
     
    17151723          IF ( boundary_flags(j,i) == 8 )  THEN
    17161724         
    1717              DO  k = nzb_s_inner(j,i) + 1, nzt
     1725             DO  k = nzb_s_inner(j,i)+1, nzt
    17181726                v_comp               = v(k,j,i) - v_gtrans
    1719                 swap_flux_y_local(k) = v_comp *                               &
     1727                swap_flux_y_local(k) = v_comp *                                &
    17201728                                       ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
    1721                 swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),    &
    1722                                                  sk(k,j,i), sk(k,j-1,i),      &
    1723                                                  sk(k,j-1,i), v_comp, 0.5, ddy)
     1729                swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),     &
     1730                                                 sk(k,j,i), sk(k,j-1,i),       &
     1731                                                 sk(k,j-1,i), v_comp, 0.5, ddy )
    17241732             ENDDO
    17251733             
    17261734          ELSE
    17271735         
    1728              DO  k = nzb_s_inner(j,i) + 1, nzt
     1736             DO  k = nzb_s_inner(j,i)+1, nzt
    17291737                v_comp               = v(k,j,i) - v_gtrans
    1730                 swap_flux_y_local(k) = v_comp * (                             &
    1731                                        37.  * ( sk(k,j,i) + sk(k,j-1,i)   )   &
    1732                                        - 8. * ( sk(k,j+1,i) + sk(k,j-2,i) )   &
    1733                                        +      ( sk(k,j+2,i) + sk(k,j-3,i) ) )
    1734                                        * adv_sca_5
    1735                 swap_diss_y_local(k)= - abs(v_comp) * (                       &
    1736                                       10.  * ( sk(k,j,i) - sk(k,j-1,i)   )    &
    1737                                       - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
    1738                                       +      ( sk(k,j+2,i)-sk(k,j-3,i) ) )    &
    1739                                       * adv_sca_5
     1738                swap_flux_y_local(k) = v_comp * (                              &
     1739                                       37.0 * ( sk(k,j,i) + sk(k,j-1,i)   )    &
     1740                                     -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
     1741                                     +        ( sk(k,j+2,i) + sk(k,j-3,i) ) )
     1742                                     * adv_sca_5
     1743                swap_diss_y_local(k)= - ABS( v_comp ) * (                      &
     1744                                      10.0 * ( sk(k,j,i) - sk(k,j-1,i)   )     &
     1745                                    -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
     1746                                    +        ( sk(k,j+2,i)-sk(k,j-3,i) ) )     &
     1747                                    * adv_sca_5
    17401748             ENDDO
    17411749             
     
    17491757               
    17501758                  CASE ( 1 )
    1751                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1759                     DO  k = nzb_s_inner(j,i)+1, nzt
    17521760                        u_comp    = u(k,j,i+1) - u_gtrans
    1753                         flux_r(k) = u_comp * (                                &
    1754                                   7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
    1755                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    1756                         diss_r(k) = - abs(u_comp) * (                         &
    1757                                   3. * ( sk(k,j,i+1) - sk(k,j,i) )            &
    1758                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     1761                        flux_r(k) = u_comp * (                                 &
     1762                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
     1763                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
     1764                                  * adv_sca_3
     1765                        diss_r(k) = - ABS( u_comp ) * (                        &
     1766                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        &
     1767                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
     1768                                  * adv_sca_3
    17591769                     ENDDO
    17601770                     
    17611771                  CASE ( 2 )
    1762                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1772                     DO  k = nzb_s_inner(j,i)+1, nzt
    17631773                        u_comp    = u(k,j,i+1) - u_gtrans
    17641774                        flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
    1765                         diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),       &
    1766                                               sk(k,j,i), sk(k,j,i-1),         &
     1775                        diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),        &
     1776                                              sk(k,j,i), sk(k,j,i-1),          &
    17671777                                              sk(k,j,i-2), u_comp, 0.5, ddx )
    17681778                     ENDDO
    17691779                     
    17701780                  CASE ( 3 )
    1771                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1781                     DO  k = nzb_s_inner(j,i)+1, nzt
    17721782                        v_comp    = v(k,j+1,i) - v_gtrans
    1773                         flux_n(k) = v_comp * (                                &
    1774                                   7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
    1775                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    1776                         diss_n(k) = - abs(v_comp) * (                         &
    1777                                   3. * ( sk(k,j+1,i) - sk(k,j,i) )            &
    1778                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     1783                        flux_n(k) = v_comp * (                                 &
     1784                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
     1785                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
     1786                                  * adv_sca_3
     1787                        diss_n(k) = - ABS( v_comp ) * (                        &
     1788                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        &
     1789                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
     1790                                  * adv_sca_3
    17791791                     ENDDO
    17801792                     
    17811793                  CASE ( 4 )
    1782                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1794                     DO  k = nzb_s_inner(j,i)+1, nzt
    17831795                        v_comp    = v(k,j+1,i) - v_gtrans
    17841796                        flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
    1785                         diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),       &
    1786                                               sk(k,j,i), sk(k,j-1,i),         &
     1797                        diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),        &
     1798                                              sk(k,j,i), sk(k,j-1,i),          &
    17871799                                              sk(k,j-2,i), v_comp, 0.5, ddy )     
    17881800                     ENDDO
    17891801                     
    17901802                  CASE ( 5 )
    1791                      DO  k = nzb_w_inner(j,i) + 1, nzt
     1803                     DO  k = nzb_w_inner(j,i)+1, nzt
    17921804                        u_comp    = u(k,j,i+1) - u_gtrans
    1793                         flux_r(k) = u_comp * (                                &
    1794                                   7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
    1795                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    1796                         diss_r(k) = - abs(u_comp) * (                         &
    1797                                   3. * ( sk(k,j,i+1) - sk(k,j,i) )            &
    1798                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     1805                        flux_r(k) = u_comp * (                                 &
     1806                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
     1807                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
     1808                                  * adv_sca_3
     1809                        diss_r(k) = - ABS( u_comp ) * (                        &
     1810                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        &
     1811                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
     1812                                  * adv_sca_3
    17991813                     ENDDO 
    18001814                     
    18011815                  CASE ( 6 )
    1802                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1816                     DO  k = nzb_s_inner(j,i)+1, nzt
    18031817                        u_comp    = u(k,j,i+1) - u_gtrans
    1804                         flux_r(k) = u_comp * (                                &
    1805                                   7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
    1806                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    1807                         diss_r(k) = - abs(u_comp) * (                         &
    1808                                   3. * ( sk(k,j,i+1) - sk(k,j,i) )            &
    1809                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     1818                        flux_r(k) = u_comp * (                                 &
     1819                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
     1820                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
     1821                                  * adv_sca_3
     1822                        diss_r(k) = - ABS( u_comp ) * (                        &
     1823                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        &
     1824                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
     1825                                  * adv_sca_3
    18101826                     ENDDO
    18111827                     
    18121828                  CASE ( 7 )
    1813                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1829                     DO  k = nzb_s_inner(j,i)+1, nzt
    18141830                        v_comp    = v(k,j+1,i) - v_gtrans
    1815                         flux_n(k) = v_comp * (                                &
    1816                                   7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
    1817                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    1818                         diss_n(k) = - abs(v_comp) * (                         &
    1819                                   3. * ( sk(k,j+1,i) - sk(k,j,i) )            &
    1820                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     1831                        flux_n(k) = v_comp * (                                 &
     1832                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
     1833                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
     1834                                  * adv_sca_3
     1835                        diss_n(k) = - ABS( v_comp ) * (                        &
     1836                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        &
     1837                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
     1838                                  * adv_sca_3
    18211839                     ENDDO
    18221840                     
    18231841                  CASE ( 8 )
    1824                      DO  k = nzb_s_inner(j,i) + 1, nzt
     1842                     DO  k = nzb_s_inner(j,i)+1, nzt
    18251843                        v_comp    = v(k,j+1,i) - v_gtrans
    1826                         flux_n(k) = v_comp * (                                &
    1827                                   7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
    1828                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    1829                         diss_n(k) = -  abs(v_comp) * (                        &
    1830                                   3. * ( sk(k,j+1,i) - sk(k,j,i) )            &
    1831                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     1844                        flux_n(k) = v_comp * (                                 &
     1845                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
     1846                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
     1847                                  * adv_sca_3
     1848                        diss_n(k) = -  ABS( v_comp ) * (                       &
     1849                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        &
     1850                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
     1851                                  * adv_sca_3
    18321852                     ENDDO 
    18331853                     
     
    18361856               END SELECT
    18371857               
    1838                IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2     &
    1839              .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6)  THEN
     1858               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
     1859                   boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )    &
     1860               THEN
    18401861             
    1841                   DO  k = nzb_s_inner(j,i) + 1, nzt
     1862                  DO  k = nzb_s_inner(j,i)+1, nzt
    18421863                     v_comp    = v(k,j+1,i) - v_gtrans
    1843                      flux_n(k) = v_comp * (                                   &
    1844                                  37.  * ( sk(k,j+1,i) + sk(k,j,i)  )          &
    1845                                  - 8. * (sk(k,j+2,i) + sk(k,j-1,i) )          &
    1846                                  +      ( sk(k,j+3,i) + sk(k,j-2,i) ) )       &
    1847                                  * adv_sca_5
    1848                      diss_n(k) = - abs(v_comp) * (                            &
    1849                                  10.  * ( sk(k,j+1,i) - sk(k,j,i)   )         &
    1850                                  - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )         &
    1851                                  +      ( sk(k,j+3,i) - sk(k,j-2,i) ) )       &
    1852                                  * adv_sca_5
     1864                     flux_n(k) = v_comp * (                                    &
     1865                                 37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )          &
     1866                               -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )          &
     1867                               +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )        &
     1868                               * adv_sca_5
     1869                     diss_n(k) = - ABS( v_comp ) * (                           &
     1870                                 10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )          &
     1871                               -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )          &
     1872                               +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )        &
     1873                               * adv_sca_5
    18531874                  ENDDO
    18541875                 
    18551876               ELSE
    18561877               
    1857                   DO  k = nzb_s_inner(j,i) + 1, nzt
     1878                  DO  k = nzb_s_inner(j,i)+1, nzt
    18581879                     u_comp    = u(k,j,i+1) - u_gtrans 
    1859                      flux_r(k) = u_comp * (                                   &
    1860                                  37.  * ( sk(k,j,i+1) + sk(k,j,i)   )         &
    1861                                  - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )         &
    1862                                  +      ( sk(k,j,i+3) + sk(k,j,i-2) ) )       &
    1863                                  * adv_sca_5
    1864                      diss_r(k) = - abs(u_comp) * (                            &
    1865                                  10.  * ( sk(k,j,i+1) - sk(k,j,i)   )         &
    1866                                  - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )         &
    1867                                  +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )       &
    1868                                  * adv_sca_5
     1880                     flux_r(k) = u_comp * (                                    &
     1881                                 37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )          &
     1882                               -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )          &
     1883                               +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )        &
     1884                               * adv_sca_5
     1885                     diss_r(k) = - ABS( u_comp ) * (                           &
     1886                                 10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )          &
     1887                               -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )          &
     1888                               +        ( sk(k,j,i+3) - sk(k,j,i-2) ) )        &
     1889                               * adv_sca_5
    18691890                  ENDDO
    18701891                 
     
    18731894            ELSE
    18741895           
    1875                DO  k = nzb_s_inner(j,i) + 1, nzt
     1896               DO  k = nzb_s_inner(j,i)+1, nzt
    18761897                  u_comp    = u(k,j,i+1) - u_gtrans
    1877                   flux_r(k) = u_comp * (                                      &
    1878                               37.  * ( sk(k,j,i+1) + sk(k,j,i)   )            &
    1879                               - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )            &
    1880                               +      ( sk(k,j,i+3) + sk(k,j,i-2) ) )          &
    1881                               * adv_sca_5
    1882                   diss_r(k) = - abs(u_comp) * (                               &
    1883                               10.  * ( sk(k,j,i+1) - sk(k,j,i)   )            &
    1884                               - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )            &
    1885                               +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )          &
    1886                               * adv_sca_5
     1898                  flux_r(k) = u_comp * (                                       &
     1899                              37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     1900                            -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )             &
     1901                            +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )           &
     1902                            * adv_sca_5
     1903                  diss_r(k) = - ABS( u_comp ) * (                              &
     1904                              10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     1905                            -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )             &
     1906                            +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )             &
     1907                            * adv_sca_5
    18871908                 
    18881909                  v_comp    = v(k,j+1,i) - v_gtrans
    1889                   flux_n(k) = v_comp * (                                      &
    1890                               37.  * ( sk(k,j+1,i) + sk(k,j,i)   )            &
    1891                               - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )            &
    1892                               +      ( sk(k,j+3,i) + sk(k,j-2,i) ) )          &
    1893                               * adv_sca_5
    1894                   diss_n(k) = - abs(v_comp) * (                               &
    1895                               10.  * ( sk(k,j+1,i) - sk(k,j,i)   )            &
    1896                               - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )            &
    1897                               +      ( sk(k,j+3,i) - sk(k,j-2,i) ) )          &
    1898                               * adv_sca_5                 
     1910                  flux_n(k) = v_comp * (                                       &
     1911                              37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     1912                            -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )             &
     1913                            +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )           &
     1914                            * adv_sca_5
     1915                  diss_n(k) = - ABS( v_comp ) * (                              &
     1916                              10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
     1917                            -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )             &
     1918                            +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )           &
     1919                            * adv_sca_5                 
    18991920               ENDDO
    19001921               
    19011922            ENDIF
    19021923           
    1903             DO  k = nzb_s_inner(j,i) + 1, nzt
     1924            DO  k = nzb_s_inner(j,i)+1, nzt
    19041925               tend(k,j,i) = tend(k,j,i) - (                                  &
    19051926                  ( flux_r(k) + diss_r(k)                                     &
    1906                 - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx &
     1927                -   swap_flux_x_local(k,j) - swap_diss_x_local(k,j)  ) * ddx &
    19071928                + ( flux_n(k) + diss_n(k)                                     &
    1908                 - ( swap_flux_y_local(k) + swap_diss_y_local(k) )     ) * ddy)
     1929                -   swap_flux_y_local(k) - swap_diss_y_local(k)       ) * ddy)
    19091930                   
    19101931                swap_flux_x_local(k,j) = flux_r(k)
     
    19231944         DO  j = nys, nyn
    19241945!
    1925 !--         2nd order scheme
    1926             k=nzb_s_inner(j,i) + 1
     1946!--         2nd order scheme (bottom)
     1947            k=nzb_s_inner(j,i)+1
    19271948!           
    19281949!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
    19291950!--         reasons the top flux at the surface should be 0.
    1930             flux_t(nzb_s_inner(j,i)) = 0.
    1931             diss_t(nzb_s_inner(j,i)) = 0.
     1951            flux_t(nzb_s_inner(j,i)) = 0.0
     1952            diss_t(nzb_s_inner(j,i)) = 0.0
    19321953            flux_d = flux_t(k-1)
    19331954            diss_d = diss_t(k-1)
    1934             flux_t(k) = w(k,j,i) * (sk(k+1,j,i) + sk(k,j,i) ) *0.5
     1955            flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
    19351956            diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i),        &
    19361957                                  sk(k,j,i), sk(k,j,i), w(k,j,i),             &
    19371958                                  0.5, ddzw(k) )
    19381959            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    1939                                       - ( flux_d + diss_d )     ) * ddzw(k)   
    1940 !
    1941 !--         WS3 as an intermediate step
    1942             k         = nzb_s_inner(j,i) + 2
     1960                                      -   flux_d    - diss_d       ) * ddzw(k)   
     1961!
     1962!--         WS3 as an intermediate step (bottom)
     1963            k         = nzb_s_inner(j,i)+2
    19431964            flux_d    = flux_t(k-1)
    19441965            diss_d    = diss_t(k-1)
    19451966            flux_t(k) = w(k,j,i) * (                                          &
    1946                         7. * ( sk(k+1,j,i) + sk(k,j,i) )                      &
     1967                        7.0 * ( sk(k+1,j,i) + sk(k,j,i) )                     &
    19471968                       - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    1948             diss_t(k) = - abs(w(k,j,i)) * (                                   &
    1949                         3. * ( sk(k+1,j,i) - sk(k,j,i) )                      &
    1950                         - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
     1969            diss_t(k) = - ABS( w(k,j,i) ) * (                                 &
     1970                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   &
     1971                      -      ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    19511972            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    1952                                       - ( flux_d + diss_d )     ) * ddzw(k)
     1973                                      -   flux_d    - diss_d       ) * ddzw(k)
    19531974!
    19541975!--         WS5
    1955             DO  k = nzb_s_inner(j,i) + 3, nzt - 2
     1976            DO  k = nzb_s_inner(j,i)+3, nzt-2
    19561977               flux_d    = flux_t(k-1)
    19571978               diss_d    = diss_t(k-1)
    19581979               flux_t(k) = w(k,j,i) * (                                       &
    1959                            37.  * ( sk(k+1,j,i) + sk(k,j,i)   )               &
    1960                            - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
    1961                            +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
    1962                diss_t(k) = - abs(w(k,j,i)) * (                                &
    1963                            10.  * ( sk(k+1,j,i) -sk(k,j,i) )                  &
    1964                            - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
    1965                            +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
     1980                           37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )               &
     1981                         -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
     1982                               ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
     1983               diss_t(k) = - ABS(w(k,j,i)) * (                                &
     1984                           10.0 * ( sk(k+1,j,i) -sk(k,j,i)    )               &
     1985                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
     1986                               ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
    19661987
    19671988               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
    1968                                       - ( flux_d + diss_d )     ) * ddzw(k)
     1989                                         -   flux_d    - diss_d     ) * ddzw(k)
    19691990            ENDDO
    19701991!
    1971 !--         WS3 as an intermediate step
     1992!--         WS3 as an intermediate step (top)
    19721993            k         = nzt - 1
    19731994            flux_d    = flux_t(k-1)
    19741995            diss_d    = diss_t(k-1)
    19751996            flux_t(k) = w(k,j,i) * (                                          &
    1976                         7. * ( sk(k+1,j,i) + sk(k,j,i)  )                     &
    1977                         - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    1978             diss_t(k) = - abs(w(k,j,i)) * (                                   &
    1979                         3. * ( sk(k+1,j,i) - sk(k,j,i)   )                    &
    1980                         -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
     1997                        7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )                   &
     1998                      -      ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
     1999            diss_t(k) = - ABS(w(k,j,i)) * (                                   &
     2000                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   &
     2001                      -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    19812002                       
    19822003            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    1983                                       - ( flux_d + diss_d )     ) * ddzw(k)
    1984 !
    1985 !--         2nd order scheme
     2004                                      -   flux_d    - diss_d       ) * ddzw(k)
     2005!
     2006!--         2nd order scheme (top)
    19862007            k         = nzt
    19872008            flux_d    = flux_t(k-1)
     
    19932014
    19942015            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    1995                                       - ( flux_d + diss_d )     ) * ddzw(k)
     2016                                      -   flux_d    - diss_d       ) * ddzw(k)
    19962017!
    19972018!--         evaluation of statistics
     
    20002021               CASE ( 'pt' )
    20012022                 DO  k = nzb_s_inner(j,i), nzt
    2002                    sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                 &
    2003                       + ( flux_t(k) + diss_t(k) )                              &
    2004                       * weight_substep(intermediate_timestep_count)            &
    2005                       * rmask(j,i,:)
     2023                   sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                &
     2024                      + ( flux_t(k) + diss_t(k) )                             &
     2025                      *   weight_substep(intermediate_timestep_count)         &
     2026                      *   rmask(j,i,:)
    20062027                 ENDDO
    20072028               CASE ( 'sa' )
    20082029                 DO  k = nzb_s_inner(j,i), nzt
    2009                    sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                 &
    2010                       + ( flux_t(k) + diss_t(k) )                              &
    2011                       * weight_substep(intermediate_timestep_count)            &
    2012                       * rmask(j,i,:)
     2030                   sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                &
     2031                      + ( flux_t(k) + diss_t(k) )                             &
     2032                      *   weight_substep(intermediate_timestep_count)         &
     2033                      *   rmask(j,i,:)
    20132034                 ENDDO
    20142035               CASE ( 'q' )
    20152036                 DO  k = nzb_s_inner(j,i), nzt
    2016                    sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                   &
    2017                       + ( flux_t(k) + diss_t(k) )                              &
    2018                       * weight_substep(intermediate_timestep_count)            &
    2019                       * rmask(j,i,:)
     2037                   sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                  &
     2038                      + ( flux_t(k) + diss_t(k) )                             &
     2039                      *   weight_substep(intermediate_timestep_count)         &
     2040                      *   rmask(j,i,:)
    20202041                 ENDDO
    20212042
     
    20282049
    20292050
    2030 !
     2051!------------------------------------------------------------------------------!
    20312052! Advection of u - Call for all grid points
    20322053!------------------------------------------------------------------------------!
     
    20592080          IF( boundary_flags(j,i) == 5 )  THEN
    20602081         
    2061              DO  k = nzb_u_inner(j,i) + 1, nzt
     2082             DO  k = nzb_u_inner(j,i)+1, nzt
    20622083                u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
    2063                 swap_flux_x_local_u(k,j) = u_comp(k) *                        &
     2084                swap_flux_x_local_u(k,j) = u_comp(k) *                         &
    20642085                                           ( u(k,j,i) + u(k,j,i-1) ) * 0.25
    2065                 swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),  &
    2066                                                      u(k,j,i), u(k,j,i-1),    &
    2067                                                      u(k,j,i-1), u_comp(k),   &
     2086                swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),   &
     2087                                                     u(k,j,i), u(k,j,i-1),     &
     2088                                                     u(k,j,i-1), u_comp(k),    &
    20682089                                                     0.25, ddx )
    20692090             ENDDO
     
    20712092          ELSE
    20722093         
    2073              DO  k = nzb_u_inner(j,i) + 1, nzt
     2094             DO  k = nzb_u_inner(j,i)+1, nzt
    20742095                u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
    2075                 swap_flux_x_local_u(k,j) = u_comp(k) * (                      &
    2076                                          37.  * ( u(k,j,i) + u(k,j,i-1)   )   &
    2077                                          - 8. * ( u(k,j,i+1) + u(k,j,i-2) )   &
    2078                                          +      (u(k,j,i+2)+u(k,j,i-3) )  )   &
     2096                swap_flux_x_local_u(k,j) = u_comp(k) * (                       &
     2097                                           37.0 * ( u(k,j,i) + u(k,j,i-1)   )  &
     2098                                         -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )  &
     2099                                         +        (u(k,j,i+2)+u(k,j,i-3) )  )  &
    20792100                                         * adv_mom_5
    2080                 swap_diss_x_local_u(k,j) = - abs(u_comp(k)) * (               &
    2081                                          10.  * ( u(k,j,i) - u(k,j,i-1)   )   &
    2082                                          - 5. * ( u(k,j,i+1) - u(k,j,i-2) )   &
    2083                                          +      ( u(k,j,i+2) - u(k,j,i-3) ) ) &
     2101                swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                &
     2102                                           10.0 * ( u(k,j,i) - u(k,j,i-1)   )  &
     2103                                         -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )  &
     2104                                         +        ( u(k,j,i+2) - u(k,j,i-3) ) )&
    20842105                                         * adv_mom_5
    20852106             ENDDO
     
    20962117!         
    20972118!--         Compute southside fluxes for the south boundary of PE domain
    2098             DO  k = nzb_u_inner(j,i) + 1, nzt
     2119            DO  k = nzb_u_inner(j,i)+1, nzt
    20992120               v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
    21002121               swap_flux_y_local_u(k) = v_comp *                              &
     
    21082129         ELSE
    21092130         
    2110             DO  k = nzb_u_inner(j,i) + 1, nzt
     2131            DO  k = nzb_u_inner(j,i)+1, nzt
    21112132               v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
    2112                swap_flux_y_local_u(k) = v_comp * (                            &
    2113                                       37.  * ( u(k,j,i) + u(k,j-1,i)   )      &
    2114                                       - 8. * ( u(k,j+1,i) + u(k,j-2,i) )      &
    2115                                       +      ( u(k,j+2,i) + u(k,j-3,i) ) )    &
     2133               swap_flux_y_local_u(k) = v_comp * (                             &
     2134                                        37.0 * ( u(k,j,i) + u(k,j-1,i)   )     &
     2135                                      -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )     &
     2136                                      +        ( u(k,j+2,i) + u(k,j-3,i) ) )   &
    21162137                                      * adv_mom_5
    2117                swap_diss_y_local_u(k) = - abs(v_comp) * (                     &
    2118                                       10.  * ( u(k,j,i) - u(k,j-1,i)   )      &
    2119                                       - 5. * ( u(k,j+1,i) - u(k,j-2,i) )      &
    2120                                       +      ( u(k,j+2,i) - u(k,j-3,i) ) )    &
     2138               swap_diss_y_local_u(k) = - ABS( v_comp ) * (                    &
     2139                                        10.0 * ( u(k,j,i) - u(k,j-1,i)   )     &
     2140                                      -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )     &
     2141                                      +        ( u(k,j+2,i) - u(k,j-3,i) ) )   &
    21212142                                      * adv_mom_5
    21222143            ENDDO
     
    21322153               
    21332154                  CASE ( 1 )
    2134                      DO  k = nzb_u_inner(j,i) + 1, nzt
     2155                     DO  k = nzb_u_inner(j,i)+1, nzt
    21352156                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2136                         flux_r(k) = ( u_comp(k) - gu ) * (                    &
    2137                                     7. * ( u(k,j,i+1) + u(k,j,i) )            &
    2138                                     - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    2139                         diss_r(k) = - abs(u_comp(k) - gu) * (                 &
    2140                                     3. * ( u(k,j,i+1) - u(k,j,i) )            &
    2141                                     - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
     2157                        flux_r(k) = ( u_comp(k) - gu ) * (                     &
     2158                                    7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
     2159                                  -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
     2160                                  * adv_mom_3
     2161                        diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
     2162                                    3.0 * ( u(k,j,i+1) - u(k,j,i) )            &
     2163                                  -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
     2164                                  * adv_mom_3
    21422165                     ENDDO
    21432166                     
    21442167                  CASE ( 2 )
    2145                      DO  k = nzb_u_inner(j,i) + 1, nzt
     2168                     DO  k = nzb_u_inner(j,i)+1, nzt
    21462169                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
    21472170                        flux_r(k) = ( u_comp(k) - gu ) *                      &
     
    21532176                     
    21542177                  CASE ( 3 )
    2155                      DO  k = nzb_u_inner(j,i) + 1, nzt
     2178                     DO  k = nzb_u_inner(j,i)+1, nzt
    21562179                        v_comp   = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2157                         flux_n(k) = v_comp * (                                &
    2158                                     7. * ( u(k,j+1,i) + u(k,j,i) )            &
    2159                                     - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    2160                         diss_n(k) = - abs(v_comp) * (                         &
    2161                                     3. * ( u(k,j+1,i) - u(k,j,i) )            &
    2162                                     - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
     2180                        flux_n(k) = v_comp * (                                 &
     2181                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
     2182                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
     2183                                  * adv_mom_3
     2184                        diss_n(k) = - ABS( v_comp ) * (                        &
     2185                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          &
     2186                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
     2187                                  * adv_mom_3
    21632188                     ENDDO
    21642189                     
    21652190                  CASE ( 4 )
    2166                      DO  k = nzb_u_inner(j,i) + 1, nzt
     2191                     DO  k = nzb_u_inner(j,i)+1, nzt
    21672192                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    21682193                        flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
     
    21732198                     
    21742199                  CASE ( 5 )
    2175                      DO  k = nzb_u_inner(j,i) + 1, nzt       
     2200                     DO  k = nzb_u_inner(j,i)+1, nzt       
    21762201                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2177                         flux_r(k) = ( u_comp(k) - gu ) * (                    &
    2178                                     7. * ( u(k,j,i+1) + u(k,j,i) )            &
    2179                                     - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    2180                         diss_r(k) = - abs(u_comp(k) - gu) * (                 &
    2181                                     3. * ( u(k,j,i+1) - u(k,j,i) )            &
    2182                                     - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
     2202                        flux_r(k) = ( u_comp(k) - gu ) * (                     &
     2203                                    7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
     2204                                  -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
     2205                                  * adv_mom_3
     2206                        diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
     2207                                    3.0 * ( u(k,j,i+1) - u(k,j,i)   )          &
     2208                                  -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
     2209                                  * adv_mom_3
    21832210                     ENDDO
    21842211                     
    21852212                  CASE ( 7 )
    2186                      DO  k = nzb_u_inner(j,i) + 1, nzt                           
     2213                     DO  k = nzb_u_inner(j,i)+1, nzt                           
    21872214                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2188                         flux_n(k) = v_comp * (                                &
    2189                                     7. * ( u(k,j+1,i) + u(k,j,i) )            &
    2190                                     - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    2191                         diss_n(k) = - abs(v_comp) * (                         &
    2192                                     3. * ( u(k,j+1,i) - u(k,j,i) )            &
    2193                                     - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
     2215                        flux_n(k) = v_comp * (                                 &
     2216                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
     2217                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
     2218                                  * adv_mom_3
     2219                        diss_n(k) = - ABS( v_comp ) * (                        &
     2220                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          &
     2221                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
     2222                                  * adv_mom_3
    21942223                     ENDDO
    21952224                     
    21962225                  CASE ( 8 )
    2197                      DO  k = nzb_u_inner(j,i) + 1, nzt
     2226                     DO  k = nzb_u_inner(j,i)+1, nzt
    21982227                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2199                         flux_n(k) = v_comp * (                                &
    2200                                     7. * ( u(k,j+1,i) + u(k,j,i) )            &
    2201                                     - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    2202                         diss_n(k) = - abs(v_comp) * (                         &
    2203                                     3. * ( u(k,j+1,i) - u(k,j,i) )            &
    2204                                     - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
     2228                        flux_n(k) = v_comp * (                                 &
     2229                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
     2230                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
     2231                                  * adv_mom_3
     2232                        diss_n(k) = - ABS( v_comp ) * (                        &
     2233                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          &
     2234                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
     2235                                  * adv_mom_3
    22052236                     ENDDO 
    22062237                     
     
    22102241!               
    22112242!--            Compute the crosswise 5th order fluxes at the outflow
    2212                IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2     &
    2213               .OR. boundary_flags(j,i) == 5)  THEN
     2243               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
     2244                    boundary_flags(j,i) == 5 )  THEN
    22142245             
    2215                   DO  k = nzb_u_inner(j,i) + 1, nzt
     2246                  DO  k = nzb_u_inner(j,i)+1, nzt
    22162247                     v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2217                      flux_n(k) = v_comp * (                                   &
    2218                                37.  * ( u(k,j+1,i) + u(k,j,i)   )             &
    2219                                - 8. * ( u(k,j+2,i) +u(k,j-1,i)  )             &
    2220                                +      ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    2221                      diss_n(k) = - abs(v_comp) * (                            &
    2222                                10.  * ( u(k,j+1,i) - u(k,j,i) )               &
    2223                                - 5. * ( u(k,j+2,i) - u(k,j-1,i) )             &
    2224                                +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     2248                     flux_n(k) = v_comp * (                                    &
     2249                                 37.0 * ( u(k,j+1,i) + u(k,j,i)   )            &
     2250                               -  8.0 * ( u(k,j+2,i) +u(k,j-1,i)  )            &
     2251                               +        ( u(k,j+3,i) + u(k,j-2,i) ) )          &
     2252                               * adv_mom_5
     2253                     diss_n(k) = - ABS( v_comp ) * (                           &
     2254                                 10.0 * ( u(k,j+1,i) - u(k,j,i)   )            &
     2255                               -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )            &
     2256                               +        ( u(k,j+3,i) - u(k,j-2,i) ) )          &
     2257                               * adv_mom_5
    22252258                  ENDDO
    22262259                 
    22272260               ELSE
    22282261               
    2229                   DO  k = nzb_u_inner(j,i) + 1, nzt
     2262                  DO  k = nzb_u_inner(j,i)+1, nzt
    22302263                     u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2231                      flux_r(k) = ( u_comp(k) - gu ) * (                       &
    2232                                37.  * ( u(k,j,i+1) + u(k,j,i)   )             &
    2233                                - 8. * ( u(k,j,i+2) + u(k,j,i-1) )             &
    2234                                +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    2235                      diss_r(k) = - abs(u_comp(k) - gu) * (                    &
    2236                                10.  * ( u(k,j,i+1) - u(k,j,i)   )             &
    2237                                - 5. * ( u(k,j,i+2) - u(k,j,i-1) )             &
    2238                                +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     2264                     flux_r(k) = ( u_comp(k) - gu ) * (                        &
     2265                                 37.0 * ( u(k,j,i+1) + u(k,j,i)   )            &
     2266                               -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )            &
     2267                               +        ( u(k,j,i+3) + u(k,j,i-2) ) )          &
     2268                               * adv_mom_5
     2269                     diss_r(k) = - ABS( u_comp(k) - gu ) * (                   &
     2270                                 10.0 * ( u(k,j,i+1) - u(k,j,i)   )            &
     2271                               -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )            &
     2272                               +        ( u(k,j,i+3) - u(k,j,i-2) ) )          &
     2273                               * adv_mom_5
    22392274                  ENDDO
    22402275                 
     
    22432278            ELSE
    22442279           
    2245                DO  k = nzb_u_inner(j,i) + 1, nzt
     2280               DO  k = nzb_u_inner(j,i)+1, nzt
    22462281                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2247                   flux_r(k) = ( u_comp(k) - gu ) * (                          &
    2248                             37.  * ( u(k,j,i+1) + u(k,j,i)   )                &
    2249                             - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                &
    2250                             +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    2251                   diss_r(k) = - abs(u_comp(k) - gu) * (                       &
    2252                             10.  * ( u(k,j,i+1) - u(k,j,i)   )                &
    2253                             - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                &
    2254                             +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     2282                  flux_r(k) = ( u_comp(k) - gu ) * (                           &
     2283                              37.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
     2284                            -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )               &
     2285                            +        ( u(k,j,i+3) + u(k,j,i-2) ) )             &
     2286                            * adv_mom_5
     2287                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
     2288                              10.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
     2289                            -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )               &
     2290                            +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    22552291
    22562292                  v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2257                   flux_n(k) = v_comp * (                                      &
    2258                             37.  * ( u(k,j+1,i) + u(k,j,i)   )                &
    2259                             - 8. * ( u(k,j+2,i) + u(k,j-1,i) )                &
    2260                             +      ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    2261                   diss_n(k) = - abs(v_comp) * (                               &
    2262                             10.  * ( u(k,j+1,i) - u(k,j,i)   )                &
    2263                             - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                &
    2264                             +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     2293                  flux_n(k) = v_comp * (                                       &
     2294                              37.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
     2295                            -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )               &
     2296                            +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     2297                  diss_n(k) = - ABS( v_comp ) * (                              &
     2298                              10.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
     2299                            -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )               &
     2300                            +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    22652301                               
    22662302               ENDDO
     
    22682304            ENDIF
    22692305           
    2270             DO  k = nzb_u_inner(j,i) + 1, nzt
     2306            DO  k = nzb_u_inner(j,i)+1, nzt
    22712307           
    2272                tend(k,j,i) = tend(k,j,i) - (                                  &
    2273             ( flux_r(k) + diss_r(k)               &
    2274             - ( swap_flux_x_local_u(k,j) + swap_diss_x_local_u(k,j) ) ) * ddx &
    2275             + ( flux_n(k) + diss_n(k)               &
    2276             - ( swap_flux_y_local_u(k) + swap_diss_y_local_u(k)     ) ) * ddy )
     2308               tend(k,j,i) = tend(k,j,i) - (                                   &
     2309              ( flux_r(k) + diss_r(k)                                          &
     2310            -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j)   ) * ddx &
     2311            + ( flux_n(k) + diss_n(k)                                          &
     2312            -   swap_flux_y_local_u(k) - swap_diss_y_local_u(k)      ) * ddy )
    22772313                 
    22782314               swap_flux_x_local_u(k,j) = flux_r(k)   
     
    22812317               swap_diss_y_local_u(k)   = diss_n(k)
    22822318                     
    2283                sums_us2_ws_l(k,:)  = sums_us2_ws_l(k,:)                       &
    2284                  + ( flux_r(k)                                                &
    2285                  * ( u_comp(k) - 2. * hom(k,1,1,:) )                          &
    2286                  / ( u_comp(k) - gu + 1.0E-20 )                               &
    2287                  + diss_r(k)                                                  &
    2288                  * abs(u_comp(k) - 2. * hom(k,1,1,:) )                        &
    2289                  / (abs(u_comp(k) - gu) + 1.0E-20) )                          &
    2290                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     2319               sums_us2_ws_l(k,:)  = sums_us2_ws_l(k,:)                        &
     2320                 + ( flux_r(k)                                                 &
     2321                 * ( u_comp(k) - 2.0 * hom(k,1,1,:) )                          &
     2322                 / ( u_comp(k) - gu + 1.0E-20 )                                &
     2323                 +   diss_r(k)                                                 &
     2324                 *   ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                     &
     2325                 / ( ABS( u_comp(k) - gu) + 1.0E-20) )                         &
     2326                 *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    22912327            ENDDO
    2292             sums_us2_ws_l(nzb_u_inner(j,i),:) =                               &
     2328            sums_us2_ws_l(nzb_u_inner(j,i),:) =                                &
    22932329                                           sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
    22942330         ENDDO
     
    23012337       DO i = nxlu, nxr
    23022338          DO  j = nys, nyn
    2303              k = nzb_u_inner(j,i) + 1
     2339             k = nzb_u_inner(j,i)+1
    23042340!             
    23052341!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
    23062342!--         reasons the top flux at the surface should be 0.
    2307             flux_t(nzb_u_inner(j,i)) = 0.
    2308             diss_t(nzb_u_inner(j,i)) = 0.
     2343            flux_t(nzb_u_inner(j,i)) = 0.0
     2344            diss_t(nzb_u_inner(j,i)) = 0.0
    23092345            flux_d = flux_t(k-1)
    23102346            diss_d = diss_t(k-1)             
    23112347!
    2312 !--         2nd order scheme
     2348!--         2nd order scheme (bottom)
    23132349             w_comp    = w(k,j,i) + w(k,j,i-1)
    23142350             flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
    2315              diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),          &
    2316                                    0., 0., w_comp, 0.25, ddzw(k) )
     2351             diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),           &
     2352                                   0.0, 0.0, w_comp, 0.25, ddzw(k) )
    23172353                                   
    23182354             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2319                                       - ( flux_d + diss_d )     ) * ddzw(k)
    2320 !
    2321 !--         WS3 as an intermediate step
    2322             k         = nzb_u_inner(j,i) + 2
     2355                                       -   flux_d    -  diss_d       ) * ddzw(k)
     2356!
     2357!--         WS3 as an intermediate step (bottom)
     2358            k         = nzb_u_inner(j,i)+2
    23232359            flux_d    = flux_t(k-1)
    23242360            diss_d    = diss_t(k-1)
    23252361            w_comp    = w(k,j,i) + w(k,j,i-1)
    23262362            flux_t(k) = w_comp * (                                            &
    2327                         7. * (u(k+1,j,i) + u(k,j,i)    )                      &
    2328                         -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    2329             diss_t(k) = - abs(w_comp) * (                                     &
    2330                         3. * ( u(k+1,j,i) - u(k,j,i) )                        &
    2331                         -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
     2363                        7.0 * ( u(k+1,j,i) + u(k,j,i)   )                     &
     2364                      -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
     2365            diss_t(k) = - ABS( w_comp ) * (                                   &
     2366                        3.0 * ( u(k+1,j,i) - u(k,j,i)   )                     &
     2367                      -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    23322368
    23332369            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2334                                       - ( flux_d + diss_d )     ) * ddzw(k)
     2370                                      -   flux_d    -  diss_d      ) * ddzw(k)
    23352371!
    23362372!WS5
    2337              DO  k = nzb_u_inner(j,i) + 3, nzt - 3
     2373             DO  k = nzb_u_inner(j,i)+3, nzt-2
    23382374
    23392375                flux_d    = flux_t(k-1)
     
    23412377                w_comp    = w(k,j,i) + w(k,j,i-1)
    23422378                flux_t(k) = w_comp * (                                        &
    2343                             37.  * ( u(k+1,j,i) + u(k,j,i)   )                &
    2344                             - 8. * ( u(k+2,j,i) + u(k-1,j,i) )                &
    2345                             +      ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
    2346                 diss_t(k) = - abs(w_comp) * (                                 &
    2347                             10.  * ( u(k+1,j,i) - u(k,j,i)   )                &
    2348                             - 5. * ( u(k+2,j,i) - u(k-1,j,i) )                &
    2349                             +      ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
     2379                            37.0 * ( u(k+1,j,i) + u(k,j,i)   )                &
     2380                          -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                &
     2381                                ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
     2382                diss_t(k) = - ABS( w_comp ) * (                               &
     2383                            10.0 * ( u(k+1,j,i) - u(k,j,i)   )                &
     2384                          -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                &
     2385                                ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
    23502386
    23512387                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
    2352                                       - ( flux_d + diss_d )     ) * ddzw(k)
     2388                                          -   flux_d    -  diss_d   ) * ddzw(k)
    23532389
    23542390             ENDDO
    23552391!
    2356 !--         WS3 as an intermediate step
    2357              DO k = nzt - 2, nzt - 1
    2358                 flux_d = flux_t(k-1)
    2359                 diss_d = diss_t(k-1)
    2360                 w_comp = w(k,j,i) + w(k,j,i-1)
    2361                 flux_t(k) = w_comp * (                                        &
    2362                             7. * ( u(k+1,j,i) + u(k,j,i) )                    &
    2363                             -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    2364                 diss_t(k) = - abs(w_comp) * (                                 &
    2365                             3. * ( u(k+1,j,i) - u(k,j,i) )                    &
    2366                             -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
     2392!--          WS3 as an intermediate step (top)
     2393             k         = nzt-1
     2394             flux_d    = flux_t(k-1)
     2395             diss_d    = diss_t(k-1)
     2396             w_comp    = w(k,j,i) + w(k,j,i-1)
     2397             flux_t(k) = w_comp * (                                           &
     2398                         7.0 * ( u(k+1,j,i) + u(k,j,i) )                      &
     2399                       -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
     2400             diss_t(k) = - ABS( w_comp ) * (                                  &
     2401                         3.0 * ( u(k+1,j,i) - u(k,j,i) )                      &
     2402                       -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    23672403                           
    2368                 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
    2369                                       - ( flux_d + diss_d )     ) * ddzw(k)
    2370              ENDDO
    2371 !
    2372 !--         2nd order scheme
     2404             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
     2405                                       -   flux_d    -  diss_d      ) * ddzw(k)
     2406!
     2407!--         2nd order scheme (top)
    23732408             k         = nzt
    23742409             flux_d    = flux_t(k-1)
     
    23812416
    23822417             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
    2383                                       - ( flux_d + diss_d )     ) * ddzw(k)
     2418                                       -   flux_d    -  diss_d      ) * ddzw(k)
    23842419!
    23852420!-- at last vertical momentum flux is accumulated
    23862421            DO  k = nzb_u_inner(j,i), nzt
    23872422               sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:)                       &
    2388                    + ( flux_t(k) + diss_t(k) )                                 &
    2389                    * weight_substep(intermediate_timestep_count)               &
    2390                    * rmask(j,i,:)
     2423                              + ( flux_t(k) + diss_t(k) )                      &
     2424                              *   weight_substep(intermediate_timestep_count)  &
     2425                              *  rmask(j,i,:)
    23912426            ENDDO
    23922427          ENDDO
     
    23972432   
    23982433   
    2399 !
     2434!------------------------------------------------------------------------------!
    24002435! Advection of v - Call for all grid points
    24012436!------------------------------------------------------------------------------!
     
    24292464       
    24302465          IF ( boundary_flags(j,i) == 6 )  THEN
    2431              DO  k = nzb_v_inner(j,i) + 1, nzt
     2466             DO  k = nzb_v_inner(j,i)+1, nzt
    24322467                u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
    24332468                swap_flux_x_local_v(k,j) = u_comp *                           &
     
    24432478             DO  k = nzb_v_inner(j,i)+1, nzt
    24442479                u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
    2445                 swap_flux_x_local_v(k,j) = u_comp * (                         &
    2446                                          37.  * ( v(k,j,i) + v(k,j,i-1)   )   &
    2447                                          - 8. * ( v(k,j,i+1) + v(k,j,i-2) )   &
    2448                                          +      ( v(k,j,i+2) + v(k,j,i-3) ) ) &
     2480                swap_flux_x_local_v(k,j) = u_comp * (                          &
     2481                                           37.0 * ( v(k,j,i) + v(k,j,i-1)   )  &
     2482                                         -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )  &
     2483                                         +        ( v(k,j,i+2) + v(k,j,i-3) ) )&
    24492484                                         * adv_mom_5
    2450                 swap_diss_x_local_v(k,j) = - abs(u_comp) * (                  &
    2451                                          10. * ( v(k,j,i) - v(k,j,i-1)   )    &
    2452                                          -5. * ( v(k,j,i+1) - v(k,j,i-2) )    &
    2453                                          +     ( v(k,j,i+2) - v(k,j,i-3) ) )  &
     2485                swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                 &
     2486                                           10.0 * ( v(k,j,i) - v(k,j,i-1)   )  &
     2487                                         -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )  &
     2488                                         +        ( v(k,j,i+2) - v(k,j,i-3) ) )&
    24542489                                         * adv_mom_5
    24552490             ENDDO
     
    24632498         IF ( boundary_flags(j,i) == 7 )  THEN
    24642499         
    2465             DO  k = nzb_v_inner(j,i) + 1, nzt
     2500            DO  k = nzb_v_inner(j,i)+1, nzt
    24662501               v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
    24672502               swap_flux_y_local_v(k) = v_comp(k) *                           &
     
    24752510         ELSE
    24762511         
    2477             DO  k = nzb_v_inner(j,i) + 1, nzt
     2512            DO  k = nzb_v_inner(j,i)+1, nzt
    24782513               v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
    2479                swap_flux_y_local_v(k) = v_comp(k) * (                         &
    2480                                       37.  * ( v(k,j,i) + v(k,j-1,i)   )      &
    2481                                       - 8. * ( v(k,j+1,i) + v(k,j-2,i) )      &
    2482                                       +      ( v(k,j+2,i) + v(k,j-3,i) ) )    &
     2514               swap_flux_y_local_v(k) = v_comp(k) * (                          &
     2515                                        37.0 * ( v(k,j,i) + v(k,j-1,i)   )     &
     2516                                      -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )     &
     2517                                      +        ( v(k,j+2,i) + v(k,j-3,i) ) )   &
    24832518                                      * adv_mom_5
    2484                swap_diss_y_local_v(k) = - abs(v_comp(k)) * (                  &
    2485                                       10.  * ( v(k,j,i) - v(k,j-1,i)   )      &
    2486                                       - 5. * ( v(k,j+1,i) - v(k,j-2,i) )      &
    2487                                       +      ( v(k,j+2,i) - v(k,j-3,i) ) )    &
     2519               swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                 &
     2520                                        10.0 * ( v(k,j,i) - v(k,j-1,i)   )     &
     2521                                      -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )     &
     2522                                      +        ( v(k,j+2,i) - v(k,j-3,i) ) )   &
    24882523                                      * adv_mom_5
    24892524            ENDDO
     
    24982533         
    24992534                  CASE ( 1 )
    2500                      DO  k = nzb_v_inner(j,i) + 1, nzt
     2535                     DO  k = nzb_v_inner(j,i)+1, nzt
    25012536                        u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    2502                         flux_r(k) = u_comp * (                                &
    2503                                  7. * (v(k,j,i+1) + v(k,j,i)    )             &
    2504                                  -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    2505                         diss_r(k) = - abs(u_comp) * (                         &
    2506                                3. * ( v(k,j,i+1) - v(k,j,i)   )               &
    2507                                -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
     2537                        flux_r(k) = u_comp * (                                 &
     2538                                    7.0 * (v(k,j,i+1) + v(k,j,i)    )          &
     2539                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )        &
     2540                                  * adv_mom_3
     2541                        diss_r(k) = - ABS( u_comp ) * (                        &
     2542                                    3.0 * ( v(k,j,i+1) - v(k,j,i)   )          &
     2543                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )        &
     2544                                  * adv_mom_3
    25082545                     ENDDO
    25092546               
    25102547                  CASE ( 2 )
    2511                      DO  k = nzb_v_inner(j,i) + 1, nzt
     2548                     DO  k = nzb_v_inner(j,i)+1, nzt
    25122549                        u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    25132550                        flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25
     
    25182555               
    25192556                  CASE ( 3 )
    2520                      DO  k = nzb_v_inner(j,i) + 1, nzt
     2557                     DO  k = nzb_v_inner(j,i)+1, nzt
    25212558                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
    25222559                        flux_n(k) = ( v_comp(k)- gv ) * (                     &
    2523                                  7. * ( v(k,j+1,i) + v(k,j,i)   )             &
    2524                                  -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    2525                         diss_n(k) = - abs(v_comp(k) - gv) * (                 &
    2526                                  3. * ( v(k,j+1,i) - v(k,j,i)   )             &
    2527                                  -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
     2560                                    7.0 * ( v(k,j+1,i) + v(k,j,i)   )         &
     2561                                  -       ( v(k,j+2,i) + v(k,j-1,i) ) )       &
     2562                                  * adv_mom_3
     2563                        diss_n(k) = - ABS(v_comp(k) - gv) * (                 &
     2564                                    3.0 * ( v(k,j+1,i) - v(k,j,i)   )         &
     2565                                  -       ( v(k,j+2,i) - v(k,j-1,i) ) )       &
     2566                                  * adv_mom_3
    25282567                     ENDDO
    25292568               
    25302569                  CASE ( 4 )
    2531                      DO  k = nzb_v_inner(j,i) + 1, nzt
     2570                     DO  k = nzb_v_inner(j,i)+1, nzt
    25322571                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
    25332572                        flux_n(k) = ( v_comp(k) - gv ) *                      &
     
    25392578               
    25402579                  CASE ( 5 )
    2541                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2580                     DO  k = nzb_w_inner(j,i)+1, nzt
    25422581                        u_comp    = u(k,j-1,i) + u(k,j,i) - gu
    25432582                        flux_r(k) = u_comp * (                                &
    2544                                7. * ( v(k,j,i+1) + v(k,j,i)   )               &
    2545                                -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    2546                         diss_r(k) = - abs(u_comp) * (                         &
    2547                                3. * ( w(k,j,i+1) - w(k,j,i)   )               &
    2548                                -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
     2583                                    7.0 * ( v(k,j,i+1) + v(k,j,i)   )         &
     2584                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )       &
     2585                                  * adv_mom_3
     2586                        diss_r(k) = - ABS (u_comp ) * (                       &
     2587                                    3.0 * ( w(k,j,i+1) - w(k,j,i)   )         &
     2588                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )       &
     2589                                  * adv_mom_3
    25492590                     ENDDO
    25502591                 
    25512592                  CASE ( 6 )
    2552                      DO  k = nzb_v_inner(j,i) + 1, nzt
     2593                     DO  k = nzb_v_inner(j,i)+1, nzt
    25532594                                 
    25542595                        u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
    25552596                        flux_r(k) = u_comp * (                                &
    2556                                  7. * ( v(k,j,i+1) + v(k,j,i)   )             &
    2557                                  -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    2558                         diss_r(k) = - abs(u_comp) * (                         &
    2559                                  3. * ( v(k,j,i+1) - v(k,j,i)   )             &
    2560                                  -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
     2597                                    7.0 * ( v(k,j,i+1) + v(k,j,i)   )         &
     2598                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )       &
     2599                                  * adv_mom_3
     2600                        diss_r(k) = - ABS( u_comp ) * (                       &
     2601                                    3.0 * ( v(k,j,i+1) - v(k,j,i)   )         &
     2602                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )       &
     2603                                  * adv_mom_3
    25612604                     ENDDO
    25622605               
    25632606                  CASE ( 7 )
    2564                      DO  k = nzb_v_inner(j,i) + 1, nzt                                 
     2607                     DO  k = nzb_v_inner(j,i)+1, nzt                                 
    25652608                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
    25662609                        flux_n(k) = ( v_comp(k) - gv ) * (                    &
    2567                                  7. * ( v(k,j+1,i) + v(k,j,i)   )             &
    2568                                  -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    2569                         diss_n(k) = - abs(v_comp(k) - gv) * (                 &
    2570                                  3. * ( v(k,j+1,i) - v(k,j,i)   )             &
    2571                                  -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
     2610                                    7.0 * ( v(k,j+1,i) + v(k,j,i)   )         &
     2611                                  -       ( v(k,j+2,i) + v(k,j-1,i) ) )       &
     2612                                  * adv_mom_3
     2613                        diss_n(k) = - ABS( v_comp(k) - gv ) * (               &
     2614                                    3.0 * ( v(k,j+1,i) - v(k,j,i)   )         &
     2615                                  -       ( v(k,j+2,i) - v(k,j-1,i) ) )       &
     2616                                  * adv_mom_3
    25722617                     ENDDO
    25732618               
     
    25762621               END SELECT
    25772622!         
    2578 !--         Compute the crosswise 5th order fluxes at the outflow
    2579                IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2    &
    2580             .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
     2623!--            Compute the crosswise 5th order fluxes at the outflow
     2624               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
     2625                    boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )   &
     2626               THEN
    25812627         
    2582                   DO  k = nzb_v_inner(j,i) + 1, nzt
     2628                  DO  k = nzb_v_inner(j,i)+1, nzt
    25832629                     v_comp(k) = v(k,j+1,i) + v(k,j,i)
    25842630                     flux_n(k) = ( v_comp(k) - gv ) * (                       &
    2585                                37.  * ( v(k,j+1,i) + v(k,j,i)   )             &
    2586                                - 8. * ( v(k,j+2,i) + v(k,j-1,i) )             &
    2587                                +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    2588                      diss_n(k) = - abs(v_comp(k) - gv ) * (                   &
    2589                                10.  * ( v(k,j+1,i) - v(k,j,i)   )             &
    2590                                - 5. * ( v(k,j+2,i) - v(k,j-1,i) )             &
    2591                                +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
     2631                                 37.0 * ( v(k,j+1,i) + v(k,j,i)   )           &
     2632                               -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )           &
     2633                               +        ( v(k,j+3,i) + v(k,j-2,i) ) )         &
     2634                               * adv_mom_5
     2635                     diss_n(k) = - ABS( v_comp(k) - gv ) * (                  &
     2636                                 10.0 * ( v(k,j+1,i) - v(k,j,i)   )           &
     2637                               -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )           &
     2638                               +        ( v(k,j+3,i) - v(k,j-2,i) ) )         &
     2639                               * adv_mom_5
    25922640                   ENDDO
    25932641             
    25942642               ELSE
    25952643           
    2596                   DO  k = nzb_v_inner(j,i) + 1, nzt
     2644                  DO  k = nzb_v_inner(j,i)+1, nzt
    25972645                     u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    25982646                     flux_r(k) = u_comp * (                                   &
    2599                                37.  * ( v(k,j,i+1) + v(k,j,i)   )             &
    2600                               - 8. * ( v(k,j,i+2) + v(k,j,i-1) )              &
    2601                             +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    2602                      diss_r(k) = - abs(u_comp) * (                            &
    2603                                10. * ( v(k,j,i+1) - v(k,j,i)   )              &
    2604                                -5. * ( v(k,j,i+2) - v(k,j,i-1) )              &
    2605                                +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     2647                                 37.0 * ( v(k,j,i+1) + v(k,j,i)   )           &
     2648                              -   8.0 * ( v(k,j,i+2) + v(k,j,i-1) )           &
     2649                              +         ( v(k,j,i+3) + v(k,j,i-2) ) )         &
     2650                              * adv_mom_5
     2651                     diss_r(k) = - ABS( u_comp ) * (                          &
     2652                                 10.0 * ( v(k,j,i+1) - v(k,j,i)   )           &
     2653                               -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )           &
     2654                               +        ( v(k,j,i+3) - v(k,j,i-2) ) )         &
     2655                               * adv_mom_5
    26062656                  ENDDO
    26072657               
     
    26112661            ELSE
    26122662         
    2613                DO  k = nzb_v_inner(j,i) + 1, nzt
     2663               DO  k = nzb_v_inner(j,i)+1, nzt
    26142664                  u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    26152665                  flux_r(k) = u_comp * (                                      &
    2616                             37.  * ( v(k,j,i+1) + v(k,j,i)   )                &
    2617                             - 8. * ( v(k,j,i+2) + v(k,j,i-1) )                &
    2618                             +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    2619                   diss_r(k) = - abs(u_comp) * (                               &
    2620                             10. * ( v(k,j,i+1) - v(k,j,i)   )                 &
    2621                             -5. * ( v(k,j,i+2) - v(k,j,i-1) )                 &
    2622                             +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     2666                              37.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
     2667                            -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )              &
     2668                            +        ( v(k,j,i+3) + v(k,j,i-2) ) )            &
     2669                            * adv_mom_5
     2670                  diss_r(k) = - ABS ( u_comp ) * (                            &
     2671                              10.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
     2672                            -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )              &
     2673                            +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    26232674
    26242675                  v_comp(k) = v(k,j+1,i) + v(k,j,i)
    26252676                  flux_n(k) = ( v_comp(k) - gv ) * (                          &
    2626                             37.  * ( v(k,j+1,i) + v(k,j,i)   )                &
    2627                             - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                &
    2628                             +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    2629                   diss_n(k) = - abs(v_comp(k) - gv ) * (                      &
    2630                             10.  * ( v(k,j+1,i) - v(k,j,i)   )                &
    2631                             - 5. * ( v(k,j+2,i) - v(k,j-1,i) )                &
    2632                             +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
     2677                              37.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
     2678                            -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )              &
     2679                            +        ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
     2680                  diss_n(k) = - ABS( v_comp(k) - gv ) * (                     &
     2681                              10.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
     2682                            -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )              &
     2683                            +        ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
    26332684
    26342685               ENDDO
    26352686            ENDIF
    26362687           
    2637             DO  k = nzb_v_inner(j,i) + 1, nzt
     2688            DO  k = nzb_v_inner(j,i)+1, nzt
    26382689               tend(k,j,i) = tend(k,j,i) - (                                  &
    26392690              ( flux_r(k) + diss_r(k)                                         &
    2640             - ( swap_flux_x_local_v(k,j) + swap_diss_x_local_v(k,j) ) ) * ddx &
     2691            -   swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j)  ) * ddx &
    26412692            + ( flux_n(k) + diss_n(k)                                         &
    2642             - ( swap_flux_y_local_v(k) + swap_diss_y_local_v(k) )     ) * ddy )
     2693            -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)       ) * ddy )
    26432694                 
    26442695               swap_flux_x_local_v(k,j) = flux_r(k)   
     
    26472698               swap_diss_y_local_v(k)   = diss_n(k)   
    26482699
    2649                sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                        &
    2650                   + (flux_n(k) *( v_comp(k) - 2. * hom(k,1,2,:) )             &
    2651                   / ( v_comp(k) - gv + 1.0E-20)                               &
    2652                   + diss_n(k) * abs(v_comp(k) - 2. * hom(k,1,2,:) )           &
    2653                   /(abs(v_comp(k) - gv) + 1.0E-20 ) )                         &
    2654                   * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     2700               sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                         &
     2701                  + ( flux_n(k) * ( v_comp(k) - 2.0 * hom(k,1,2,:) )           &
     2702                  / ( v_comp(k) - gv + 1.0E-20 )                               &
     2703                  +   diss_n(k) * ABS( v_comp(k) - 2.0 * hom(k,1,2,:) )        &
     2704                  / ( ABS( v_comp(k) - gv ) + 1.0E-20 ) )                      &
     2705                  *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    26552706            ENDDO
    2656             sums_vs2_ws_l(nzb_v_inner(j,i),:) =                             
     2707            sums_vs2_ws_l(nzb_v_inner(j,i),:) =                               
    26572708                                           sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
    26582709         ENDDO
     
    26672718!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
    26682719!--         reasons the top flux at the surface should be 0.
    2669              flux_t(nzb_v_inner(j,i)) = 0.
    2670              diss_t(nzb_v_inner(j,i)) = 0.
    2671              k      = nzb_v_inner(j,i) + 1             
     2720             flux_t(nzb_v_inner(j,i)) = 0.0
     2721             diss_t(nzb_v_inner(j,i)) = 0.0
     2722             k      = nzb_v_inner(j,i)+1             
    26722723             flux_d = flux_t(k-1)
    26732724             diss_d = diss_t(k-1)
    26742725!
    2675 !--          2nd order scheme
     2726!--          2nd order scheme (bottom)
    26762727             w_comp    = w(k,j-1,i) + w(k,j,i)
    26772728             flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
    26782729             diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i),          &
    2679                                    0., 0., w_comp, 0.25, ddzw(k) )
     2730                                   0.0, 0.0, w_comp, 0.25, ddzw(k) )
    26802731
    26812732             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
    2682                                       - ( flux_d + diss_d )     ) * ddzw(k)
    2683 !
    2684 !--         WS3 as an intermediate step
     2733                                       -   flux_d    - diss_d       ) * ddzw(k)
     2734!
     2735!--          WS3 as an intermediate step (bottom)
    26852736             k          = nzb_v_inner(j,i) + 2
    26862737             flux_d     = flux_t(k-1)
     
    26882739             w_comp     = w(k,j-1,i) + w(k,j,i)
    26892740             flux_t(k) = w_comp * (                                           &
    2690                          7. * ( v(k+1,j,i) + v(k,j,i) )                       &
    2691                          -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
    2692              diss_t(k) = - abs(w_comp) * (                                    &
    2693                          3. * ( v(k+1,j,i) - v(k,j,i) )                       &
    2694                          -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
     2741                         7.0 * ( v(k+1,j,i) + v(k,j,i) )                      &
     2742                       -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
     2743             diss_t(k) = - ABS( w_comp ) * (                                  &
     2744                         3.0 * ( v(k+1,j,i) - v(k,j,i) )                      &
     2745                       -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
    26952746                         
    26962747             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
    2697                                       - ( flux_d + diss_d )     ) * ddzw(k)
     2748                                       -   flux_d    - diss_d       ) * ddzw(k)
    26982749
    26992750!
    27002751!-- WS5
    2701              DO  k = nzb_v_inner(j,i) + 3, nzt - 3
     2752             DO  k = nzb_v_inner(j,i)+3, nzt-2
    27022753                flux_d = flux_t(k-1)
    27032754                diss_d = diss_t(k-1)
    27042755                w_comp = w(k,j-1,i) + w(k,j,i)
    27052756                flux_t(k) = w_comp * (                                        &
    2706                             37.  * ( v(k+1,j,i) + v(k,j,i)   )                &
    2707                             - 8. * ( v(k+2,j,i) + v(k-1,j,i) )                &
    2708                             +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
    2709                 diss_t(k) = - abs(w_comp) * (                                 &
    2710                             10. * ( v(k+1,j,i) - v(k,j,i)   )                 &
    2711                             -5. * ( v(k+2,j,i) - v(k-1,j,i) )                 &
    2712                             +     ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
     2757                            37.0 * ( v(k+1,j,i) + v(k,j,i)   )                &
     2758                          -  8.0 * ( v(k+2,j,i) + v(k-1,j,i) )                &
     2759                                ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
     2760                diss_t(k) = - ABS( w_comp ) * (                               &
     2761                            10.0 * ( v(k+1,j,i) - v(k,j,i)   )                &
     2762                          -  5.0 * ( v(k+2,j,i) - v(k-1,j,i) )                &
     2763                          +        ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
    27132764                           
    2714                 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
    2715                                       - ( flux_d + diss_d )     ) * ddzw(k)
     2765                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
     2766                                          -   flux_d    - diss_d    ) * ddzw(k)
    27162767            ENDDO
    27172768!
    2718 !--         WS3 as an intermediate step
    2719             DO k = nzt - 2, nzt - 1
    2720                flux_d     = flux_t(k-1)
    2721                diss_d     = diss_t(k-1)
    2722                w_comp     = w(k,j-1,i) + w(k,j,i)
    2723                flux_t(k) = w_comp * (                                         &
    2724                            7. * ( v(k+1,j,i) + v(k,j,i) )                     &
    2725                            -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
    2726                diss_t(k) = - abs(w_comp) * (                                  &
    2727                            3. * ( v(k+1,j,i) - v(k,j,i) )                     &
    2728                            -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
     2769!--         WS3 as an intermediate step (top)
     2770            k     = nzt-1
     2771            flux_d    = flux_t(k-1)
     2772            diss_d    = diss_t(k-1)
     2773            w_comp    = w(k,j-1,i) + w(k,j,i)
     2774            flux_t(k) = w_comp * (                                            &
     2775                        7.0 * ( v(k+1,j,i) + v(k,j,i)  )                     &
     2776                      -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
     2777            diss_t(k) = - ABS( w_comp ) * (                                   &
     2778                        3.0 * ( v(k+1,j,i) - v(k,j,i)  )                     &
     2779                      -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
    27292780                         
    2730              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
    2731                                       - ( flux_d + diss_d )     ) * ddzw(k)
    2732             ENDDO
    2733 !
    2734 !--          2nd order scheme
     2781            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
     2782                                      -   flux_d    - diss_d     ) * ddzw(k)
     2783!
     2784!--         2nd order scheme (top)
    27352785            k         = nzt
    27362786            flux_d    = flux_t(k-1)
     
    27382788            w_comp    = w(k,j-1,i) + w(k,j,i)
    27392789            flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
    2740             diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i),       &
    2741                                   v(k-1,j,i), v(k-2,j,i), w_comp, 0.25,ddzw(k))
    2742 
    2743             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2744                                       - ( flux_d + diss_d )     ) * ddzw(k)
    2745 !
    2746 !- at last vertical momentum flux is accumulated
     2790            diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i),        &
     2791                                  v(k-1,j,i), v(k-2,j,i), w_comp, 0.25, ddzw(k))
     2792
     2793            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                &
     2794                                      -   flux_d    - diss_d       ) * ddzw(k)
     2795!
     2796!-          At last vertical momentum flux is accumulated.
    27472797            DO  k = nzb_v_inner(j,i), nzt
    27482798               sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:)                       &
    2749                    + ( flux_t(k) + diss_t(k) )                                 &
    2750                    * weight_substep(intermediate_timestep_count)              &
    2751                    * rmask(j,i,:)
     2799                               + ( flux_t(k) + diss_t(k) )                     &
     2800                               *   weight_substep(intermediate_timestep_count) &
     2801                               *  rmask(j,i,:)
    27522802            ENDDO
    27532803            sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                &
    2754                    sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
     2804                                             sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
    27552805          ENDDO
    27562806       ENDDO
     
    27592809   
    27602810   
    2761 !
     2811!------------------------------------------------------------------------------!
    27622812! Advection of w - Call for all grid points
    27632813!------------------------------------------------------------------------------!
     
    27902840          IF ( boundary_flags(j,i) == 6 )  THEN
    27912841         
    2792              DO  k = nzb_v_inner(j,i) + 1, nzt
     2842             DO  k = nzb_v_inner(j,i)+1, nzt
    27932843                u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
    27942844                swap_flux_x_local_w(k,j) = u_comp *                           &
     
    28022852          ELSE
    28032853         
    2804              DO  k = nzb_v_inner(j,i) + 1, nzt
     2854             DO  k = nzb_v_inner(j,i)+1, nzt
    28052855                u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
    2806                 swap_flux_x_local_w(k,j) = u_comp * (                         &
    2807                                          37.  * ( w(k,j,i) + w(k,j,i-1)   )   &
    2808                                          - 8. * ( w(k,j,i+1) + w(k,j,i-2) )   &
    2809                                          +      ( w(k,j,i+2) + w(k,j,i-3) ) ) &
     2856                swap_flux_x_local_w(k,j) = u_comp * (                          &
     2857                                           37.0 * ( w(k,j,i) + w(k,j,i-1)   )  &
     2858                                         -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )  &
     2859                                         +        ( w(k,j,i+2) + w(k,j,i-3) ) )&
    28102860                                         * adv_mom_5
    2811                 swap_diss_x_local_w(k,j) = - abs(u_comp) * (                  &
    2812                                          10.  * ( w(k,j,i) - w(k,j,i-1)   )   &
    2813                                          - 5. * ( w(k,j,i+1) - w(k,j,i-2) )   &
    2814                                          +      ( w(k,j,i+2) - w(k,j,i-3) ) ) &
     2861                swap_diss_x_local_w(k,j) = - ABS( u_comp ) * (                 &
     2862                                           10.0 * ( w(k,j,i) - w(k,j,i-1)   )  &
     2863                                         -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )  &
     2864                                         +        ( w(k,j,i+2) - w(k,j,i-3) ) )&
    28152865                                         * adv_mom_5
    28162866             ENDDO
     
    28242874         IF ( boundary_flags(j,i) == 8 )  THEN
    28252875         
    2826             DO  k = nzb_v_inner(j,i) + 1, nzt
     2876            DO  k = nzb_v_inner(j,i)+1, nzt
    28272877               v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
    2828                swap_flux_y_local_w(k) = v_comp *                              &
     2878               swap_flux_y_local_w(k) = v_comp *                               &
    28292879                                        ( w(k,j,i) + w(k,j-1,i) ) * 0.25
    2830                swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i),     &
    2831                                                   w(k,j,i), w(k,j-1,i),       &
    2832                                                   w(k,j-1,i), v_comp, 0.25,ddy)
     2880               swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i),      &
     2881                                                  w(k,j,i), w(k,j-1,i),        &
     2882                                                  w(k,j-1,i), v_comp, 0.25, ddy)
    28332883            ENDDO
    28342884           
    28352885         ELSE
    28362886         
    2837             DO  k = nzb_v_inner(j,i) + 1, nzt           
     2887            DO  k = nzb_v_inner(j,i)+1, nzt           
    28382888               v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
    2839                swap_flux_y_local_w(k) = v_comp * (                            &
    2840                                       37. *  ( w(k,j,i) + w(k,j-1,i)  )       &
    2841                                       - 8. * ( w(k,j+1,i) +w(k,j-2,i) )       &
    2842                                       +      ( w(k,j+2,i) +w(k,j-3,i) ) )     &
     2889               swap_flux_y_local_w(k) = v_comp * (                             &
     2890                                        37.0 * ( w(k,j,i) + w(k,j-1,i)  )      &
     2891                                      -  8.0 * ( w(k,j+1,i) +w(k,j-2,i) )      &
     2892                                      +        ( w(k,j+2,i) +w(k,j-3,i) ) )    &
    28432893                                      * adv_mom_5
    2844                swap_diss_y_local_w(k) = - abs(v_comp) * (                     &
    2845                                       10.  * ( w(k,j,i) - w(k,j-1,i)   )      &
    2846                                       - 5. * ( w(k,j+1,i) - w(k,j-2,i) )      &
    2847                                       +      ( w(k,j+2,i) - w(k,j-3,i) ) )    &
     2894               swap_diss_y_local_w(k) = - ABS( v_comp ) * (                    &
     2895                                        10.0 * ( w(k,j,i) - w(k,j-1,i)   )     &
     2896                                      -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )     &
     2897                                      +        ( w(k,j+2,i) - w(k,j-3,i) ) )   &
    28482898                                      * adv_mom_5
    28492899            ENDDO
     
    28592909         
    28602910                  CASE ( 1 )
    2861                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2911                     DO  k = nzb_w_inner(j,i)+1, nzt
    28622912                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    28632913                        flux_r(k) = u_comp * (                                &
    2864                               7. * ( w(k,j,i+1) + w(k,j,i)   )                &
    2865                               -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
    2866                         diss_r(k) = -abs(u_comp) * (                          &
    2867                               3. * ( w(k,j,i+1) - w(k,j,i)   )                &
    2868                               -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
     2914                                    7.0 * ( w(k,j,i+1) + w(k,j,i)   )         &
     2915                                  -       ( w(k,j,i+2) + w(k,j,i-1) ) )       &
     2916                                  * adv_mom_3
     2917                        diss_r(k) = -ABS( u_comp ) * (                        &
     2918                                    3.0 * ( w(k,j,i+1) - w(k,j,i)   )         &
     2919                                  -       ( w(k,j,i+2) - w(k,j,i-1) ) )       &
     2920                                  * adv_mom_3 
    28692921                     ENDDO
    28702922               
    28712923                  CASE ( 2 )
    2872                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2924                     DO  k = nzb_w_inner(j,i)+1, nzt
    28732925                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    28742926                        flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25
     
    28792931               
    28802932                  CASE ( 3 )
    2881                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2933                     DO  k = nzb_w_inner(j,i)+1, nzt
    28822934                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    28832935                        flux_n(k) = v_comp * (                                &
    2884                               7. * ( w(k,j+1,i) + w(k,j,i)   )                &
    2885                               -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    2886                         diss_n(k) = -abs(v_comp) * (                          &
    2887                               3. * ( w(k,j+1,i) - w(k,j,i)   )                &
    2888                               -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
     2936                                     7.0 * ( w(k,j+1,i) + w(k,j,i)   )        &
     2937                                  -        ( w(k,j+2,i) + w(k,j-1,i) ) )      &
     2938                                  * adv_mom_3
     2939                        diss_n(k) = -ABS( v_comp ) * (                        &
     2940                                    3.0 * ( w(k,j+1,i) - w(k,j,i)   )         &
     2941                                  -       ( w(k,j+2,i) - w(k,j-1,i) ) )       &
     2942                                  * adv_mom_3
    28892943                     ENDDO
    28902944               
    28912945                  CASE ( 4 )
    2892                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2946                     DO  k = nzb_w_inner(j,i)+1, nzt
    28932947                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    28942948                        flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25
     
    28992953               
    29002954                  CASE ( 5 )
    2901                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2955                     DO  k = nzb_w_inner(j,i)+1, nzt
    29022956                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    29032957                        flux_r(k) = u_comp * (                                &
    2904                               7. * ( w(k,j,i+1) + w(k,j,i)   )                &
    2905                               -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
    2906                         diss_r(k) = - abs(u_comp) * (                         &
    2907                               3. * ( w(k,j,i+1) - w(k,j,i) )                  &
    2908                               -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
     2958                                    7.0 * ( w(k,j,i+1) + w(k,j,i)   )         &
     2959                                  -       ( w(k,j,i+2) + w(k,j,i-1) ) )       &
     2960                                  * adv_mom_3
     2961                        diss_r(k) = - ABS( u_comp ) * (                       &
     2962                                    3.0 * ( w(k,j,i+1) - w(k,j,i) )           &
     2963                                  -       ( w(k,j,i+2) - w(k,j,i-1) ) )       &
     2964                                  * adv_mom_3
    29092965                     ENDDO
    29102966               
    29112967                  CASE ( 6 )
    2912                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2968                     DO  k = nzb_w_inner(j,i)+1, nzt
    29132969                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    29142970                        flux_r(k) = u_comp *(                                 &
    2915                               7. * ( w(k,j,i+1) + w(k,j,i)   )                &
    2916                               -    ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
    2917                         diss_r(k) = - abs(u_comp) * (                         &
    2918                               3. * ( w(k,j,i+1) - w(k,j,i) )                  &
    2919                               -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
     2971                                    7.0 * ( w(k,j,i+1) + w(k,j,i)   )         &
     2972                                  -       ( w(k,j,i+2) + w(k,j,i-1) ) )       &
     2973                                  * adv_mom_3
     2974                        diss_r(k) = - ABS( u_comp ) * (                       &
     2975                                    3.0 * ( w(k,j,i+1) - w(k,j,i) )           &
     2976                                  -       ( w(k,j,i+2) - w(k,j,i-1) ) )       &
     2977                                  * adv_mom_3
    29202978                     ENDDO
    29212979               
    29222980                  CASE ( 7 )
    2923                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2981                     DO  k = nzb_w_inner(j,i)+1, nzt
    29242982                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    29252983                        flux_n(k) = v_comp *(                                 &
    2926                               7. * ( w(k,j+1,i) + w(k,j,i)   )                &
    2927                               -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    2928                         diss_n(k) = - abs(v_comp) * (                         &
    2929                               3. * ( w(k,j+1,i) - w(k,j,i)   )                &
    2930                               -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
     2984                                    7.0 * ( w(k,j+1,i) + w(k,j,i)   )         &
     2985                                  -       ( w(k,j+2,i) + w(k,j-1,i) ) )       &
     2986                                  * adv_mom_3
     2987                        diss_n(k) = - ABS( v_comp ) * (                       &
     2988                                    3.0 * ( w(k,j+1,i) - w(k,j,i)   )         &
     2989                                  -       ( w(k,j+2,i) - w(k,j-1,i) ) )       &
     2990                                  * adv_mom_3
    29312991                      ENDDO
    29322992               
    29332993                  CASE ( 8 )
    2934                      DO  k = nzb_w_inner(j,i) + 1, nzt
     2994                     DO  k = nzb_w_inner(j,i)+1, nzt
    29352995                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    29362996                        flux_n(k) = v_comp * (                                &
    2937                               7. * ( w(k,j+1,i) + w(k,j,i)   )                &
    2938                              -     ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    2939                         diss_n(k) = - abs(v_comp) * (                         &
    2940                               3. * ( w(k,j+1,i) - w(k,j,i) )                  &
    2941                               -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
     2997                                    7.0 * ( w(k,j+1,i) + w(k,j,i)   )         &
     2998                                  -       ( w(k,j+2,i) + w(k,j-1,i) ) )       &
     2999                                  * adv_mom_3
     3000                        diss_n(k) = - ABS( v_comp ) * (                       &
     3001                                    3.0 * ( w(k,j+1,i) - w(k,j,i) )           &
     3002                                  -       ( w(k,j+2,i) - w(k,j-1,i) ) )       &
     3003                                  * adv_mom_3
    29423004               
    29433005                     ENDDO 
     
    29473009               END SELECT
    29483010!         
    2949 !--      Compute the crosswise 5th order fluxes at the outflow
    2950                IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2    &
    2951             .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
     3011!--            Compute the crosswise 5th order fluxes at the outflow
     3012               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
     3013                    boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )   &
     3014               THEN
    29523015         
    2953                   DO  k = nzb_w_inner(j,i) + 1, nzt
     3016                  DO  k = nzb_w_inner(j,i)+1, nzt
    29543017                     v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    29553018                     flux_n(k) = v_comp * (                                   &
    2956                            37.  * ( w(k,j+1,i) + w(k,j,i)   )                 &
    2957                            - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                 &
    2958                            +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    2959                      diss_n(k) = - abs(v_comp) * (                            &
    2960                            10.  * ( w(k,j+1,i) - w(k,j,i)   )                 &
    2961                            - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                 &
    2962                            +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     3019                                 37.0 * ( w(k,j+1,i) + w(k,j,i)   )           &
     3020                               -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )           &
     3021                               +        ( w(k,j+3,i) + w(k,j-2,i) ) )         &
     3022                               * adv_mom_5
     3023                     diss_n(k) = - ABS( v_comp ) * (                          &
     3024                                 10.0 * ( w(k,j+1,i) - w(k,j,i)   )           &
     3025                               -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )           &
     3026                               +        ( w(k,j+3,i) - w(k,j-2,i) ) )         &
     3027                               * adv_mom_5
    29633028                  ENDDO
    29643029           
    29653030               ELSE
    29663031         
    2967                   DO  k = nzb_w_inner(j,i) + 1, nzt         
     3032                  DO  k = nzb_w_inner(j,i)+1, nzt         
    29683033                     u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    29693034                     flux_r(k) = u_comp * (                                   &
    2970                            37.  * ( w(k,j,i+1) + w(k,j,i) )                   &
    2971                            - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                 &
    2972                            +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    2973                      diss_r(k) = - abs(u_comp) * (                            &
    2974                            10.  * ( w(k,j,i+1) - w(k,j,i)   )                 &
    2975                            - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                 &
    2976                             +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     3035                                 37.0 * ( w(k,j,i+1) + w(k,j,i) )             &
     3036                               -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )           &
     3037                               +        ( w(k,j,i+3) + w(k,j,i-2) ) )         &
     3038                               * adv_mom_5
     3039                     diss_r(k) = - ABS( u_comp ) * (                          &
     3040                                 10.0 * ( w(k,j,i+1) - w(k,j,i)   )           &
     3041                               -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )           &
     3042                               +        ( w(k,j,i+3) - w(k,j,i-2) ) )         &
     3043                               * adv_mom_5
    29773044                  ENDDO
    29783045           
     
    29823049!       
    29833050!--            Compute the fifth order fluxes for the interior of PE domain.               
    2984                DO  k = nzb_w_inner(j,i) + 1, nzt
     3051               DO  k = nzb_w_inner(j,i)+1, nzt
    29853052                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    29863053                  flux_r(k) = u_comp * (                                      &
    2987                         37.  * ( w(k,j,i+1) + w(k,j,i)   )                    &
    2988                         - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                    &
    2989                         +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    2990                   diss_r(k) = - abs(u_comp) * (                               &
    2991                         10.  * ( w(k,j,i+1) - w(k,j,i)   )                    &
    2992                         - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                    &
    2993                         +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     3054                              37.0 * ( w(k,j,i+1) + w(k,j,i)   )              &
     3055                            -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )              &
     3056                            +        ( w(k,j,i+3) + w(k,j,i-2) ) )            &
     3057                            * adv_mom_5
     3058                  diss_r(k) = - ABS( u_comp ) * (                             &
     3059                              10.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
     3060                            -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )              &
     3061                            +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    29943062                 
    29953063                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    29963064                  flux_n(k) = v_comp * (                                      &
    2997                         37.  * ( w(k,j+1,i) + w(k,j,i)   )                    &
    2998                         - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                    &
    2999                         +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    3000                   diss_n(k) = - abs(v_comp) * (                               &
    3001                         10.  * ( w(k,j+1,i) - w(k,j,i)   )                    &
    3002                         - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                    &
    3003                         +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     3065                              37.0 * ( w(k,j+1,i) + w(k,j,i)   )              &
     3066                            -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )              &
     3067                            +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
     3068                  diss_n(k) = - ABS( v_comp ) * (                             &
     3069                              10.0 * ( w(k,j+1,i) - w(k,j,i)   )              &
     3070                            -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )              &
     3071                            +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    30043072               ENDDO
    30053073         
     
    30093077               tend(k,j,i) = tend(k,j,i) - (                                  &
    30103078              ( flux_r(k) + diss_r(k)                                         &
    3011             - ( swap_flux_x_local_w(k,j) + swap_diss_x_local_w(k,j) ) ) * ddx &
     3079            -   swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j)  ) * ddx &
    30123080            + ( flux_n(k) + diss_n(k)                                         &
    3013             - ( swap_flux_y_local_w(k) + swap_diss_y_local_w(k)     ) ) * ddy )
     3081            -   swap_flux_y_local_w(k) - swap_diss_y_local_w(k)      ) * ddy )
    30143082               
    30153083               swap_flux_x_local_w(k,j) = flux_r(k)
     
    30273095       DO i = nxl, nxr
    30283096          DO  j = nys, nyn
    3029             k      = nzb_w_inner(j,i) + 1
     3097            k      = nzb_w_inner(j,i)+1
    30303098            flux_d = w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i))
    30313099            flux_t(k-1) = flux_d
     
    30333101            diss_t(k-1) = diss_d
    30343102!
    3035 !--         2nd order scheme            
     3103!--         2nd order scheme (bottom)           
    30363104            w_comp = w(k+1,j,i) + w(k,j,i)
    30373105            flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
    3038             diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0.,  &
     3106            diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0, &
    30393107                                  w_comp, 0.25, ddzu(k+1) )
    30403108
    30413109            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    3042                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
     3110                                      -   flux_d    - diss_d      ) * ddzu(k+1)
    30433111!                                 
    3044 !--        WS3 as an intermediate step
     3112!--        WS3 as an intermediate step (bottom)
    30453113            k         = nzb_w_inner(j,i) + 2
     3114            flux_d    = flux_t(k-1)
     3115            diss_d    = diss_t(k-1)
     3116            w_comp    = w(k+1,j,i) + w(k,j,i)
     3117            flux_t(k) = w_comp * (                                            &
     3118                        7.0 * ( w(k+1,j,i) + w(k,j,i) )                       &
     3119                      -       ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
     3120            diss_t(k) = - ABS( w_comp ) * (                                   &
     3121                        3.0 * ( w(k+1,j,i) - w(k,j,i) )                       &
     3122                      -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
     3123                       
     3124            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
     3125                                      -   flux_d    - diss_d      ) * ddzu(k+1)
     3126!
     3127!--         WS5
     3128            DO  k = nzb_v_inner(j,i)+3, nzt-2
     3129               flux_d = flux_t(k-1)
     3130               diss_d = diss_t(k-1)
     3131
     3132               w_comp = w(k+1,j,i) + w(k,j,i)
     3133               flux_t(k) = w_comp * (                                          &
     3134                           37.0 * ( w(k+1,j,i) + w(k,j,i)   )                  &
     3135                         -  8.0 * ( w(k+2,j,i) + w(k-1,j,i) )                  &
     3136                         +        ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
     3137               diss_t(k) = - ABS( w_comp ) * (                                 &
     3138                           10.0 * ( w(k+1,j,i) - w(k,j,i)   )                  &
     3139                         -  5.0 * ( w(k+2,j,i) - w(k-1,j,i) )                  &
     3140                         +        ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
     3141
     3142               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)             &
     3143                                         -   flux_d    - diss_d    ) * ddzu(k+1)
     3144            ENDDO
     3145!                                 
     3146!--         WS3 as an intermediate step (top)
     3147            k         = nzt-1
    30463148            flux_d    = flux_t(k-1)
    30473149            diss_d    = diss_t(k-1)
    30483150            w_comp    = w(k+1,j,i)+w(k,j,i)
    30493151            flux_t(k) = w_comp * (                                            &
    3050                         7. * ( w(k+1,j,i) + w(k,j,i) )                        &
    3051                         - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
    3052             diss_t(k) = - abs(w_comp) * (                                     &
    3053                         3. * ( w(k+1,j,i) - w(k,j,i) )                        &
    3054                         - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
     3152                        7.0 * ( w(k+1,j,i) + w(k,j,i)   )                     &
     3153                      -      ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
     3154            diss_t(k) = - ABS( w_comp ) * (                                   &
     3155                          3.0 * ( w(k+1,j,i) - w(k,j,i)   )                   &
     3156                      -      ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
    30553157                       
    30563158            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    3057                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
    3058 !
    3059 !--         WS5
    3060             DO  k = nzb_v_inner(j,i) + 3, nzt
    3061                flux_d = flux_t(k-1)
    3062                diss_d = diss_t(k-1)
    3063 
    3064                w_comp = w(k+1,j,i) + w(k,j,i)
    3065                flux_t(k) = w_comp * (                                         &
    3066                            37.  * ( w(k+1,j,i) + w(k,j,i)   )                 &
    3067                            - 8. * ( w(k+2,j,i) + w(k-1,j,i) )                 &
    3068                            +      ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
    3069                diss_t(k) = - abs(w_comp) * (                                  &
    3070                            10.  * ( w(k+1,j,i) - w(k,j,i)   )                 &
    3071                            - 5. * ( w(k+2,j,i) - w(k-1,j,i) )                 &
    3072                            +      ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
    3073 
    3074                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
    3075                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
    3076             ENDDO
    3077 !                                 
    3078 !--        WS3 as an intermediate step
    3079             DO k = nzt - 2, nzt - 1
    3080                flux_d    = flux_t(k-1)
    3081                diss_d    = diss_t(k-1)
    3082                w_comp    = w(k+1,j,i)+w(k,j,i)
    3083                flux_t(k) = w_comp * (                                         &
    3084                            7. * ( w(k+1,j,i) + w(k,j,i) )                     &
    3085                            - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
    3086                diss_t(k) = - abs(w_comp) * (                                  &
    3087                            3. * ( w(k+1,j,i) - w(k,j,i) )                     &
    3088                            - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
    3089                        
    3090                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
    3091                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
    3092             ENDDO
    3093 !
    3094 !--         2nd order scheme 
     3159                                      -   flux_d    - diss_d      ) * ddzu(k+1)
     3160!
     3161!--         2nd order scheme (top)
    30953162            k = nzt
    30963163            flux_d = flux_t(k-1)
     
    31033170
    31043171            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    3105                                  - ( flux_d + diss_d )     ) * ddzu(k+1)
     3172                                      -   flux_d    - diss_d      ) * ddzu(k+1)
    31063173!
    31073174!--         at last vertical momentum flux is accumulated
    31083175            DO  k = nzb_w_inner(j,i), nzt
    31093176               sums_ws2_ws_l(k,:)  = sums_ws2_ws_l(k,:)                        &
    3110                    + ( flux_t(k) + diss_t(k) )                                 &
    3111                    * weight_substep(intermediate_timestep_count)              &
    3112                    * rmask(j,i,:)
     3177                               + ( flux_t(k) + diss_t(k) )                     &
     3178                               *   weight_substep(intermediate_timestep_count) &
     3179                               *  rmask(j,i,:)
    31133180            ENDDO
    31143181
     
    31313198!   
    31323199!-- Computation of 2nd order dissipation. This numerical dissipation is
    3133 !-- necessary to keep a stable numerical solution in region where the
    3134 !-- order of the schemes degraded.
     3200!-- necessary to keep a stable numerical solution in regions where the
     3201!-- order of the schemes is degraded.
    31353202 
    31363203     REAL FUNCTION diss_2nd( v2, v1, v0, vm1, vm2, vel_comp, factor, grid )   &
    3137                             RESULT(dissip)
     3204                            RESULT( dissip )
    31383205       
    31393206       IMPLICIT NONE
     
    31443211               value_min_p, value_max_p, diss_m, diss_0, diss_p
    31453212             
    3146        value_min_m = MIN(v0, vm1, vm2)
    3147        value_max_m = MAX(v0, vm1, vm2)
    3148        value_min = MIN(v1, v0, vm2)
    3149        value_max = MAX(v1, v0, vm2)
    3150        value_min_p = MIN(v2, v1, v0)
    3151        value_max_p = MAX(v2, v1, v0)
     3213       value_min_m = MIN( v0, vm1, vm2 )
     3214       value_max_m = MAX( v0, vm1, vm2 )
     3215       value_min   = MIN( v1, v0,  vm2 )
     3216       value_max   = MAX( v1, v0,  vm2 )
     3217       value_min_p = MIN( v2, v1,  v0  )
     3218       value_max_p = MAX( v2, v1,  v0  )
    31523219       
    3153        diss_m = MAX(0.,vm1 - value_max_m) + MIN(0.,vm1 - value_min_m)
    3154        diss_0 = MAX(0.,v0 - value_max) + MIN(0.,v0 - value_min)
    3155        diss_p = MAX(0.,v1 - value_max_p) + MIN(0.,v1 - value_min_p)
     3220       diss_m = MAX( 0.0, vm1 - value_max_m ) + MIN( 0.0, vm1 - value_min_m )
     3221       diss_0 = MAX( 0.0, v0  - value_max   ) + MIN( 0.0, v0  - value_min   )
     3222       diss_p = MAX( 0.0, v1  - value_max_p ) + MIN( 0.0, v1  - value_min_p )
    31563223       
    3157        dissip = abs(vel_comp) * (diss_p - 2.*diss_0 + diss_m)         &
    3158                          * grid**2 * factor
     3224       dissip = ABS( vel_comp ) * ( diss_p - 2.0 * diss_0 + diss_m )          &
     3225                                * grid**2 * factor
    31593226             
    31603227    END FUNCTION diss_2nd
  • palm/trunk/SOURCE/init_3d_model.f90

    r710 r713  
    77! Current revisions:
    88! -----------------
    9 !
     9! ! Reformulate weight_substep and weight_pres as broken numbers.
    1010!
    1111! Former revisions:
     
    13531353    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
    13541354
    1355        weight_substep(1) = 0.166666666666666
    1356        weight_substep(2) = 0.3
    1357        weight_substep(3) = 0.533333333333333
    1358 
    1359        weight_pres(1)    = 0.333333333333333
    1360        weight_pres(2)    = 0.416666666666666
    1361        weight_pres(3)    = 0.25
     1355       weight_substep(1) = 1./6.
     1356       weight_substep(2) = 3./10.
     1357       weight_substep(3) = 8./15.
     1358
     1359       weight_pres(1)    = 1./3.
     1360       weight_pres(2)    = 5./12.
     1361       weight_pres(3)    = 1./4.
    13621362
    13631363    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
    13641364
    1365        weight_substep(1) = 0.5
    1366        weight_substep(2) = 0.5
     1365       weight_substep(1) = 1./2.
     1366       weight_substep(2) = 1./2.
    13671367         
    1368        weight_pres(1)    = 0.5
    1369        weight_pres(2)    = 0.5         
     1368       weight_pres(1)    = 1./2.
     1369       weight_pres(2)    = 1./2.       
    13701370
    13711371    ELSE                                     ! for Euler- and leapfrog-method
Note: See TracChangeset for help on using the changeset viewer.