Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Mar 26, 2012 2:18:34 PM (12 years ago)
Author:
suehring
Message:

WS5 is available in combination with topography. Version number changed from 3.8 to 3.8a. Modification in init_pt_anomaly.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/advec_ws.f90

    r857 r861  
    44! Current revisions:
    55! ------------------
     6! ws-scheme also work with topography in combination with vector version.
     7! ws-scheme also work with outflow boundaries in combination with
     8! vector version.
     9
     10! Degradation of the applied order of scheme is now steered by multiplying with
     11! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each
     12! grid point and mulitplied with the appropriate flag.
     13! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
     14! term derived according to the 4th and 6th order terms is applied. It turns
     15! out that diss_2nd does not provide sufficient dissipation near walls.
     16! Therefore, the function diss_2nd is removed.
     17
     18! Near walls a divergence correction is necessary to overcome numerical
     19! instabilities due to too less divergence reduction of the velocity field.
     20!
     21! boundary_flags and logicals steering the degradation are removed.
     22! Empty SUBROUTINE local_diss removed.
    623!
    724! Former revisions:
    825! -----------------
    926! $Id$
    10 !
    11 ! 856 2012-03-20 18:29:17Z suehring
    12 ! Bug concerning numerical dissipation at the first grid level. Function
    13 ! diss_2nd always returns zero.  Because diss_2nd will be replaced by numerical
    14 ! dissipation which is more consistent to the WS-schemes, the return value
    15 ! is set to zero for now.
    1627!
    1728! 801 2012-01-10 17:30:36Z suehring
    … …  
    6172! routine using for initialisation and steering of the statical evaluation.
    6273! The computation of turbulent fluxes takes place inside the advection
    63 ! routines.
    64 ! In case of vector architectures Dirichlet and Radiation boundary conditions
    65 ! are outstanding and not available.
    66 ! A further routine local_diss_ij is available (next weeks) and is used if a
    67 ! control of dissipative fluxes is desired.
    68 !
    69 ! OUTSTANDING: - Dirichlet and Radiation boundary conditions for
    70 !                vector architectures
    71 !              - dissipation control for cache architectures ( next weeks )
    72 !              - Topography ( next weeks )
     74! routines.
     75! Near non-cyclic boundaries the order of the applied advection scheme is
     76! degraded.
     77! A divergence correction is applied. It is necessary for topography, since
     78! the divergence is not sufficiently reduced, resulting in erroneous fluxes and
     79! partly numerical instabilities.
    7380!-----------------------------------------------------------------------------!
    7481
    7582    PRIVATE
    7683    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
    77              local_diss, ws_init, ws_statistics
     84             ws_init, ws_statistics
    7885
    7986    INTERFACE ws_init
    … …  
    104111       MODULE PROCEDURE advec_w_ws_ij
    105112    END INTERFACE advec_w_ws
    106 
    107     INTERFACE local_diss
    108        MODULE PROCEDURE local_diss
    109        MODULE PROCEDURE local_diss_ij
    110     END INTERFACE local_diss
    111113
    112114 CONTAINS
    … …  
    126128
    127129!
    128 !--    Allocate arrays needed for dissipation control.
    129        IF ( dissipation_control )  THEN
    130 !          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),        &
    131 !                   var_y(nzb+1:nzt,nys:nyn,nxl:nxr),        &
    132 !                   var_z(nzb+1:nzt,nys:nyn,nxl:nxr),        &
    133 !                   gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    134 !                   gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    135 !                   gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
    136        ENDIF
    137 
    138 !
    139130!--    Set the appropriate factors for scalar and momentum advection.
    140131       adv_sca_5 = 1./60.
    141132       adv_sca_3 = 1./12.
     133       adv_sca_1 = 1./2.
    142134       adv_mom_5 = 1./120.
    143135       adv_mom_3 = 1./24.
    144 
     136       adv_mom_1 = 1./4.
    145137!         
    146138!--    Arrays needed for statical evaluation of fluxes.
    … …  
    231223       ENDIF
    232224
    233 !
    234 !--    Determine the flags where the order of the scheme for horizontal
    235 !--    advection has to be degraded.
    236        ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
    237        DO  i = nxl, nxr
    238           DO  j = nys, nyn
    239 
    240              boundary_flags(j,i) = 0
    241              IF ( outflow_l )  THEN
    242                 IF ( i == nxlu  )  boundary_flags(j,i) = 5
    243                 IF ( i == nxl   )  boundary_flags(j,i) = 6
    244              ELSEIF ( outflow_r )  THEN
    245                 IF ( i == nxr-1 )  boundary_flags(j,i) = 1
    246                 IF ( i == nxr   )  boundary_flags(j,i) = 2
    247              ELSEIF ( outflow_n )  THEN
    248                 IF ( j == nyn-1 )  boundary_flags(j,i) = 3
    249                 IF ( j == nyn   )  boundary_flags(j,i) = 4
    250              ELSEIF ( outflow_s )  THEN
    251                 IF ( j == nysv  )  boundary_flags(j,i) = 7
    252                 IF ( j == nys   )  boundary_flags(j,i) = 8
    253              ENDIF
    254 
    255           ENDDO
    256        ENDDO
    257        
    258225    END SUBROUTINE ws_init
    259226
    … …  
    271238!       
    272239!--    The arrays needed for statistical evaluation are set to to 0 at the
    273 !--    begin of prognostic_equations.
     240!--    beginning of prognostic_equations.
    274241       IF ( ws_scheme_mom )  THEN
    275242          sums_wsus_ws_l = 0.0
    … …  
    307274       IMPLICIT NONE
    308275
    309        INTEGER ::  i, i_omp, j, k, tn
    310        LOGICAL ::  degraded_l, degraded_s
    311        REAL    ::  flux_d, diss_d, u_comp, v_comp
     276       INTEGER ::  i, i_omp, j, k, tn, k_ppp, k_pp, k_mm
     277       REAL    ::  flux_d, diss_d, u_comp, v_comp, div
    312278       REAL, DIMENSION(:,:,:), POINTER    :: sk
    313279       REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
    … …  
    320286       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
    321287
    322 
    323        degraded_l = .FALSE.
    324        degraded_s = .FALSE.
    325 
    326        IF ( boundary_flags(j,i) /= 0 )  THEN
    327 !       
    328 !--       Degrade the order for Dirichlet bc. at the outflow boundary
    329           SELECT CASE ( boundary_flags(j,i) )
    330 
    331              CASE ( 1 )
    332 
    333                 DO  k = nzb_s_inner(j,i)+1, nzt
    334                    u_comp    = u(k,j,i+1) - u_gtrans
    335                    flux_r(k) = u_comp * (                                      &
    336                                7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
    337                                    - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    338                    diss_r(k) = -ABS( u_comp ) * (                              &
    339                                3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
    340                                    - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    341                 ENDDO
    342 
    343              CASE ( 2 )
    344 
    345                 DO  k = nzb_s_inner(j,i)+1, nzt
    346                    u_comp    = u(k,j,i+1) - u_gtrans
    347                    flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
    348                    diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
    349                                          sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
    350                                          0.5, ddx )
    351                 ENDDO
    352 
    353              CASE ( 3 )
    354 
    355                 DO  k = nzb_s_inner(j,i)+1, nzt
    356                    v_comp = v(k,j+1,i) - v_gtrans
    357                    flux_n(k) = v_comp * (                                      &
    358                                7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
    359                                    - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    360                    diss_n(k) = -ABS( v_comp ) * (                              &
    361                                3.0 * ( sk(k,j+1,i) - sk(k,j,i)     )           &
    362                                    - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    363                ENDDO
    364 
    365              CASE ( 4 )
    366 
    367                 DO  k = nzb_s_inner(j,i)+1, nzt
    368                    v_comp    = v(k,j+1,i) - v_gtrans
    369                    flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
    370                    diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
    371                                          sk(k,j-1,i), sk(k,j-2,i), v_comp,    &
    372                                          0.5, ddy )     
    373                 ENDDO
    374 
    375              CASE ( 5 )
    376 
    377                 DO  k = nzb_w_inner(j,i)+1, nzt
    378                    u_comp    = u(k,j,i+1) - u_gtrans
    379                    flux_r(k) = u_comp  * (                                     &
    380                                7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
    381                                    - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    382                    diss_r(k) = -ABS( u_comp ) * (                              &
    383                                3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
    384                                    - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    385                 ENDDO 
    386 
    387              CASE ( 6 )
    388 
    389                 DO  k = nzb_s_inner(j,i)+1, nzt
    390                    u_comp    = u(k,j,i+1) - u_gtrans
    391                    flux_r(k) = u_comp * (                                      &
    392                                7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
    393                                    - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    394                    diss_r(k) = -ABS( u_comp ) * (                              &
    395                                3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
    396                                    - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    397 !                             
    398 !--                Compute leftside fluxes for the left boundary of PE domain
    399                    u_comp                 = u(k,j,i) - u_gtrans
    400                    swap_flux_x_local(k,j,tn) = u_comp * (                         &
    401                                             sk(k,j,i) + sk(k,j,i-1) ) * 0.5
    402                    swap_diss_x_local(k,j,tn) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
    403                                                       sk(k,j,i), sk(k,j,i-1),  &
    404                                                       sk(k,j,i-1), u_comp,     &
    405                                                       0.5, ddx )
    406                 ENDDO
    407                 degraded_l = .TRUE.
    408 
    409              CASE ( 7 )
    410 
    411                 DO  k = nzb_s_inner(j,i)+1, nzt
    412                    v_comp    = v(k,j+1,i)-v_gtrans
    413                    flux_n(k) = v_comp * (                                      &
    414                                7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
    415                                    - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    416                    diss_n(k) = -ABS( v_comp ) * (                              &
    417                                3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
    418                                    - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    419                 ENDDO
    420 
    421              CASE ( 8 )
    422 
    423                 DO  k = nzb_s_inner(j,i)+1, nzt
    424                    v_comp    = v(k,j+1,i) - v_gtrans
    425                    flux_n(k) = v_comp * (                                      &
    426                                7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
    427                                    - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    428                    diss_n(k) = -ABS( v_comp ) * (                              &
    429                                3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
    430                                    - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    431 !                             
    432 !--                Compute southside fluxes for the south boundary of PE domain
    433                    v_comp               = v(k,j,i) - v_gtrans
    434                    swap_flux_y_local(k,tn) = v_comp *                             &
    435                                           ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
    436                    swap_diss_y_local(k,tn) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
    437                                                     sk(k,j,i), sk(k,j-1,i),    &
    438                                                     sk(k,j-1,i), v_comp,       &
    439                                                     0.5, ddy )
    440                 ENDDO
    441                 degraded_s = .TRUE.
    442                
    443              CASE DEFAULT
    444            
    445           END SELECT
    446 
    447 !         
    448 !--       Compute the crosswise 5th order fluxes at the outflow
    449           IF ( boundary_flags(j,i) == 1  .OR.  boundary_flags(j,i) == 2  .OR. &
    450                boundary_flags(j,i) == 5  .OR.  boundary_flags(j,i) == 6 )  THEN
    451 
    452              DO  k = nzb_s_inner(j,i)+1, nzt
    453                 v_comp    = v(k,j+1,i) - v_gtrans
    454                 flux_n(k) = v_comp * (                                         &
    455                             37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )               &
    456                           -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
    457                           +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    458                 diss_n(k) = -ABS( v_comp ) * (                                 &
    459                             10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )               &
    460                           -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
    461                           +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    462              ENDDO
    463            
    464           ELSE
    465          
    466              DO  k = nzb_s_inner(j,i)+1, nzt
    467                 u_comp    = u(k,j,i+1) - u_gtrans 
    468                 flux_r(k) = u_comp * (                                         &
    469                             37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )               &
    470                           -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
    471                           +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    472                 diss_r(k) = -ABS( u_comp ) * (                                 &
    473                             10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )               &
    474                           -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
    475                           +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    476              ENDDO
    477            
    478           ENDIF
    479 
    480        ELSE
    481                
    482 !
    483 !--       Compute the fifth order fluxes for the interior of PE domain.
    484           DO  k = nzb_u_inner(j,i)+1, nzt
    485              u_comp    = u(k,j,i+1) - u_gtrans
    486              flux_r(k) = u_comp * (                                           &
    487                          37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
    488                        -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
    489                        +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    490              diss_r(k) = -ABS( u_comp ) * (                                   &
    491                          10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
    492                        -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
    493                        +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    494 
    495              v_comp    = v(k,j+1,i) - v_gtrans
    496              flux_n(k) = v_comp * (                                           &
    497                          37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
    498                        -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
    499                        +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    500              diss_n(k) = -ABS( v_comp ) * (                                   &
    501                          10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
    502                        -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
    503                        +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    504           ENDDO
    505          
     288!
     289!--    Compute southside fluxes of the respective PE bounds.
     290       IF ( j == nys )  THEN
     291!
     292!--       Up to the top of the highest topography.
     293          DO  k = nzb+1, nzb_max
     294
     295             v_comp                  = v(k,j,i) - v_gtrans
     296             swap_flux_y_local(k,tn) = v_comp * (                             &
     297                         ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     298                      +     7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     299                      +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1  &
     300                         ) *                                                  &
     301                                     ( sk(k,j,i)  + sk(k,j-1,i)     )         &
     302                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     303                      +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     304                         ) *                                                  &
     305                                     ( sk(k,j+1,i) + sk(k,j-2,i)    )         &
     306                  +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     307                         ) *                                                  &
     308                                     ( sk(k,j+2,i) + sk(k,j-3,i)    )         &
     309                                                )
     310
     311             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
     312                         ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     313                      +     3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     314                      +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1  &
     315                         ) *                                                  &
     316                                     ( sk(k,j,i)   - sk(k,j-1,i)    )         &
     317                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     318                      +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     319                         ) *                                                  &
     320                                     ( sk(k,j+1,i) - sk(k,j-2,i)  )           &
     321                  +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     322                         ) *                                                  &
     323                                     ( sk(k,j+2,i) - sk(k,j-3,i) )            &
     324                                                        )
     325
     326          ENDDO
     327!
     328!--       Above to the top of the highest topography. No degradation necessary.
     329          DO  k = nzb_max+1, nzt
     330
     331             v_comp               = v(k,j,i) - v_gtrans
     332             swap_flux_y_local(k,tn) = v_comp * (                            &
     333                                    37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
     334                                  -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
     335                                  +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
     336                                             ) * adv_sca_5
     337              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                   &
     338                                    10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
     339                                  -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
     340                                  +          sk(k,j+2,i) - sk(k,j-3,i)       &
     341                                                        ) * adv_sca_5
     342
     343          ENDDO
     344
    506345       ENDIF
    507 !       
    508 !--    Compute left- and southside fluxes of the respective PE bounds.       
    509        IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
    510        
    511            DO  k = nzb_s_inner(j,i)+1, nzt
    512               v_comp               = v(k,j,i) - v_gtrans
    513               swap_flux_y_local(k,tn) = v_comp * (                               &
    514                                      37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
    515                                    -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
    516                                    +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
    517                                               ) * adv_sca_5
    518               swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                       &
    519                                      10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
    520                                    -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
    521                                    +          sk(k,j+2,i) - sk(k,j-3,i)       &
    522                                                       ) * adv_sca_5
    523            ENDDO           
    524          
    525         ENDIF
     346!
     347!--    Compute leftside fluxes of the respective PE bounds.
     348       IF ( i == i_omp )  THEN
    526349       
    527         IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
    528        
    529            DO  k = nzb_s_inner(j,i)+1, nzt
    530               u_comp                 = u(k,j,i) - u_gtrans
    531               swap_flux_x_local(k,j,tn) = u_comp * (                             &
    532                                        37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
    533                                      -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
    534                                      +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
    535                                                 ) * adv_sca_5
    536               swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                     &
    537                                        10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
    538                                      -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
    539                                      +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
    540                                                         ) * adv_sca_5
    541            ENDDO
     350          DO  k = nzb+1, nzb_max
     351
     352             u_comp                 = u(k,j,i) - u_gtrans
     353             swap_flux_x_local(k,j,tn) = u_comp * (                           &
     354                         ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     355                      +     7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     356                      +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1  &
     357                         ) *                                                  &
     358                                     ( sk(k,j,i)   + sk(k,j,i-1)    )         &
     359                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     360                      +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     361                         ) *                                                  &
     362                                     ( sk(k,j,i+1) + sk(k,j,i-2)    )         &
     363                  +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     364                         ) *         ( sk(k,j,i+2) + sk(k,j,i-3)    )         &
     365                                                  )
     366
     367              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
     368                         ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     369                      +     3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     370                      +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1  &
     371                         ) *                                                  &
     372                                     ( sk(k,j,i)   - sk(k,j,i-1)    )         &
     373                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     374                      +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     375                         ) *                                                  &
     376                                     ( sk(k,j,i+1) - sk(k,j,i-2)  )           &
     377                  +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     378                         ) *                                                  &
     379                                     ( sk(k,j,i+2) - sk(k,j,i-3) )            &
     380                                                           )
     381
     382          ENDDO
     383
     384          DO  k = nzb_max+1, nzt
     385
     386             u_comp                 = u(k,j,i) - u_gtrans
     387             swap_flux_x_local(k,j,tn) = u_comp * (                          &
     388                                      37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
     389                                    -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
     390                                    +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
     391                                                  ) * adv_sca_5
     392
     393             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
     394                                      10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
     395                                    -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
     396                                    +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
     397                                                          ) * adv_sca_5
     398
     399          ENDDO
    542400           
    543         ENDIF
    544 
     401       ENDIF
     402
     403       flux_t(0) = 0.0
     404       diss_t(0) = 0.0
     405       flux_d    = 0.0
     406       diss_d    = 0.0
    545407!       
    546 !--    Now compute the tendency terms for the horizontal parts
    547        DO  k = nzb_s_inner(j,i)+1, nzt
    548          
    549           tend(k,j,i) = tend(k,j,i) - (                                      &
    550                           ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
    551                             swap_diss_x_local(k,j,tn) ) * ddx                   &
    552                         + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    553                             swap_diss_y_local(k,tn)   ) * ddy                   &
    554                                       )
     408!--    Now compute the fluxes and tendency terms for the horizontal and
     409!--    vertical parts up to the top of the highest topography.
     410       DO  k = nzb+1, nzb_max
     411!
     412!--       Note: It is faster to conduct all multiplications explicitly, e.g.
     413!--       * adv_sca_5 ... than to determine a factor and multiplicate the
     414!--       flux at the end. Indeed, per loop the first requires 6
     415!--       multiplications more, but the the second requires 9 more internal
     416!--       calls of IBITS().
     417          u_comp    = u(k,j,i+1) - u_gtrans
     418          flux_r(k) = u_comp * (                                            &
     419                     ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     420                  +     7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     421                  +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1    &
     422                     ) *                                                    &
     423                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
     424              -      (  8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     425                  +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     426                     ) *                                                    &
     427                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
     428              +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     429                     ) *     ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
     430                               )
     431
     432          diss_r(k) = -ABS( u_comp ) * (                                    &
     433                     ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     434                  +     3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     435                  +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1    &
     436                     ) *                                                    &
     437                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
     438              -      (  5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     439                  +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     440                     ) *                                                    &
     441                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
     442              +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     443                     ) *                                                    &
     444                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
     445                                       )
     446
     447          v_comp    = v(k,j+1,i) - v_gtrans
     448          flux_n(k) = v_comp * (                                            &
     449                     ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     450                  +     7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     451                  +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1    &
     452                     ) *                                                    &
     453                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
     454              -      (  8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     455                  +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     456                     ) *                                                    &
     457                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
     458              +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     459                     ) *                                                    &
     460                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
     461                               )
     462
     463          diss_n(k) = -ABS( v_comp ) * (                                    &
     464                     ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     465                  +     3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     466                  +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1    &
     467                     ) *                                                    &
     468                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
     469              -      (  5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     470                  +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     471                     ) *                                                    &
     472                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
     473              +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     474                     ) *                                                    &
     475                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
     476                                       )
     477!
     478!--       k index has to be modified near bottom and top, else array
     479!--       subscripts will be exceeded.
     480          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
     481          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1)  )
     482          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
     483
     484
     485          flux_t(k) = w(k,j,i) * (                                          &
     486                     ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     487                  +     7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     488                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     489                     ) *                                                    &
     490                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
     491              -      (  8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     492                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     493                     ) *                                                    &
     494                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
     495              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     496                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
     497                                 )
     498
     499          diss_t(k) = -ABS( w(k,j,i) ) * (                                  &
     500                     ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     501                  +     3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     502                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     503                     ) *                                                    &
     504                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
     505              -      (  5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     506                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     507                     ) *                                                    &
     508                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
     509              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     510                     ) *                                                    &
     511                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
     512                                         )
     513!
     514!--       Calculate the divergence of the velocity field. A respective
     515!--       correction is needed to overcome numerical instabilities caused
     516!--       by a not sufficient reduction of divergences near topography.
     517          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
     518                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
     519                        + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     520
     521          tend(k,j,i) = tend(k,j,i) - (                                       &
     522                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
     523                          swap_diss_x_local(k,j,tn)            ) * ddx        &
     524                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
     525                          swap_diss_y_local(k,tn)              ) * ddy        &
     526                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
     527                                                               ) * ddzw(k)    &
     528                                      ) + sk(k,j,i) * div
    555529
    556530          swap_flux_y_local(k,tn)   = flux_n(k)
    … …  
    558532          swap_flux_x_local(k,j,tn) = flux_r(k)
    559533          swap_diss_x_local(k,j,tn) = diss_r(k)
    560          
     534          flux_d                    = flux_t(k)
     535          diss_d                    = diss_t(k)
     536
    561537       ENDDO
    562 
    563 !
    564 !--    Vertical advection, degradation of order near bottom and top.
    565 !--    The fluxes flux_d and diss_d at the surface are 0. Due to later
    566 !--    calculation of statistics the top flux at the surface should be 0.
    567        flux_t(nzb_s_inner(j,i)) = 0.0
    568        diss_t(nzb_s_inner(j,i)) = 0.0
    569 
    570 !       
    571 !--    2nd-order scheme (bottom)
    572        k         = nzb_s_inner(j,i)+1
    573        flux_d    = flux_t(k-1)
    574        diss_d    = diss_t(k-1)
    575        flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
    576 
    577 !       
    578 !--    sk(k,j,i) is referenced three times to avoid an access below surface
    579        diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i),  &
    580                              sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
    581                    
    582        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
    583                                    * ddzw(k)
    584 !
    585 !--    WS3 as an intermediate step (bottom)
    586        k         = nzb_s_inner(j,i) + 2
    587        flux_d    = flux_t(k-1)
    588        diss_d    = diss_t(k-1)
    589        flux_t(k) = w(k,j,i) * (                                               &
    590                                 7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )           &
    591                                     - ( sk(k+2,j,i) + sk(k-1,j,i) )           &
    592                               ) * adv_sca_3
    593        diss_t(k) = -ABS( w(k,j,i) ) * (                                       &
    594                                 3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )           &
    595                                     - ( sk(k+2,j,i) - sk(k-1,j,i) )           &
    596                                       ) * adv_sca_3
    597 
    598        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
    599                                    * ddzw(k)
    600 !
    601 !--    WS5
    602        DO  k = nzb_s_inner(j,i)+3, nzt-2
    603        
    604           flux_d    = flux_t(k-1)
    605           diss_d    = diss_t(k-1)         
    606           flux_t(k) = w(k,j,i) * (                                            &
    607                                    37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )       &
    608                                  -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )       &
    609                                  +        ( sk(k+3,j,i) + sk(k-2,j,i) )       &
    610                                  ) * adv_sca_5
    611           diss_t(k) = -ABS( w(k,j,i) ) * (                                     &
    612                                            10.0 * ( sk(k+1,j,i) - sk(k,j,i)   )&
    613                                          -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )&
    614                                          +        ( sk(k+3,j,i) - sk(k-2,j,i) )&
    615                                          ) * adv_sca_5
    616 
    617           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    618                                     - ( flux_d + diss_d ) ) * ddzw(k)
     538!
     539!--    Now compute the fluxes and tendency terms for the horizontal and
     540!--    vertical parts above the top of the highest topography. No degradation
     541!--    for the horizontal parts, but for the vertical it is stell needed.
     542       DO  k = nzb_max+1, nzt
     543
     544          u_comp    = u(k,j,i+1) - u_gtrans
     545          flux_r(k) = u_comp * (                                            &
     546                      37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                  &
     547                    -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
     548                    +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     549          diss_r(k) = -ABS( u_comp ) * (                                    &
     550                      10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                  &
     551                    -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
     552                    +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     553
     554          v_comp    = v(k,j+1,i) - v_gtrans
     555          flux_n(k) = v_comp * (                                            &
     556                      37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                  &
     557                    -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
     558                    +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     559          diss_n(k) = -ABS( v_comp ) * (                                    &
     560                      10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                  &
     561                    -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                  &
     562                    +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     563!
     564!--       k index has to be modified near bottom and top, else array
     565!--       subscripts will be exceeded.
     566          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
     567          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1)  )
     568          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
     569
     570          flux_t(k) = w(k,j,i) * (                                          &
     571                     ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     572                  +     7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     573                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     574                     ) *                                                    &
     575                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
     576              -      (  8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     577                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     578                     ) *                                                    &
     579                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
     580              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     581                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
     582                                 )
     583
     584          diss_t(k) = -ABS( w(k,j,i) ) * (                                  &
     585                     ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     586                  +     3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     587                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     588                     ) *                                                    &
     589                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
     590              -      (  5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     591                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     592                     ) *                                                    &
     593                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
     594              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     595                     ) *                                                    &
     596                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
     597                                         )
     598!
     599!--       Calculate the divergence of the velocity field. A respective
     600!--       correction is needed to overcome numerical instabilities introduced
     601!--       by a not sufficient reduction of divergences near topography.
     602          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
     603                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
     604                        + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     605
     606          tend(k,j,i) = tend(k,j,i) - (                                       &
     607                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
     608                          swap_diss_x_local(k,j,tn)            ) * ddx        &
     609                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
     610                          swap_diss_y_local(k,tn)              ) * ddy        &
     611                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
     612                                                               ) * ddzw(k)    &
     613                                      ) + sk(k,j,i) * div
     614
     615          swap_flux_y_local(k,tn)   = flux_n(k)
     616          swap_diss_y_local(k,tn)   = diss_n(k)
     617          swap_flux_x_local(k,j,tn) = flux_r(k)
     618          swap_diss_x_local(k,j,tn) = diss_r(k)
     619          flux_d                    = flux_t(k)
     620          diss_d                    = diss_t(k)
    619621
    620622       ENDDO
    621623
    622624!
    623 !--    WS3 as an intermediate step (top)
    624        k         = nzt - 1
    625        flux_d    = flux_t(k-1)
    626        diss_d    = diss_t(k-1)
    627        flux_t(k) = w(k,j,i) * (                                                &
    628                                 7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )            &
    629                                     - ( sk(k+2,j,i) + sk(k-1,j,i) )            &
    630                               ) * adv_sca_3
    631        diss_t(k) = -ABS( w(k,j,i) ) * (                                        &
    632                                         3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )    &
    633                                             - ( sk(k+2,j,i) - sk(k-1,j,i) )    &
    634                                       ) * adv_sca_3
    635 
    636        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
    637                                    * ddzw(k)
    638 !                                 
    639 !--    2nd-order scheme (top)
    640        k         = nzt
    641        flux_d    = flux_t(k-1)
    642        diss_d    = diss_t(k-1)
    643        flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
    644 
    645 !       
    646 !--    sk(k+1) is referenced two times to avoid a segmentation fault at top
    647        diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i), &
    648                              sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
    649 
    650        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
    651                                    * ddzw(k)
    652 !       
    653625!--    Evaluation of statistics
    654626       SELECT CASE ( sk_char )
    … …  
    656628          CASE ( 'pt' )
    657629
    658              DO  k = nzb_s_inner(j,i), nzt
    659                 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                      &
     630             DO  k = nzb, nzt
     631                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                &
    660632                                       ( flux_t(k) + diss_t(k) )               &
    661633                                 * weight_substep(intermediate_timestep_count)
    … …  
    664636          CASE ( 'sa' )
    665637
    666              DO  k = nzb_s_inner(j,i), nzt
    667                 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                      &
     638             DO  k = nzb, nzt
     639                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                &
    668640                                       ( flux_t(k) + diss_t(k) )               &
    669641                                 * weight_substep(intermediate_timestep_count)
    … …  
    672644          CASE ( 'q' )
    673645
    674              DO  k = nzb_s_inner(j,i), nzt
    675                 sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                       &
     646             DO  k = nzb, nzt
     647                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                 &
    676648                                      ( flux_t(k) + diss_t(k) )                &
    677649                                 * weight_substep(intermediate_timestep_count)
    … …  
    699671       IMPLICIT NONE
    700672
    701        INTEGER ::  i, i_omp, j, k, tn
    702        LOGICAL ::  degraded_l, degraded_s
    703        REAL    ::  gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp
     673       INTEGER ::  i, i_omp, j, k, tn, k_ppp, k_pp, k_mm
     674       REAL    ::  gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp, div
    704675       REAL, DIMENSION(nzb:nzt+1) ::  flux_t, diss_t, flux_r, diss_r,        &
    705676                                      flux_n, diss_n, u_comp
    706                                        
    707        degraded_l = .FALSE.
    708        degraded_s = .FALSE.
    709677
    710678       gu = 2.0 * u_gtrans
    711679       gv = 2.0 * v_gtrans
    712          
    713        IF ( boundary_flags(j,i) /= 0 )  THEN
    714 !       
    715 !--      Degrade the order for Dirichlet bc. at the outflow boundary
    716          SELECT CASE ( boundary_flags(j,i) )
    717          
    718             CASE ( 1 )
    719                DO  k = nzb_u_inner(j,i)+1, nzt
    720                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
    721                   flux_r(k) = ( u_comp(k) - gu ) * (                          &
    722                               7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
    723                             -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    724                   diss_r(k) = - ABS( u_comp(k) - gu ) * (                     &
    725                               3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
    726                             -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
    727               ENDDO
    728              
    729             CASE ( 2 )
    730                DO  k = nzb_u_inner(j,i)+1, nzt
    731                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
    732                   flux_r(k) = ( u_comp(k) - gu ) * (                          &
    733                                 u(k,j,i+1) + u(k,j,i) ) * 0.25
    734                   diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i),     &
    735                                         u(k,j,i-1), u(k,j,i-2), u_comp(k),    &
    736                                         0.25, ddx )
    737                ENDDO
    738                
    739             CASE ( 3 )
    740                DO  k = nzb_u_inner(j,i)+1, nzt
    741                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    742                   flux_n(k) = v_comp * (                                      &
    743                               7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
    744                             -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    745                   diss_n(k) = - ABS( v_comp ) * (                             &
    746                               3.0 * ( u(k,j+1,i) - u(k,j,i)     )             &
    747                             -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    748                ENDDO
    749                
    750             CASE ( 4 )
    751                DO  k = nzb_u_inner(j,i)+1, nzt
    752                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    753                   flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
    754                   diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), u(k,j,i),     &
    755                                         u(k,j-1,i), u(k,j-2,i), v_comp,       &
    756                                         0.25, ddy )
    757                ENDDO
    758                
    759             CASE ( 5 )
    760                DO  k = nzb_u_inner(j,i)+1, nzt
    761 !               
    762 !--               Compute leftside fluxes for the left boundary of PE domain
    763                   u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
    764                   flux_l_u(k,j,tn) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
    765                   diss_l_u(k,j,tn) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
    766                                             u(k,j,i-1), u(k,j,i-1), u_comp(k),&
    767                                             0.25, ddx )
    768                  
    769                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
    770                   flux_r(k) = ( u_comp(k) - gu ) * (                          &
    771                               7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
    772                             -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    773                   diss_r(k) = - ABS( u_comp(k) -gu ) * (                      &
    774                               3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
    775                             -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
    776                ENDDO
    777                degraded_l = .TRUE.
    778                
    779             CASE ( 7 )
    780                DO  k = nzb_u_inner(j,i)+1, nzt                           
    781                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    782                   flux_n(k) = v_comp * (                                      &
    783                               7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
    784                             -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    785                   diss_n(k) = - ABS( v_comp ) * (                             &
    786                               3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
    787                             -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    788                ENDDO
    789                
    790             CASE ( 8 )
    791                DO  k = nzb_u_inner(j,i)+1, nzt
    792 !               
    793 !--              Compute southside fluxes for the south boundary of PE domain
    794                   v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    795                   flux_s_u(k,tn) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25
    796                   diss_s_u(k,tn) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
    797                                           u(k,j-1,i), u(k,j-1,i), v_comp,     &
    798                                           0.25, ddy )
    799                                  
    800                   v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    801                   flux_n(k) = v_comp * (                                      &
    802                               7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
    803                             -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    804                   diss_n(k) = - ABS( v_comp ) * (                             &
    805                               3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
    806                             -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    807                ENDDO 
    808                degraded_s = .TRUE.
    809                
    810             CASE DEFAULT
    811            
    812          END SELECT
    813 !         
    814 !--      Compute the crosswise 5th order fluxes at the outflow
    815          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
    816               boundary_flags(j,i) == 5 )  THEN
    817          
    818              DO  k = nzb_u_inner(j,i)+1, nzt
    819                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    820                flux_n(k) = v_comp * (                                         &
    821                            37.0 * ( u(k,j+1,i) + u(k,j,i)   )                 &
    822                          -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                 &
    823                          +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    824                diss_n(k) = - ABS ( v_comp ) * (                               &
    825                            10.0 * ( u(k,j+1,i) - u(k,j,i)   )                 &
    826                          -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                 &
    827                          +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    828              ENDDO
    829              
    830          ELSE
    831          
    832             DO  k = nzb_u_inner(j,i)+1, nzt
    833                u_comp(k) = u(k,j,i+1) + u(k,j,i)
    834                flux_r(k) = ( u_comp(k) - gu ) * (                             &
    835                            37.0 * ( u(k,j,i+1) + u(k,j,i)   )                 &
    836                          -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                 &
    837                          +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    838                diss_r(k) = - ABS( u_comp(k) - gu ) * (                        &
    839                            10.0 * ( u(k,j,i+1) - u(k,j,i) )                   &
    840                          -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                 &
    841                          +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    842             ENDDO
    843            
    844          ENDIF
    845                  
    846        ELSE
    847 !       
    848 !--       Compute the fifth order fluxes for the interior of PE domain.
    849           DO  k = nzb_u_inner(j,i)+1, nzt
    850              u_comp(k) = u(k,j,i+1) + u(k,j,i)
    851              flux_r(k) = ( u_comp(k) - gu ) * (                               &
     680!
     681!--    Compute southside fluxes for the respective boundary of PE
     682       IF ( j == nys  )  THEN
     683       
     684          DO  k = nzb+1, nzb_max
     685
     686             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
     687             flux_s_u(k,tn) = v_comp * (                                      &
     688                         ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     689                      +     7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     690                      +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
     691                         ) *                                                  &
     692                                     ( u(k,j,i)   + u(k,j-1,i) )              &
     693                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     694                      +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     695                         ) *                                                  &
     696                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
     697                  +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     698                         ) *                                                  &
     699                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
     700                                       )
     701
     702             diss_s_u(k,tn) = - ABS ( v_comp ) * (                            &
     703                         ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     704                      +     3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     705                      +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
     706                         ) *                                                  &
     707                                     ( u(k,j,i)   - u(k,j-1,i) )              &
     708                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     709                      +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     710                         ) *                                                  &
     711                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
     712                  +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     713                         ) *                                                  &
     714                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
     715                                                 )
     716
     717          ENDDO
     718
     719          DO  k = nzb_max+1, nzt
     720
     721             v_comp         = v(k,j,i) + v(k,j,i-1) - gv
     722             flux_s_u(k,tn) = v_comp * (                                      &
     723                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
     724                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
     725                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
     726             diss_s_u(k,tn) = - ABS(v_comp) * (                               &
     727                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
     728                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
     729                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
     730
     731          ENDDO
     732         
     733       ENDIF
     734!
     735!--    Compute leftside fluxes for the respective boundary of PE
     736       IF ( i == i_omp )  THEN
     737       
     738          DO  k = nzb+1, nzb_max
     739
     740             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     741             flux_l_u(k,j,tn) = u_comp_l * (                                   &
     742                          ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     743                       +     7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     744                       +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1 &
     745                          ) *                                                  &
     746                                     ( u(k,j,i)   + u(k,j,i-1) )               &
     747                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     748                       +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     749                          ) *                                                  &
     750                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
     751                   +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     752                          ) *                                                  &
     753                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
     754                                           )
     755
     756              diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                         &
     757                          ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     758                       +     3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     759                       +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1 &
     760                          ) *                                                  &
     761                                     ( u(k,j,i)   - u(k,j,i-1) )               &
     762                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     763                       +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     764                          ) *                                                  &
     765                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
     766                   +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     767                          ) *                                                  &
     768                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
     769                                                     )
     770
     771          ENDDO
     772
     773          DO  k = nzb_max+1, nzt
     774
     775             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     776             flux_l_u(k,j,tn) = u_comp_l * (                                   &
     777                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )                &
     778                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )                &
     779                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
     780             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                            &
     781                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )                &
     782                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )                &
     783                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
     784
     785          ENDDO
     786         
     787       ENDIF
     788
     789       flux_t(0) = 0.0
     790       diss_t(0) = 0.0
     791       flux_d    = 0.0
     792       diss_d    = 0.0
     793!
     794!--    Now compute the fluxes tendency terms for the horizontal and
     795!--    vertical parts.
     796       DO  k = nzb+1, nzb_max
     797
     798          u_comp(k) = u(k,j,i+1) + u(k,j,i)
     799          flux_r(k) = ( u_comp(k) - gu ) * (                                 &
     800                     ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     801                  +     7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     802                  +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1    &
     803                     ) *                                                     &
     804                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
     805              -      (  8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     806                  +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     807                     ) *                                                     &
     808                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
     809              +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     810                     ) *                                                     &
     811                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
     812                                           )
     813
     814          diss_r(k) = - ABS( u_comp(k) - gu ) * (                            &
     815                     ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     816                  +     3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     817                  +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1    &
     818                     ) *                                                     &
     819                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
     820              -      (  5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     821                  +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     822                     ) *                                                     &
     823                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
     824              +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     825                     ) *                                                     &
     826                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
     827                                                )
     828
     829          v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
     830          flux_n(k) = v_comp * (                                             &
     831                     ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     832                  +     7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     833                  +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1    &
     834                     ) *                                                     &
     835                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
     836              -      (  8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     837                  +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     838                     ) *                                                     &
     839                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
     840              +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     841                     ) *                                                     &
     842                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
     843                                           )
     844
     845          diss_n(k) = - ABS ( v_comp ) * (                                   &
     846                     ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     847                  +     3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     848                  +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1    &
     849                     ) *                                                     &
     850                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
     851              -      (  5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     852                  +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     853                     ) *                                                     &
     854                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
     855              +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     856                     ) *                                                     &
     857                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
     858                                                )
     859!
     860!--       k index has to be modified near bottom and top, else array
     861!--       subscripts will be exceeded.
     862          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
     863          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1)  )
     864          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
     865
     866          w_comp    = w(k,j,i) + w(k,j,i-1)
     867          flux_t(k) = w_comp  * (                                            &
     868                     ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     869                  +     7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     870                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     871                     ) *                                                     &
     872                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
     873              -      (  8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     874                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     875                     ) *                                                     &
     876                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
     877              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     878                     ) *                                                     &
     879                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     880                                 )
     881
     882          diss_t(k) = - ABS( w_comp ) * (                                    &
     883                     ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     884                  +     3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     885                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     886                     ) *                                                     &
     887                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
     888              -      (  5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     889                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     890                     ) *                                                     &
     891                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
     892              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     893                     ) *                                                     &
     894                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     895                                         )
     896!
     897!--       Calculate the divergence of the velocity field. A respective
     898!--       correction is needed to overcome numerical instabilities introduced
     899!--       by a not sufficient reduction of divergences near topography.
     900          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx      &
     901               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy      &
     902               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)  &
     903                ) * 0.5
     904
     905          tend(k,j,i) = tend(k,j,i) - (                                       &
     906                            ( flux_r(k) + diss_r(k)                           &
     907                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
     908                          + ( flux_n(k) + diss_n(k)                           &
     909                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
     910                          + ( flux_t(k) + diss_t(k)                           &
     911                          -   flux_d    - diss_d                  ) * ddzw(k) &
     912                                       ) + div * u(k,j,i)
     913
     914           flux_l_u(k,j,tn) = flux_r(k)
     915           diss_l_u(k,j,tn) = diss_r(k)
     916           flux_s_u(k,tn)   = flux_n(k)
     917           diss_s_u(k,tn)   = diss_n(k)
     918           flux_d           = flux_t(k)
     919           diss_d           = diss_t(k)
     920!
     921!--        Statistical Evaluation of u'u'. The factor has to be applied for
     922!--        right evaluation when gallilei_trans = .T. .
     923           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
     924                              + ( flux_r(k) *                                 &
     925                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
     926                              / ( u_comp(k) - gu + 1.0E-20      )             &
     927                              +   diss_r(k) *                                 &
     928                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
     929                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
     930                              *   weight_substep(intermediate_timestep_count)
     931!
     932!--        Statistical Evaluation of w'u'.
     933           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                        &
     934                              + ( flux_t(k) + diss_t(k) )                     &
     935                              *   weight_substep(intermediate_timestep_count)
     936       ENDDO
     937
     938       DO  k = nzb_max+1, nzt
     939
     940          u_comp(k) = u(k,j,i+1) + u(k,j,i)
     941          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
    852942                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
    853943                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
    … …  
    867957                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
    868958                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    869           ENDDO
    870          
    871        ENDIF
    872 !       
    873 !--    Compute left- and southside fluxes for the respective boundary of PE
    874        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
    875        
    876           DO  k = nzb_u_inner(j,i)+1, nzt
    877              v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    878              flux_s_u(k,tn) = v_comp * (                                         &
    879                            37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
    880                          -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
    881                          +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    882              diss_s_u(k,tn) = - ABS(v_comp) * (                                  &
    883                            10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
    884                          -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
    885                          +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
    886           ENDDO
    887          
    888        ENDIF
    889        
    890        IF ( i == i_omp .AND. ( .NOT. degraded_l ) )  THEN
    891        
    892           DO  k = nzb_u_inner(j,i)+1, nzt
    893              u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
    894              flux_l_u(k,j,tn) = u_comp_l * (                                     &
    895                              37.0 * ( u(k,j,i) + u(k,j,i-1)   )               &
    896                            -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )               &
    897                            +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
    898              diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                              &
    899                              10.0 * ( u(k,j,i) - u(k,j,i-1)   )               &
    900                            -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )               &
    901                            +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
    902           ENDDO
    903          
    904        ENDIF
    905 !       
    906 !--    Now compute the tendency terms for the horizontal parts.
    907        DO  k = nzb_u_inner(j,i)+1, nzt                   
     959!
     960!--       k index has to be modified near bottom and top, else array
     961!--       subscripts will be exceeded.
     962          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
     963          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1)  )
     964          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
     965
     966          w_comp    = w(k,j,i) + w(k,j,i-1)
     967          flux_t(k) = w_comp  * (                                            &
     968                     ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     969                  +     7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     970                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     971                     ) *                                                     &
     972                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
     973              -      (  8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     974                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     975                     ) *                                                     &
     976                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
     977              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     978                     ) *                                                     &
     979                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     980                                 )
     981
     982          diss_t(k) = - ABS( w_comp ) * (                                    &
     983                     ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     984                  +     3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     985                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     986                     ) *                                                     &
     987                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
     988              -      (  5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     989                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     990                     ) *                                                     &
     991                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
     992              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     993                     ) *                                                     &
     994                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     995                                         )
     996!
     997!--       Calculate the divergence of the velocity field. A respective
     998!--       correction is needed to overcome numerical instabilities introduced
     999!--       by a not sufficient reduction of divergences near topography.
     1000          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx      &
     1001               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy      &
     1002               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)  &
     1003                ) * 0.5
     1004
    9081005          tend(k,j,i) = tend(k,j,i) - (                                       &
    9091006                            ( flux_r(k) + diss_r(k)                           &
    910                           -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx           &
     1007                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
    9111008                          + ( flux_n(k) + diss_n(k)                           &
    912                           -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy  )
     1009                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
     1010                          + ( flux_t(k) + diss_t(k)                           &
     1011                          -   flux_d    - diss_d                  ) * ddzw(k) &
     1012                                       ) + div * u(k,j,i)
    9131013
    9141014           flux_l_u(k,j,tn) = flux_r(k)
    … …  
    9161016           flux_s_u(k,tn)   = flux_n(k)
    9171017           diss_s_u(k,tn)   = diss_n(k)
    918 !
    919 !--        Statistical Evaluation of u'u'. The factor has to be applied for
     1018           flux_d           = flux_t(k)
     1019           diss_d           = diss_t(k)
     1020!
     1021!--        Statistical Evaluation of u'u'. The factor has to be applied for
    9201022!--        right evaluation when gallilei_trans = .T. .
    921            sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                                &
    922              + ( flux_r(k) *                                                  &
    923                ( u_comp(k) - 2.0 * hom(k,1,1,0) )                             &
    924              / ( u_comp(k) - gu + 1.0E-20      )                              &
    925              +   diss_r(k) *                                                  &
    926                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )                        &
    927              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )                          &
    928              *   weight_substep(intermediate_timestep_count)
    929        ENDDO
    930        sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
    931                                            
    932 
    933 !
    934 !--    Vertical advection, degradation of order near surface and top.
    935 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    936 !--    statistical evaluation the top flux at the surface should be 0
    937        flux_t(nzb_u_inner(j,i)) = 0. !statistical reasons
    938        diss_t(nzb_u_inner(j,i)) = 0.
    939 !
    940 !--    2nd order scheme (bottom)
    941        k         = nzb_u_inner(j,i)+1
    942        flux_d    = flux_t(k-1)
    943        diss_d    = diss_t(k-1)
    944        w_comp    = w(k,j,i) + w(k,j,i-1)
    945        flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) *0.25
    946        diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), 0., 0.,        &
    947                              w_comp, 0.25, ddzw(k) )
    948              
    949        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    950                                  -   flux_d    - diss_d      ) * ddzw(k)
    951 !
    952 !--    WS3 as an intermediate step (bottom)
    953        k         = nzb_u_inner(j,i)+2
    954        flux_d    = flux_t(k-1)
    955        diss_d    = diss_t(k-1)
    956        w_comp    = w(k,j,i) + w(k,j,i-1)
    957        flux_t(k) = w_comp * (                                                 &
    958                    7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
    959                  -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    960        diss_t(k) = -ABS( w_comp) * (                                          &
    961                    3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          &
    962                  -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    963 
    964        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    965                                  -   flux_d    - diss_d      ) * ddzw(k)
    966 !
    967 !--    WS5
    968        DO  k = nzb_u_inner(j,i)+3, nzt-2
    969           flux_d    = flux_t(k-1)
    970           diss_d    = diss_t(k-1)
    971           w_comp    = w(k,j,i) + w(k,j,i-1)
    972           flux_t(k) = w_comp * (                                              &
    973                       37.0 * ( u(k+1,j,i) + u(k,j,i)   )                      &
    974                     -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                      &
    975                     +        ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
    976           diss_t(k) = - ABS( w_comp ) * (                                     &
    977                       10.0 *  ( u(k+1,j,i) - u(k,j,i)   )                     &
    978                     -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                      &
    979                     +        ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
    980 
    981           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    982                                     -   flux_d    - diss_d      ) * ddzw(k)
     1023           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
     1024                              + ( flux_r(k) *                                 &
     1025                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
     1026                              / ( u_comp(k) - gu + 1.0E-20      )             &
     1027                              +   diss_r(k) *                                 &
     1028                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
     1029                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
     1030                              *   weight_substep(intermediate_timestep_count)
     1031!
     1032!--        Statistical Evaluation of w'u'.
     1033           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                        &
     1034                              + ( flux_t(k) + diss_t(k) )                     &
     1035                              *   weight_substep(intermediate_timestep_count)
    9831036       ENDDO
    984 !
    985 !--    WS3 as an intermediate step (top)
    986        k         = nzt - 1
    987        flux_d    = flux_t(k-1)
    988        diss_d    = diss_t(k-1)
    989        w_comp    = w(k,j,i) + w(k,j,i-1)
    990        flux_t(k) = w_comp * (                                                 &
    991                    7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
    992                  -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    993        diss_t(k) = - ABS( w_comp ) * (                                        &
    994                    3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          &
    995                  -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    996                    
    997        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    998                                  -   flux_d    - diss_d      ) * ddzw(k)
    999        
    1000 !
    1001 !--    2nd order scheme (top)
    1002        k         = nzt
    1003        flux_d    = flux_t(k-1)
    1004        diss_d    = diss_t(k-1)
    1005        w_comp    = w(k,j,i)+w(k,j,i-1)
    1006        flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
    1007        diss_t(k) = diss_2nd( u(k+1,j,i), u(k+1,j,i), u(k,j,i), u(k-1,j,i),    &
    1008                              u(k-2,j,i), w_comp, 0.25, ddzw(k) )
    1009 
    1010        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1011                                  -   flux_d    - diss_d      ) * ddzw(k)
    1012 
    1013 !
    1014 !--    sum up the vertical momentum fluxes
    1015        DO  k = nzb_u_inner(j,i), nzt
    1016           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                               &
    1017               + ( flux_t(k) + diss_t(k) )                                     &
    1018               *   weight_substep(intermediate_timestep_count)
    1019        ENDDO
     1037
     1038       sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
     1039
    10201040
    10211041
    10221042    END SUBROUTINE advec_u_ws_ij
    1023 
    10241043
    10251044
    … …  
    10391058       IMPLICIT NONE
    10401059
    1041        INTEGER  ::  i, i_omp, j, k, tn
    1042        LOGICAL  ::  degraded_l, degraded_s
    1043        REAL     ::  gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp
     1060       INTEGER  ::  i, i_omp, j, k, tn, k_ppp, k_pp, k_mm
     1061       REAL     ::  gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp, div
    10441062       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_n,                &
    10451063                                       diss_n, flux_r, diss_r, v_comp
    1046                                        
    1047        degraded_l = .FALSE.
    1048        degraded_s = .FALSE.
    1049      
     1064
    10501065       gu = 2.0 * u_gtrans
    10511066       gv = 2.0 * v_gtrans
    1052        
    1053        IF ( boundary_flags(j,i) /= 0 )  THEN
     1067
    10541068!       
    1055 !--       Degrade the order for Dirichlet bc. at the outflow boundary
    1056           SELECT CASE ( boundary_flags(j,i) )
    1057          
    1058              CASE ( 1 )
    1059                 DO  k = nzb_v_inner(j,i)+1, nzt
    1060                    u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1061                    flux_r(k) = u_comp * (                                     &
    1062                                7.0 * (v(k,j,i+1) + v(k,j,i)    )              &
    1063                              -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    1064                    diss_r(k) = - ABS( u_comp ) * (                            &
    1065                                3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
    1066                              -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    1067                 ENDDO
    1068                
    1069              CASE ( 2 )
    1070                 DO  k = nzb_v_inner(j,i)+1, nzt
    1071                    u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1072                    flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25
    1073                    diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), v(k,j,i),    &
    1074                                          v(k,j,i-1), v(k,j,i-2), u_comp,      &
    1075                                          0.25, ddx )
    1076                 ENDDO
    1077                
    1078              CASE ( 3 )
    1079                 DO  k = nzb_v_inner(j,i)+1, nzt
    1080                    v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1081                    flux_n(k) = ( v_comp(k)- gv ) * (                          &
    1082                                7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
    1083                              -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    1084                    diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
    1085                                3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
    1086                              -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
    1087                 ENDDO
    1088                
    1089              CASE ( 4 )
    1090                 DO  k = nzb_v_inner(j,i)+1, nzt
    1091                    v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1092                    flux_n(k) = ( v_comp(k) - gv ) *                           &
    1093                                ( v(k,j+1,i) + v(k,j,i) ) * 0.25
    1094                    diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), v(k,j,i),    &
    1095                                          v(k,j-1,i), v(k,j-2,i), v_comp(k),   &
    1096                                          0.25, ddy )
    1097                 ENDDO
    1098                
    1099              CASE ( 5 )
    1100                 DO  k = nzb_w_inner(j,i)+1, nzt
    1101                    u_comp    = u(k,j-1,i) + u(k,j,i) - gu
    1102                    flux_r(k) = u_comp * (                                     &
    1103                                7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
    1104                              -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    1105                    diss_r(k) = - ABS( u_comp ) * (                            &
    1106                                3.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
    1107                              -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    1108                 ENDDO
    1109                
    1110              CASE ( 6 )
    1111                 DO  k = nzb_v_inner(j,i)+1, nzt
    1112                    u_comp        = u(k,j-1,i) + u(k,j,i) - gu
    1113                    flux_l_v(k,j,tn) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
    1114                    diss_l_v(k,j,tn) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
    1115                                              v(k,j,i-1), v(k,j,i-1), u_comp,  &
    1116                                              0.25, ddx )
    1117                                  
    1118                    u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1119                    flux_r(k) = u_comp * (                                     &
    1120                                7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
    1121                              -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    1122                    diss_r(k) = - ABS( u_comp ) * (                            &
    1123                                3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
    1124                              -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    1125                 ENDDO
    1126                 degraded_l = .TRUE.
    1127                
    1128              CASE ( 7 )
    1129                 DO  k = nzb_v_inner(j,i)+1, nzt
    1130                    v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
    1131                    flux_s_v(k,tn) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
    1132                    diss_s_v(k,tn) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
    1133                                            v(k,j-1,i), v(k,j-1,i), v_comp(k), &
    1134                                            0.25, ddy )
    1135                                  
    1136                    v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1137                    flux_n(k) = ( v_comp(k) - gv ) * (                         &
    1138                                7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
    1139                              -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    1140                    diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
    1141                                3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
    1142                              -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
    1143                 ENDDO
    1144                 degraded_s = .TRUE.
    1145                
    1146               CASE DEFAULT
    1147              
    1148           END SELECT
    1149 !         
    1150 !--       Compute the crosswise 5th order fluxes at the outflow
    1151           IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.   &
    1152                boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
    1153          
    1154              DO  k = nzb_v_inner(j,i)+1, nzt
    1155                 v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1156                 flux_n(k) = ( v_comp(k) - gv ) * (                            &
    1157                             37.0 * ( v(k,j+1,i) + v(k,j,i)   )                &
    1158                           -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                &
    1159                           +        ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    1160                 diss_n(k) = - ABS( v_comp(k) - gv ) * (                       &
    1161                              10.0 * ( v(k,j+1,i) - v(k,j,i)   )               &
    1162                            -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )               &
    1163                            +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    1164               ENDDO
    1165              
    1166           ELSE
    1167          
    1168              DO  k = nzb_v_inner(j,i)+1, nzt
    1169                 u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1170                 flux_r(k) = u_comp * (                                        &
    1171                             37.0 * ( v(k,j,i+1) + v(k,j,i)   )                &
    1172                           -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                &
    1173                           +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    1174                 diss_r(k) = - ABS( u_comp ) * (                               &
    1175                             10.0 * ( v(k,j,i+1) - v(k,j,i)   )                &
    1176                           -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                &
    1177                           +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    1178              ENDDO
    1179              
    1180           ENDIF
    1181          
    1182        ELSE
    1183 !       
    1184 !--       Compute the fifth order fluxes for the interior of PE domain.                 
    1185           DO  k = nzb_v_inner(j,i)+1, nzt
    1186              u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1187              flux_r(k) = u_comp * (                                           &
    1188                          37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
    1189                        -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                   &
    1190                        +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    1191              diss_r(k) = - ABS( u_comp ) * (                                  &
    1192                          10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
    1193                        -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
    1194                        +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    1195 
    1196              v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1197              flux_n(k) = ( v_comp(k) - gv ) * (                               &
    1198                          37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
    1199                        -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                   &
    1200                          +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    1201              diss_n(k) = - ABS( v_comp(k) - gv ) * (                          &
    1202                          10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
    1203                        -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
    1204                        +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    1205           ENDDO
    1206          
    1207        ENDIF
    1208 !       
    1209 !--    Compute left- and southside fluxes for the respective boundary       
    1210        IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
    1211 
    1212           DO  k = nzb_v_inner(j,i)+1, nzt
    1213              u_comp        = u(k,j-1,i) + u(k,j,i) - gu
    1214              flux_l_v(k,j,tn) = u_comp * (                                       &
     1069!--    Compute leftside fluxes for the respective boundary.
     1070       IF ( i == i_omp )  THEN
     1071
     1072          DO  k = nzb+1, nzb_max
     1073
     1074             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
     1075             flux_l_v(k,j,tn) = u_comp * (                                     &
     1076                          ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
     1077                       +     7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
     1078                       +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 &
     1079                          ) *                                                  &
     1080                                     ( v(k,j,i)   + v(k,j,i-1) )               &
     1081                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
     1082                       +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
     1083                          ) *                                                  &
     1084                                     ( v(k,j,i+1) + v(k,j,i-2) )               &
     1085                   +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
     1086                          ) *                                                  &
     1087                                     ( v(k,j,i+2) + v(k,j,i-3) )               &
     1088                                         )
     1089
     1090              diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
     1091                          ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
     1092                       +     3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
     1093                       +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 &
     1094                          ) *                                                  &
     1095                                     ( v(k,j,i)   - v(k,j,i-1) )               &
     1096                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
     1097                       +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
     1098                          ) *                                                  &
     1099                                     ( v(k,j,i+1) - v(k,j,i-2) )               &
     1100                   +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
     1101                          ) *                                                  &
     1102                                     ( v(k,j,i+2) - v(k,j,i-3) )               &
     1103                                                   )
     1104
     1105          ENDDO
     1106
     1107          DO  k = nzb_max+1, nzt
     1108
     1109             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
     1110             flux_l_v(k,j,tn) = u_comp * (                                    &
    12151111                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
    12161112                           -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
    12171113                           +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    1218              diss_l_v(k,j,tn) = - ABS( u_comp ) * (                              &
     1114             diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
    12191115                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
    12201116                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
    12211117                           +        ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
     1118
    12221119          ENDDO
    12231120         
    12241121       ENDIF
     1122!
     1123!--    Compute southside fluxes for the respective boundary.
     1124       IF ( j == nysv )  THEN
    12251125       
    1226        IF ( j == nysv  .AND.  .NOT. degraded_s )  THEN
    1227        
    1228           DO  k = nzb_v_inner(j,i)+1, nzt
    1229              v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
    1230              flux_s_v(k,tn) = v_comp_l * (                                       &
     1126          DO  k = nzb+1, nzb_max
     1127
     1128             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     1129             flux_s_v(k,tn) = v_comp_l * (                                    &
     1130                         ( 37.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
     1131                      +     7.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
     1132                      +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 &
     1133                         ) *                                                  &
     1134                                     ( v(k,j,i)   + v(k,j-1,i) )              &
     1135                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
     1136                      +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
     1137                         ) *                                                  &
     1138                                     ( v(k,j+1,i) + v(k,j-2,i) )              &
     1139                  +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
     1140                         ) *                                                  &
     1141                                     ( v(k,j+2,i) + v(k,j-3,i) )              &
     1142                                         )
     1143
     1144             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
     1145                         ( 10.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
     1146                      +     3.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
     1147                      +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 &
     1148                         ) *                                                  &
     1149                                     ( v(k,j,i)   - v(k,j-1,i) )              &
     1150                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
     1151                      +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
     1152                         ) *                                                  &
     1153                                     ( v(k,j+1,i) - v(k,j-2,i) )              &
     1154                  +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
     1155                         ) *                                                  &
     1156                                     ( v(k,j+2,i) - v(k,j-3,i) )              &
     1157                                                 )
     1158
     1159          ENDDO
     1160
     1161          DO  k = nzb_max+1, nzt
     1162
     1163             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     1164             flux_s_v(k,tn) = v_comp_l * (                                    &
    12311165                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
    12321166                         -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
    12331167                         +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
    1234              diss_s_v(k,tn) = - ABS( v_comp_l ) * (                              &
     1168             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    12351169                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
    12361170                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
    12371171                         +        ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
     1172
    12381173          ENDDO
    12391174         
    12401175       ENDIF
    1241 !       
    1242 !--    Now compute the tendency terms for the horizontal parts.         
    1243        DO  k = nzb_v_inner(j,i)+1, nzt                 
     1176
     1177       flux_t(0) = 0.0
     1178       diss_t(0) = 0.0
     1179       flux_d    = 0.0
     1180       diss_d    = 0.0
     1181!
     1182!--    Now compute the fluxes and tendency terms for the horizontal and
     1183!--    verical parts.
     1184       DO  k = nzb+1, nzb_max
     1185
     1186          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
     1187          flux_r(k) = u_comp * (                                             &
     1188                     ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
     1189                  +     7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
     1190                  +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1    &
     1191                     ) *                                                     &
     1192                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
     1193              -      (  8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
     1194                  +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
     1195                     ) *                                                     &
     1196                                 ( v(k,j,i+2) + v(k,j,i-1) )                 &
     1197              +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
     1198                     ) *                                                     &
     1199                                 ( v(k,j,i+3) + v(k,j,i-2) )                 &
     1200                               )
     1201
     1202          diss_r(k) = - ABS( u_comp ) * (                                    &
     1203                     ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
     1204                  +     3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
     1205                  +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1    &
     1206                     ) *                                                     &
     1207                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
     1208              -      (  5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
     1209                  +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
     1210                     ) *                                                     &
     1211                                 ( v(k,j,i+2) - v(k,j,i-1) )                 &
     1212              +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
     1213                     ) *                                                     &
     1214                                 ( v(k,j,i+3) - v(k,j,i-2) )                 &
     1215                                        )
     1216
     1217          v_comp(k) = v(k,j+1,i) + v(k,j,i)
     1218          flux_n(k) = ( v_comp(k) - gv ) * (                                 &
     1219                     ( 37.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
     1220                  +     7.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
     1221                  +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1    &
     1222                     ) *                                                     &
     1223                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
     1224              -      (  8.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
     1225                  +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
     1226                     ) *                                                     &
     1227                                 ( v(k,j+2,i) + v(k,j-1,i) )                 &
     1228              +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
     1229                     ) *                                                     &
     1230                                 ( v(k,j+3,i) + v(k,j-2,i) )                 &
     1231                               )
     1232
     1233          diss_n(k) = - ABS( v_comp(k) - gv ) * (                            &
     1234                     ( 10.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
     1235                  +     3.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
     1236                  +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1    &
     1237                     ) *                                                     &
     1238                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
     1239              -      (  5.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
     1240                  +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
     1241                     ) *                                                     &
     1242                                 ( v(k,j+2,i) - v(k,j-1,i) )                 &
     1243              +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
     1244                     ) *                                                     &
     1245                                 ( v(k,j+3,i) - v(k,j-2,i) )                 &
     1246                                        )
     1247!
     1248!--       k index has to be modified near bottom and top, else array
     1249!--       subscripts will be exceeded.
     1250          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),26,1)
     1251          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),24,1)  )
     1252          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),26,1)
     1253
     1254          w_comp    = w(k,j-1,i) + w(k,j,i)
     1255          flux_t(k) = w_comp  * (                                            &
     1256                     ( 37.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1257                  +     7.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1258                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
     1259                     ) *                                                     &
     1260                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
     1261              -      (  8.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1262                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1263                     ) *                                                     &
     1264                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
     1265              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1266                     ) *                                                     &
     1267                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
     1268                                 )
     1269
     1270          diss_t(k) = - ABS( w_comp ) * (                                    &
     1271                     ( 10.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1272                  +     3.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1273                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
     1274                     ) *                                                     &
     1275                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
     1276              -      (  5.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1277                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1278                     ) *                                                     &
     1279                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
     1280              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1281                     ) *                                                     &
     1282                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
     1283                                         )
     1284!
     1285!--       Calculate the divergence of the velocity field. A respective
     1286!--       correction is needed to overcome numerical instabilities introduced
     1287!--       by a not sufficient reduction of divergences near topography.
     1288          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx      &
     1289               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy      &
     1290               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)  &
     1291                ) * 0.5
     1292
    12441293          tend(k,j,i) = tend(k,j,i) - (                                       &
    12451294                         ( flux_r(k) + diss_r(k)                              &
    1246                        -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx            &
     1295                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
    12471296                       + ( flux_n(k) + diss_n(k)                              &
    1248                        -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy     )
    1249        
     1297                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
     1298                       + ( flux_t(k) + diss_t(k)                              &
     1299                       -   flux_d    - diss_d                    ) * ddzw(k)  &
     1300                                      ) + v(k,j,i) * div
     1301
    12501302           flux_l_v(k,j,tn) = flux_r(k)
    12511303           diss_l_v(k,j,tn) = diss_r(k)
    12521304           flux_s_v(k,tn)   = flux_n(k)
    12531305           diss_s_v(k,tn)   = diss_n(k)
    1254 
    1255 !
    1256 !--        Statistical Evaluation of v'v'. The factor has to be applied for
     1306           flux_d           = flux_t(k)
     1307           diss_d           = diss_t(k)
     1308
     1309!
     1310!--        Statistical Evaluation of v'v'. The factor has to be applied for
    12571311!--        right evaluation when gallilei_trans = .T. .
    1258 
    1259            sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                                &
     1312           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
    12601313             + ( flux_n(k)                                                    &
    12611314             * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
    … …  
    12651318             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
    12661319             *   weight_substep(intermediate_timestep_count)
     1320!
     1321!--        Statistical Evaluation of w'v'.
     1322           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                        &
     1323                              + ( flux_t(k) + diss_t(k) )                     &
     1324                              *   weight_substep(intermediate_timestep_count)
    12671325
    12681326       ENDDO
    1269        sums_vs2_ws_l(nzb_v_inner(j,i),tn) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,tn) 
    1270                                            
    1271 !
    1272 !--    Vertical advection, degradation of order near surface and top.
    1273 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    1274 !--    statistical evaluation the top flux at the surface should be 0
    1275        flux_t(nzb_v_inner(j,i)) = 0.0  !statistical reasons
    1276        diss_t(nzb_v_inner(j,i)) = 0.0
    1277 !
    1278 !--    2nd order scheme (bottom)     
    1279        k         = nzb_v_inner(j,i)+1
    1280        flux_d    = flux_t(k-1)
    1281        diss_d    = diss_t(k-1)
    1282        w_comp    = w(k,j-1,i) + w(k,j,i)
    1283        flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
    1284        diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0.0, 0.0,      &
    1285                              w_comp, 0.25, ddzw(k) )
    1286 
    1287        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1288                                  -   flux_d    - diss_d       ) * ddzw(k)
    1289        
    1290 !
    1291 !--    WS3 as an intermediate step (bottom)
    1292        k         = nzb_v_inner(j,i)+2
    1293        flux_d    = flux_t(k-1)
    1294        diss_d    = diss_t(k-1)
    1295        w_comp    = w(k,j-1,i) + w(k,j,i)
    1296        flux_t(k) = w_comp * (                                                 &
    1297                    7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &
    1298                  -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
    1299        diss_t(k) = - ABS( w_comp ) * (                                        &
    1300                    3.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 
    1303        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1304                                  -   flux_d    - diss_d       ) * ddzw(k)
    1305 !                                 
    1306 !--    WS5
    1307        DO  k = nzb_v_inner(j,i)+3, nzt-2
    1308           flux_d    = flux_t(k-1)
    1309           diss_d    = diss_t(k-1)
     1327
     1328       DO  k = nzb_max+1, nzt
     1329
     1330          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
     1331          flux_r(k) = u_comp * (                                           &
     1332                      37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
     1333                    -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                   &
     1334                    +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
     1335
     1336          diss_r(k) = - ABS( u_comp ) * (                                  &
     1337                      10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
     1338                    -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
     1339                    +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     1340
     1341
     1342          v_comp(k) = v(k,j+1,i) + v(k,j,i)
     1343          flux_n(k) = ( v_comp(k) - gv ) * (                               &
     1344                      37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
     1345                    -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                   &
     1346                      +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
     1347
     1348          diss_n(k) = - ABS( v_comp(k) - gv ) * (                          &
     1349                      10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
     1350                    -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
     1351                    +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
     1352!
     1353!--       k index has to be modified near bottom and top, else array
     1354!--       subscripts will be exceeded.
     1355          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),26,1)
     1356          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),24,1)  )
     1357          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),26,1)
     1358
    13101359          w_comp    = w(k,j-1,i) + w(k,j,i)
    1311           flux_t(k) = w_comp * (                                              &
    1312                       37.0 * ( v(k+1,j,i) + v(k,j,i)   )                      &           
    1313                     -  8.0 * ( v(k+2,j,i) + v(k-1,j,i) )                      &
    1314                     +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
    1315           diss_t(k) = - ABS( w_comp ) * (                                     &
    1316                       10.0 * ( v(k+1,j,i) - v(k,j,i)   )                      &
    1317                     -  5.0 * ( v(k+2,j,i) - v(k-1,j,i) )                      &
    1318                     +        ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
    1319 
    1320           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    1321                                  -   flux_d    - diss_d       ) * ddzw(k)
     1360          flux_t(k) = w_comp  * (                                            &
     1361                     ( 37.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1362                  +     7.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1363                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
     1364                     ) *                                                     &
     1365                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
     1366              -      (  8.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1367                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1368                     ) *                                                     &
     1369                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
     1370              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1371                     ) *                                                     &
     1372                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
     1373                                 )
     1374
     1375          diss_t(k) = - ABS( w_comp ) * (                                    &
     1376                     ( 10.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1377                  +     3.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1378                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
     1379                     ) *                                                     &
     1380                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
     1381              -      (  5.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1382                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
     1383                     ) *                                                     &
     1384                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
     1385              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
     1386                     ) *                                                     &
     1387                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
     1388                                         )
     1389!
     1390!--       Calculate the divergence of the velocity field. A respective
     1391!--       correction is needed to overcome numerical instabilities introduced
     1392!--       by a not sufficient reduction of divergences near topography.
     1393          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx      &
     1394               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy      &
     1395               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)  &
     1396                ) * 0.5
     1397
     1398          tend(k,j,i) = tend(k,j,i) - (                                       &
     1399                         ( flux_r(k) + diss_r(k)                              &
     1400                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
     1401                       + ( flux_n(k) + diss_n(k)                              &
     1402                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
     1403                       + ( flux_t(k) + diss_t(k)                              &
     1404                       -   flux_d    - diss_d                    ) * ddzw(k)  &
     1405                                      ) + v(k,j,i) * div
     1406
     1407           flux_l_v(k,j,tn) = flux_r(k)
     1408           diss_l_v(k,j,tn) = diss_r(k)
     1409           flux_s_v(k,tn)   = flux_n(k)
     1410           diss_s_v(k,tn)   = diss_n(k)
     1411           flux_d           = flux_t(k)
     1412           diss_d           = diss_t(k)
     1413
     1414!
     1415!--        Statistical Evaluation of v'v'. The factor has to be applied for
     1416!--        right evaluation when gallilei_trans = .T. .
     1417           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
     1418             + ( flux_n(k)                                                    &
     1419             * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
     1420             / ( v_comp(k) - gv + 1.0E-20 )                                   &
     1421             +   diss_n(k)                                                    &
     1422             *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
     1423             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
     1424             *   weight_substep(intermediate_timestep_count)
     1425!
     1426!--        Statistical Evaluation of w'v'.
     1427           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
     1428                              + ( flux_t(k) + diss_t(k) )                      &
     1429                              *   weight_substep(intermediate_timestep_count)
     1430
    13221431       ENDDO
    1323 !
    1324 !--    WS3 as an intermediate step (top)
    1325        k         = nzt - 1
    1326        flux_d    = flux_t(k-1)
    1327        diss_d    = diss_t(k-1)
    1328        w_comp    = w(k,j-1,i) + w(k,j,i)
    1329        flux_t(k) = w_comp * (                                                 &
    1330                    7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &   
    1331                  -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
    1332        diss_t(k) = - ABS( w_comp ) * (                                        &
    1333                    3.0 * ( v(k+1,j,i) - v(k,j,i) )                            &
    1334                  -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
    1335        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1336                                  -   flux_d    - diss_d       ) * ddzw(k)
    1337 !
    1338 !--    2nd order scheme (top)
    1339        k         = nzt
    1340        flux_d    = flux_t(k-1)
    1341        diss_d    = diss_t(k-1)
    1342        w_comp    = w(k,j-1,i)+w(k,j,i)
    1343        flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
    1344        diss_t(k) = diss_2nd( v(k+1,j,i), v(k+1,j,i), v(k,j,i), v(k-1,j,i),    &
    1345                              v(k-2,j,i), w_comp, 0.25, ddzw(k) )
    1346  
    1347        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1348                                  -   flux_d    - diss_d       ) * ddzw(k)
    1349              
    1350        DO  k = nzb_v_inner(j,i), nzt
    1351           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                               &
    1352                  + ( flux_t(k) + diss_t(k) )                                  &
    1353                  *   weight_substep(intermediate_timestep_count)
    1354        ENDDO
     1432       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
     1433
    13551434
    13561435    END SUBROUTINE advec_v_ws_ij
    … …  
    13721451       IMPLICIT NONE
    13731452
    1374        INTEGER ::  i, i_omp, j, k, tn
    1375        LOGICAL ::  degraded_l, degraded_s
    1376        REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
     1453       INTEGER ::  i, i_omp, j, k, tn, k_ppp, k_pp, k_mm
     1454       REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp, div
    13771455       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_r, diss_r, flux_n, &
    13781456                                       diss_n
    1379                                        
    1380        degraded_l = .FALSE. 
    1381        degraded_s = .FALSE.
    13821457
    13831458       gu = 2.0 * u_gtrans
    13841459       gv = 2.0 * v_gtrans
    13851460
    1386              
    1387        IF ( boundary_flags(j,i) /= 0 )  THEN
    1388 !       
    1389 !--      Degrade the order for Dirichlet bc. at the outflow boundary
    1390          SELECT CASE ( boundary_flags(j,i) )
    1391          
    1392             CASE ( 1 )
    1393                DO  k = nzb_w_inner(j,i)+1, nzt
    1394                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1395                   flux_r(k) = u_comp * (                                      &
    1396                               7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
    1397                             -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
    1398                   diss_r(k) = -ABS( u_comp ) * (                              &
    1399                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )               &
    1400                             -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
    1401                ENDDO
    1402                
    1403             CASE ( 2 )
    1404                DO  k = nzb_w_inner(j,i)+1, nzt
    1405                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1406                   flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25
    1407                   diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1), w(k,j,i),     &
    1408                                         w(k,j,i-1), w(k,j,i-2), u_comp,       &
    1409                                         0.25, ddx )
    1410                ENDDO
    1411                
    1412             CASE ( 3 )
    1413                DO  k = nzb_w_inner(j,i)+1, nzt
    1414                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1415                   flux_n(k) = v_comp * (                                      &
    1416                               7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
    1417                             -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    1418                   diss_n(k) = -ABS( v_comp ) * (                              &
    1419                               3.0 * ( w(k,j+1,i) - w(k,j,i)   )               &
    1420                             -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
    1421                ENDDO
    1422                
    1423             CASE ( 4 )
    1424                DO  k = nzb_w_inner(j,i)+1, nzt
    1425                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1426                   flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25
    1427                   diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i), w(k,j,i),     &
    1428                                         w(k,j-1,i), w(k,j-2,i), v_comp,       &
    1429                                         0.25, ddy )
    1430                ENDDO
    1431                
    1432             CASE ( 5 )
    1433                DO  k = nzb_w_inner(j,i)+1, nzt
    1434                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1435                   flux_r(k) = u_comp * (                                      &
    1436                               7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
    1437                             -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
    1438                   diss_r(k) = - ABS( u_comp ) * (                             &
    1439                               3.0 * ( w(k,j,i+1) - w(k,j,i) )                 &
    1440                             -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
    1441                ENDDO
    1442                
    1443             CASE ( 6 )
    1444                DO  k = nzb_w_inner(j,i)+1, nzt
    1445 !               
    1446 !--               Compute leftside fluxes for the left boundary of PE domain
    1447                   u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1448                   flux_r(k) = u_comp *(                                       &
    1449                               7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
    1450                             -       ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
    1451                   diss_r(k) = - ABS( u_comp ) * (                             &
    1452                               3.0 * ( w(k,j,i+1) - w(k,j,i) )                 &
    1453                             -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
    1454                  
    1455                   u_comp        = u(k+1,j,i) + u(k,j,i) - gu
    1456                   flux_l_w(k,j,tn) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
    1457                   diss_l_w(k,j,tn) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
    1458                                             w(k,j,i-1), w(k,j,i-1), u_comp,   &
    1459                                             0.25, ddx )
    1460                ENDDO
    1461                degraded_l = .TRUE.
    1462                
    1463             CASE ( 7 )
    1464                DO  k = nzb_w_inner(j,i)+1, nzt
    1465                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1466                   flux_n(k) = v_comp *(                                       &
    1467                               7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
    1468                             -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    1469                   diss_n(k) = - ABS( v_comp ) * (                             &
    1470                               3.0 * ( w(k,j+1,i) - w(k,j,i)   )               &
    1471                             -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
    1472                 ENDDO
    1473                
    1474             CASE ( 8 )
    1475                DO  k = nzb_w_inner(j,i)+1, nzt
    1476                   v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1477                   flux_n(k) = v_comp * (                                      &
    1478                               7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
    1479                             -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
    1480                   diss_n(k) = - ABS( v_comp ) * (                             &
    1481                               3.0 * ( w(k,j+1,i) - w(k,j,i) )                 &
    1482                             -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
    1483 !                   
    1484 !--              Compute southside fluxes for the south boundary of PE domain           
    1485                   v_comp      = v(k+1,j,i) + v(k,j,i) - gv
    1486                   flux_s_w(k,tn) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
    1487                   diss_s_w(k,tn) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
    1488                                           w(k,j-1,i), w(k,j-1,i), v_comp,     &
    1489                                           0.25, ddy )                 
    1490                ENDDO 
    1491                degraded_s = .TRUE.
    1492                
    1493             CASE DEFAULT
    1494            
    1495          END SELECT
    1496 !         
    1497 !--      Compute the crosswise 5th order fluxes at the outflow
    1498          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
    1499               boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
    1500          
    1501             DO  k = nzb_w_inner(j,i)+1, nzt
    1502                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1503                flux_n(k) = v_comp * (                                         &
    1504                            37.0 * ( w(k,j+1,i) + w(k,j,i)   )                 &
    1505                          -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                 &
    1506                          +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    1507                diss_n(k) = - ABS( v_comp ) * (                                &
    1508                            10.0 * ( w(k,j+1,i) - w(k,j,i)   )                 &
    1509                          -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                 &
    1510                          +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    1511             ENDDO
    1512            
    1513          ELSE
    1514          
    1515             DO  k = nzb_w_inner(j,i)+1, nzt         
    1516                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1517                flux_r(k) = u_comp * (                                         &
    1518                            37.0 * ( w(k,j,i+1) + w(k,j,i) )                   &
    1519                          -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                 &
    1520                          +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    1521                diss_r(k) = - ABS( u_comp ) * (                                &
    1522                            10.0 * ( w(k,j,i+1) - w(k,j,i)   )                 &
    1523                          -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                 &
    1524                          +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    1525             ENDDO
    1526            
    1527          ENDIF
    1528          
    1529        ELSE
    1530 !       
    1531 !--      Compute the fifth order fluxes for the interior of PE domain.               
    1532          DO  k = nzb_w_inner(j,i)+1, nzt
    1533             u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1534             flux_r(k) = u_comp * (                                            &
    1535                         37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
    1536                       -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
    1537                       +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    1538             diss_r(k) = - ABS( u_comp ) * (                                   &
    1539                         10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
    1540                       -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
    1541                       +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    1542                  
    1543             v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1544             flux_n(k) = v_comp * (                                            &
    1545                         37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
    1546                       -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
    1547                       +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    1548             diss_n(k) = - ABS( v_comp ) * (                                   &
    1549                         10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
    1550                       -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
    1551                       +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    1552          ENDDO
    1553          
    1554        ENDIF
    1555 !       
    1556 !--    Compute left- and southside fluxes for the respective boundary     
    1557        IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
    1558        
    1559           DO  k = nzb_w_inner(j,i)+1, nzt
     1461!
     1462!--    Compute southside fluxes for the respective boundary.
     1463       IF ( j == nys )  THEN
     1464
     1465          DO  k = nzb+1, nzb_max
     1466
    15601467             v_comp      = v(k+1,j,i) + v(k,j,i) - gv
    1561              flux_s_w(k,tn) = v_comp * (                                         &
     1468             flux_s_w(k,tn) = v_comp * (                                      &
     1469                         ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
     1470                      +     7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
     1471                      +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 &
     1472                         ) *                                                  &
     1473                                     ( w(k,j,i)   + w(k,j-1,i) )              &
     1474                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
     1475                      +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
     1476                         ) *                                                  &
     1477                                     ( w(k,j+1,i) + w(k,j-2,i) )              &
     1478                  +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
     1479                         ) *                                                  &
     1480                                     ( w(k,j+2,i) + w(k,j-3,i) )              &
     1481                                         )
     1482
     1483             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
     1484                         ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
     1485                      +     3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
     1486                      +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 &
     1487                         ) *                                                  &
     1488                                     ( w(k,j,i)   - w(k,j-1,i) )              &
     1489                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
     1490                      +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
     1491                         ) *                                                  &
     1492                                     ( w(k,j+1,i) - w(k,j-2,i) )              &
     1493                  +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
     1494                         ) *                                                  &
     1495                                     ( w(k,j+2,i) - w(k,j-3,i) )              &
     1496                                                 )
     1497
     1498          ENDDO
     1499
     1500          DO  k = nzb_max+1, nzt
     1501
     1502             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
     1503             flux_s_w(k,tn) = v_comp * (                                      &
    15621504                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
    15631505                         -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
    15641506                         +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    1565              diss_s_w(k,tn) = - ABS( v_comp ) * (                                &
     1507             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
    15661508                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
    15671509                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
    15681510                         +        ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
    1569           ENDDO
    1570          
     1511
     1512          ENDDO
     1513
    15711514       ENDIF
    1572        
    1573        IF ( i == i_omp  .AND.   .NOT. degraded_l ) THEN
    1574          
    1575           DO  k = nzb_w_inner(j,i)+1, nzt
    1576             u_comp        = u(k+1,j,i) + u(k,j,i) - gu
    1577             flux_l_w(k,j,tn) = u_comp * (                                        &
     1515!
     1516!--    Compute leftside fluxes for the respective boundary.
     1517       IF ( i == i_omp ) THEN
     1518
     1519          DO  k = nzb+1, nzb_max
     1520
     1521             u_comp         = u(k+1,j,i) + u(k,j,i) - gu
     1522             flux_l_w(k,j,tn) = u_comp * (                                     &
     1523                          ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
     1524                       +     7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
     1525                       +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 &
     1526                          ) *                                                  &
     1527                                     ( w(k,j,i)   + w(k,j,i-1) )               &
     1528                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
     1529                       +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
     1530                          ) *                                                  &
     1531                                     ( w(k,j,i+1) + w(k,j,i-2) )               &
     1532                   +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
     1533                          ) *                                                  &
     1534                                     ( w(k,j,i+2) + w(k,j,i-3) )               &
     1535                                         )
     1536
     1537               diss_l_w(k,j,tn) = - ABS( u_comp ) * (                          &
     1538                          ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
     1539                       +     3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
     1540                       +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 &
     1541                          ) *                                                  &
     1542                                     ( w(k,j,i)   - w(k,j,i-1) )               &
     1543                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
     1544                       +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
     1545                          ) *                                                  &
     1546                                     ( w(k,j,i+1) - w(k,j,i-2) )               &
     1547                   +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
     1548                          ) *                                                  &
     1549                                     ( w(k,j,i+2) - w(k,j,i-3) )               &
     1550                                                    )
     1551
     1552          ENDDO
     1553
     1554          DO  k = nzb_max+1, nzt
     1555
     1556             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
     1557             flux_l_w(k,j,tn) = u_comp * (                                    &
    15781558                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
    15791559                          -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
    15801560                          +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    1581             diss_l_w(k,j,tn) = - ABS( u_comp ) * (                               &
     1561             diss_l_w(k,j,tn) = - ABS( u_comp ) * (                           &
    15821562                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
    15831563                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
    15841564                          +        ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
    1585           ENDDO
    1586          
     1565
     1566          ENDDO
     1567
    15871568       ENDIF
    15881569
    1589 !       
    1590 !--    Now compute the tendency terms for the horizontal parts.         
    1591        DO  k = nzb_w_inner(j,i)+1, nzt
     1570       flux_t(0) = 0.0
     1571       diss_t(0) = 0.0
     1572       flux_d    = 0.0
     1573       diss_d    = 0.0
     1574!
     1575!--    Now compute the fluxes and tendency terms for the horizontal
     1576!--    and vertical parts.
     1577       DO  k = nzb+1, nzb_max
     1578
     1579          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
     1580          flux_r(k) = u_comp * (                                             &
     1581                     ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
     1582                  +     7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
     1583                  +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1    &
     1584                     ) *                                                     &
     1585                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
     1586              -      (  8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
     1587                  +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
     1588                     ) *                                                     &
     1589                                 ( w(k,j,i+2) + w(k,j,i-1) )                 &
     1590              +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
     1591                     ) *                                                     &
     1592                                 ( w(k,j,i+3) + w(k,j,i-2) )                 &
     1593                               )
     1594
     1595          diss_r(k) = - ABS( u_comp ) * (                                    &
     1596                     ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
     1597                  +     3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
     1598                  +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1    &
     1599                     ) *                                                     &
     1600                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
     1601              -      (  5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
     1602                  +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
     1603                     ) *                                                     &
     1604                                 ( w(k,j,i+2) - w(k,j,i-1) )                 &
     1605              +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
     1606                     ) *                                                     &
     1607                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
     1608                                        )
     1609
     1610          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
     1611          flux_n(k) = v_comp * (                                             &
     1612                     ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
     1613                  +     7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
     1614                  +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1    &
     1615                     ) *                                                     &
     1616                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
     1617              -      (  8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
     1618                  +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
     1619                     ) *                                                     &
     1620                                 ( w(k,j+2,i) + w(k,j-1,i) )                 &
     1621              +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
     1622                     ) *                                                     &
     1623                                 ( w(k,j+3,i) + w(k,j-2,i) )                 &
     1624                               )
     1625
     1626          diss_n(k) = - ABS( v_comp ) * (                                    &
     1627                     ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
     1628                  +     3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
     1629                  +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1    &
     1630                     ) *                                                     &
     1631                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
     1632              -      (  5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
     1633                  +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
     1634                     ) *                                                     &
     1635                                 ( w(k,j+2,i) - w(k,j-1,i) )                 &
     1636              +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
     1637                     ) *                                                     &
     1638                                 ( w(k,j+3,i) - w(k,j-2,i) )                 &
     1639                                        )
     1640!
     1641!--       k index has to be modified near bottom and top, else array
     1642!--       subscripts will be exceeded.
     1643          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
     1644          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1)  )
     1645          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
     1646
     1647          w_comp    = w(k+1,j,i) + w(k,j,i)
     1648          flux_t(k) = w_comp  * (                                            &
     1649                     ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1650                  +     7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1651                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
     1652                     ) *                                                     &
     1653                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
     1654              -      (  8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1655                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1656                     ) *                                                     &
     1657                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
     1658              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1659                     ) *                                                     &
     1660                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
     1661                                 )
     1662
     1663          diss_t(k) = - ABS( w_comp ) * (                                    &
     1664                     ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1665                  +     3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1666                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
     1667                     ) *                                                     &
     1668                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
     1669              -      (  5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1670                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1671                     ) *                                                     &
     1672                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
     1673              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1674                     ) *                                                     &
     1675                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
     1676                                         )
     1677!
     1678!--       Calculate the divergence of the velocity field. A respective
     1679!--       correction is needed to overcome numerical instabilities introduced
     1680!--       by a not sufficient reduction of divergences near topography.
     1681          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
     1682              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
     1683              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
     1684                ) * 0.5
     1685
    15921686          tend(k,j,i) = tend(k,j,i) - (                                       &
    15931687                      ( flux_r(k) + diss_r(k)                                 &
    1594                     -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx               &
     1688                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
    15951689                    + ( flux_n(k) + diss_n(k)                                 &
    1596                     -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy  )
     1690                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
     1691                    + ( flux_t(k) + diss_t(k)                                 &
     1692                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
     1693                                      ) + div * w(k,j,i)
    15971694
    15981695          flux_l_w(k,j,tn) = flux_r(k)
    15991696          diss_l_w(k,j,tn) = diss_r(k)
    1600           flux_s_w(k,tn) = flux_n(k)
    1601           diss_s_w(k,tn) = diss_n(k)
     1697          flux_s_w(k,tn)   = flux_n(k)
     1698          diss_s_w(k,tn)   = diss_n(k)
     1699          flux_d           = flux_t(k)
     1700          diss_d           = diss_t(k)
     1701!
     1702!--        Statistical Evaluation of w'w'.
     1703          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
     1704                             + ( flux_t(k) + diss_t(k) )                     &
     1705                             *   weight_substep(intermediate_timestep_count)
     1706
    16021707       ENDDO
    16031708
    1604 !
    1605 !--    Vertical advection, degradation of order near surface and top.
    1606 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    1607 !--    statistical evaluation the top flux at the surface should be 0
    1608        flux_t(nzb_w_inner(j,i)) = 0.0  !statistical reasons
    1609        diss_t(nzb_w_inner(j,i)) = 0.0
    1610 !
    1611 !--    2nd order scheme (bottom)       
    1612        k         = nzb_w_inner(j,i)+1
    1613        flux_d    = flux_t(k-1)
    1614        diss_d    = diss_t(k-1)
    1615        w_comp    = w(k+1,j,i) + w(k,j,i)
    1616        flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
    1617        diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0,        &
    1618                              w_comp, 0.25, ddzu(k+1) )
    1619                      
    1620        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1621                                  -   flux_d    - diss_d      ) * ddzu(k+1)               
    1622 !
    1623 !--    WS3 as an intermediate step (bottom)
    1624        k         = nzb_w_inner(j,i)+2
    1625        flux_d    = flux_t(k-1)
    1626        diss_d    = diss_t(k-1)
    1627        w_comp    = w(k+1,j,i) + w(k,j,i)
    1628        flux_t(k) = w_comp * (                                                 &
    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
    1634 
    1635        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1636                                  -   flux_d    - diss_d      ) * ddzu(k+1)
    1637 !                                 
    1638 !--    WS5
    1639        DO  k = nzb_w_inner(j,i)+3, nzt-2
    1640           flux_d    = flux_t(k-1)
    1641           diss_d    = diss_t(k-1)
     1709       DO  k = nzb_max+1, nzt
     1710
     1711          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
     1712          flux_r(k) = u_comp * (                                            &
     1713                      37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
     1714                    -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
     1715                    +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
     1716
     1717          diss_r(k) = - ABS( u_comp ) * (                                   &
     1718                      10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
     1719                    -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
     1720                    +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     1721
     1722          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
     1723          flux_n(k) = v_comp * (                                            &
     1724                      37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
     1725                    -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
     1726                    +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
     1727
     1728          diss_n(k) = - ABS( v_comp ) * (                                   &
     1729                      10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
     1730                    -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
     1731                    +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     1732!
     1733!--       k index has to be modified near bottom and top, else array
     1734!--       subscripts will be exceeded.
     1735          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
     1736          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1)  )
     1737          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
     1738
    16421739          w_comp    = w(k+1,j,i) + w(k,j,i)
    1643           flux_t(k) = w_comp * (                                              &
    1644                       37.0 * ( w(k+1,j,i) + w(k,j,i)   )                      &
    1645                     -  8.0 * ( w(k+2,j,i) + w(k-1,j,i) )                      &
    1646                     +        ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
    1647           diss_t(k) = - ABS( w_comp ) * (                                     &
    1648                       10.0 * ( w(k+1,j,i) - w(k,j,i)   )                      &
    1649                     -  5.0 * ( w(k+2,j,i) - w(k-1,j,i) )                      &
    1650                     +        ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
    1651 
    1652           tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    1653                                     -   flux_d    - diss_d      ) * ddzu(k+1)
     1740          flux_t(k) = w_comp  * (                                            &
     1741                     ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1742                  +     7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1743                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
     1744                     ) *                                                     &
     1745                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
     1746              -      (  8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1747                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1748                     ) *                                                     &
     1749                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
     1750              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1751                     ) *                                                     &
     1752                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
     1753                                 )
     1754
     1755          diss_t(k) = - ABS( w_comp ) * (                                    &
     1756                     ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1757                  +     3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1758                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
     1759                     ) *                                                     &
     1760                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
     1761              -      (  5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1762                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
     1763                     ) *                                                     &
     1764                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
     1765              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
     1766                     ) *                                                     &
     1767                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
     1768                                         )
     1769!
     1770!--       Calculate the divergence of the velocity field. A respective
     1771!--       correction is needed to overcome numerical instabilities introduced
     1772!--       by a not sufficient reduction of divergences near topography.
     1773          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
     1774              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
     1775              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
     1776                ) * 0.5
     1777
     1778          tend(k,j,i) = tend(k,j,i) - (                                       &
     1779                      ( flux_r(k) + diss_r(k)                                 &
     1780                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
     1781                    + ( flux_n(k) + diss_n(k)                                 &
     1782                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
     1783                    + ( flux_t(k) + diss_t(k)                                 &
     1784                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
     1785                                      ) + div * w(k,j,i)
     1786
     1787          flux_l_w(k,j,tn) = flux_r(k)
     1788          diss_l_w(k,j,tn) = diss_r(k)
     1789          flux_s_w(k,tn)   = flux_n(k)
     1790          diss_s_w(k,tn)   = diss_n(k)
     1791          flux_d           = flux_t(k)
     1792          diss_d           = diss_t(k)
     1793!
     1794!--        Statistical Evaluation of w'w'.
     1795          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
     1796                             + ( flux_t(k) + diss_t(k) )                     &
     1797                             *   weight_substep(intermediate_timestep_count)
     1798
    16541799       ENDDO
    1655 !--    WS3 as an intermediate step (top)
    1656        k         = nzt - 1
    1657        flux_d    = flux_t(k-1)
    1658        diss_d    = diss_t(k-1)
    1659        w_comp    = w(k+1,j,i) + w(k,j,i)
    1660        flux_t(k) = w_comp * (                                                 &
    1661                    7.0 * ( w(k+1,j,i) + w(k,j,i)   )                          &
    1662                  -       ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
    1663        diss_t(k) = - ABS( w_comp ) * (                                        &
    1664                    3.0 * ( w(k+1,j,i) - w(k,j,i)     )                        &
    1665                  -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
    1666                    
    1667        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1668                                  -   flux_d    - diss_d      ) * ddzu(k+1)
    1669 !
    1670 !--    2nd order scheme (top)
    1671        k         = nzt
    1672        flux_d    = flux_t(k-1)
    1673        diss_d    = diss_t(k-1)
    1674        w_comp    = w(k+1,j,i) + w(k,j,i)
    1675        flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
    1676        diss_t(k) = diss_2nd( w(k+1,j,i), w(k+1,j,i), w(k,j,i), w(k-1,j,i),    &
    1677                              w(k-2,j,i), w_comp, 0.25, ddzu(k+1) )
    1678 
    1679        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    1680                                  -   flux_d    - diss_d      ) * ddzu(k+1)
    1681        
    1682        DO  k = nzb_w_inner(j,i), nzt
    1683            sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                               &
    1684                  + ( flux_t(k) + diss_t(k) )                                  &
    1685                  *   weight_substep(intermediate_timestep_count)
    1686        ENDDO
     1800
    16871801
    16881802    END SUBROUTINE advec_w_ws_ij
    … …  
    17031817       IMPLICIT NONE
    17041818
    1705        INTEGER ::  i, j, k, tn = 0
     1819       INTEGER ::  i, j, k, tn = 0, k_ppp, k_pp, k_mm
    17061820       REAL, DIMENSION(:,:,:), POINTER ::  sk
    1707        REAL    :: flux_d, diss_d, u_comp, v_comp
    1708        REAL, DIMENSION(nzb:nzt+1)  ::  flux_r, diss_r, flux_n, diss_n
    1709        REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local,    &
    1710                                      flux_t, diss_t
     1821       REAL    :: flux_d, diss_d, u_comp, v_comp, div
     1822       REAL, DIMENSION(nzb:nzt)  ::  flux_r, diss_r, flux_n, diss_n, flux_t,  &
     1823                                     diss_t
     1824       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local
    17111825       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
    17121826                                             swap_diss_x_local
    … …  
    17171831       i = nxl
    17181832       DO  j = nys, nyn
    1719           IF ( boundary_flags(j,i) == 6 )  THEN
    1720          
    1721              DO  k = nzb_s_inner(j,i)+1, nzt
    1722                 u_comp                 = u(k,j,i) - u_gtrans
    1723                 swap_flux_x_local(k,j) = u_comp * (                           &
    1724                                          sk(k,j,i) + sk(k,j,i-1)) * 0.5
    1725                 swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2), sk(k,j,i+1),  &
    1726                                                    sk(k,j,i), sk(k,j,i-1),    &
    1727                                                    sk(k,j,i-1), u_comp,       &
    1728                                                    0.5, ddx )
    1729              ENDDO
    1730              
    1731           ELSE
    1732          
    1733              DO  k = nzb_s_inner(j,i)+1, nzt
    1734                 u_comp                 = u(k,j,i) - u_gtrans
    1735                 swap_flux_x_local(k,j) = u_comp*(                              &
    1736                                          37.0 * ( sk(k,j,i)+sk(k,j,i-1)     )  &
    1737                                        -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
    1738                                        +        ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
    1739                                        * adv_sca_5
    1740                 swap_diss_x_local(k,j) = - ABS( u_comp ) * (                   &
    1741                                          10.0 * (sk(k,j,i) - sk(k,j,i-1)    )  &
    1742                                        -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
    1743                                        +        ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
    1744                                        * adv_sca_5
    1745              ENDDO
    1746           ENDIF
     1833
     1834          DO  k = nzb+1, nzb_max
     1835
     1836             u_comp                 = u(k,j,i) - u_gtrans
     1837             swap_flux_x_local(k,j) = u_comp * (                              &
     1838                         ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     1839                      +     7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     1840                      +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1  &
     1841                         ) *                                                  &
     1842                                     ( sk(k,j,i)   + sk(k,j,i-1)    )         &
     1843                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     1844                      +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     1845                         ) *                                                  &
     1846                                     ( sk(k,j,i+1) + sk(k,j,i-2)    )         &
     1847                  +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     1848                         ) *         ( sk(k,j,i+2) + sk(k,j,i-3)    )         &
     1849                                                  )
     1850
     1851              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
     1852                         ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     1853                      +     3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     1854                      +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1  &
     1855                         ) *                                                  &
     1856                                     ( sk(k,j,i)   - sk(k,j,i-1)    )         &
     1857                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     1858                      +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
     1859                         ) *                                                  &
     1860                                     ( sk(k,j,i+1) - sk(k,j,i-2)  )           &
     1861                  +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
     1862                         ) *                                                  &
     1863                                     ( sk(k,j,i+2) - sk(k,j,i-3) )            &
     1864                                                           )
     1865
     1866          ENDDO
     1867
     1868          DO  k = nzb_max+1, nzt
     1869
     1870             u_comp                 = u(k,j,i) - u_gtrans
     1871             swap_flux_x_local(k,j) = u_comp * (                             &
     1872                                      37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
     1873                                    -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
     1874                                    +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
     1875                                                  ) * adv_sca_5
     1876
     1877             swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
     1878                                      10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
     1879                                    -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
     1880                                    +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
     1881                                                          ) * adv_sca_5
     1882
     1883          ENDDO
     1884
    17471885       ENDDO
    1748 !
    1749 !--    The following loop computes the horizontal fluxes for the interior of the
    1750 !--    processor domain plus south boundary points. Furthermore tendency terms
    1751 !--    are computed.
     1886
    17521887       DO  i = nxl, nxr
     1888
    17531889          j = nys
    1754           IF ( boundary_flags(j,i) == 8 )  THEN
    1755          
    1756              DO  k = nzb_s_inner(j,i)+1, nzt
    1757                 v_comp               = v(k,j,i) - v_gtrans
    1758                 swap_flux_y_local(k) = v_comp *                                &
    1759                                        ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
    1760                 swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),     &
    1761                                                  sk(k,j,i), sk(k,j-1,i),       &
    1762                                                  sk(k,j-1,i), v_comp, 0.5, ddy )
    1763              ENDDO
    1764              
    1765           ELSE
    1766          
    1767              DO  k = nzb_s_inner(j,i)+1, nzt
    1768                 v_comp               = v(k,j,i) - v_gtrans
    1769                 swap_flux_y_local(k) = v_comp * (                              &
    1770                                        37.0 * ( sk(k,j,i) + sk(k,j-1,i)   )    &
    1771                                      -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
    1772                                      +        ( sk(k,j+2,i) + sk(k,j-3,i) ) )  & 
    1773                                      * adv_sca_5
    1774                 swap_diss_y_local(k)= - ABS( v_comp ) * (                      &
    1775                                       10.0 * ( sk(k,j,i) - sk(k,j-1,i)   )     &
    1776                                     -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
    1777                                     +        ( sk(k,j+2,i)-sk(k,j-3,i) ) )     &
    1778                                     * adv_sca_5
    1779              ENDDO
    1780              
    1781           ENDIF
     1890          DO  k = nzb+1, nzb_max
     1891
     1892             v_comp                  = v(k,j,i) - v_gtrans
     1893             swap_flux_y_local(k) = v_comp * (                                &
     1894                         ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     1895                      +     7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     1896                      +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1  &
     1897                         ) *                                                  &
     1898                                     ( sk(k,j,i)  + sk(k,j-1,i)     )         &
     1899                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     1900                      +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     1901                         ) *                                                  &
     1902                                     ( sk(k,j+1,i) + sk(k,j-2,i)    )         &
     1903                  +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     1904                         ) *                                                  &
     1905                                     ( sk(k,j+2,i) + sk(k,j-3,i)    )         &
     1906                                             )
     1907
     1908             swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
     1909                         ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     1910                      +     3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     1911                      +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1  &
     1912                         ) *                                                  &
     1913                                     ( sk(k,j,i)   - sk(k,j-1,i)    )         &
     1914                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     1915                      +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
     1916                         ) *                                                  &
     1917                                     ( sk(k,j+1,i) - sk(k,j-2,i)  )           &
     1918                  +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
     1919                         ) *                                                  &
     1920                                     ( sk(k,j+2,i) - sk(k,j-3,i) )            &
     1921                                                     )
     1922
     1923          ENDDO
     1924!
     1925!--       Above to the top of the highest topography. No degradation necessary.
     1926          DO  k = nzb_max+1, nzt
     1927
     1928             v_comp               = v(k,j,i) - v_gtrans
     1929             swap_flux_y_local(k) = v_comp * (                               &
     1930                                    37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
     1931                                  -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
     1932                                  +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
     1933                                             ) * adv_sca_5
     1934              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
     1935                                    10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
     1936                                  -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
     1937                                  +          sk(k,j+2,i) - sk(k,j-3,i)       &
     1938                                                      ) * adv_sca_5
     1939
     1940          ENDDO
    17821941
    17831942          DO  j = nys, nyn
    1784             IF ( boundary_flags(j,i) /= 0 )  THEN
    1785 !
    1786 !--            Degrade the order for Dirichlet bc. at the outflow boundary
    1787                SELECT CASE ( boundary_flags(j,i) )
    1788                
    1789                   CASE ( 1 )
    1790                      DO  k = nzb_s_inner(j,i)+1, nzt
    1791                         u_comp    = u(k,j,i+1) - u_gtrans
    1792                         flux_r(k) = u_comp * (                                 &
    1793                                     7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
    1794                                   -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
    1795                                   * adv_sca_3
    1796                         diss_r(k) = - ABS( u_comp ) * (                        &
    1797                                     3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        &
    1798                                   -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
    1799                                   * adv_sca_3
    1800                      ENDDO
    1801                      
    1802                   CASE ( 2 )
    1803                      DO  k = nzb_s_inner(j,i)+1, nzt
    1804                         u_comp    = u(k,j,i+1) - u_gtrans
    1805                         flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
    1806                         diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),        &
    1807                                               sk(k,j,i), sk(k,j,i-1),          &
    1808                                               sk(k,j,i-2), u_comp, 0.5, ddx )
    1809                      ENDDO
    1810                      
    1811                   CASE ( 3 )
    1812                      DO  k = nzb_s_inner(j,i)+1, nzt
    1813                         v_comp    = v(k,j+1,i) - v_gtrans
    1814                         flux_n(k) = v_comp * (                                 &
    1815                                     7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
    1816                                   -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
    1817                                   * adv_sca_3
    1818                         diss_n(k) = - ABS( v_comp ) * (                        &
    1819                                     3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        &
    1820                                   -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
    1821                                   * adv_sca_3
    1822                      ENDDO
    1823                      
    1824                   CASE ( 4 )
    1825                      DO  k = nzb_s_inner(j,i)+1, nzt
    1826                         v_comp    = v(k,j+1,i) - v_gtrans
    1827                         flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
    1828                         diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),        &
    1829                                               sk(k,j,i), sk(k,j-1,i),          &
    1830                                               sk(k,j-2,i), v_comp, 0.5, ddy )     
    1831                      ENDDO
    1832                      
    1833                   CASE ( 5 )
    1834                      DO  k = nzb_w_inner(j,i)+1, nzt
    1835                         u_comp    = u(k,j,i+1) - u_gtrans
    1836                         flux_r(k) = u_comp * (                                 &
    1837                                     7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
    1838                                   -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
    1839                                   * adv_sca_3
    1840                         diss_r(k) = - ABS( u_comp ) * (                        &
    1841                                     3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        &
    1842                                   -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
    1843                                   * adv_sca_3
    1844                      ENDDO 
    1845                      
    1846                   CASE ( 6 )
    1847                      DO  k = nzb_s_inner(j,i)+1, nzt
    1848                         u_comp    = u(k,j,i+1) - u_gtrans
    1849                         flux_r(k) = u_comp * (                                 &
    1850                                     7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
    1851                                   -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
    1852                                   * adv_sca_3
    1853                         diss_r(k) = - ABS( u_comp ) * (                        &
    1854                                     3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        &
    1855                                   -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
    1856                                   * adv_sca_3
    1857                      ENDDO
    1858                      
    1859                   CASE ( 7 )
    1860                      DO  k = nzb_s_inner(j,i)+1, nzt
    1861                         v_comp    = v(k,j+1,i) - v_gtrans
    1862                         flux_n(k) = v_comp * (                                 &
    1863                                     7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
    1864                                   -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
    1865                                   * adv_sca_3
    1866                         diss_n(k) = - ABS( v_comp ) * (                        &
    1867                                     3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        &
    1868                                   -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
    1869                                   * adv_sca_3
    1870                      ENDDO
    1871                      
    1872                   CASE ( 8 )
    1873                      DO  k = nzb_s_inner(j,i)+1, nzt
    1874                         v_comp    = v(k,j+1,i) - v_gtrans
    1875                         flux_n(k) = v_comp * (                                 &
    1876                                     7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
    1877                                   -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
    1878                                   * adv_sca_3
    1879                         diss_n(k) = -  ABS( v_comp ) * (                       &
    1880                                     3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        &
    1881                                   -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
    1882                                   * adv_sca_3
    1883                      ENDDO 
    1884                      
    1885                   CASE DEFAULT
    1886                  
    1887                END SELECT
    1888                
    1889                IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
    1890                    boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )    &
    1891                THEN
    1892              
    1893                   DO  k = nzb_s_inner(j,i)+1, nzt
    1894                      v_comp    = v(k,j+1,i) - v_gtrans
    1895                      flux_n(k) = v_comp * (                                    &
    1896                                  37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )          &
    1897                                -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )          &
    1898                                +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )        &
    1899                                * adv_sca_5
    1900                      diss_n(k) = - ABS( v_comp ) * (                           &
    1901                                  10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )          &
    1902                                -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )          &
    1903                                +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )        &
    1904                                * adv_sca_5
    1905                   ENDDO
    1906                  
    1907                ELSE
    1908                
    1909                   DO  k = nzb_s_inner(j,i)+1, nzt
    1910                      u_comp    = u(k,j,i+1) - u_gtrans 
    1911                      flux_r(k) = u_comp * (                                    &
    1912                                  37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )          &
    1913                                -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )          &
    1914                                +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )        &
    1915                                * adv_sca_5
    1916                      diss_r(k) = - ABS( u_comp ) * (                           &
    1917                                  10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )          &
    1918                                -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )          &
    1919                                +        ( sk(k,j,i+3) - sk(k,j,i-2) ) )        &
    1920                                * adv_sca_5
    1921                   ENDDO
    1922                  
    1923                ENDIF     
    1924            
    1925             ELSE
    1926            
    1927                DO  k = nzb_s_inner(j,i)+1, nzt
    1928                   u_comp    = u(k,j,i+1) - u_gtrans
    1929                   flux_r(k) = u_comp * (                                       &
    1930                               37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
    1931                             -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )             &
    1932                             +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )           &
    1933                             * adv_sca_5
    1934                   diss_r(k) = - ABS( u_comp ) * (                              &
    1935                               10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
    1936                             -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )             &
    1937                             +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )             &
    1938                             * adv_sca_5
    1939                  
    1940                   v_comp    = v(k,j+1,i) - v_gtrans
    1941                   flux_n(k) = v_comp * (                                       &
    1942                               37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
    1943                             -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )             &
    1944                             +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )           &
    1945                             * adv_sca_5
    1946                   diss_n(k) = - ABS( v_comp ) * (                              &
    1947                               10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
    1948                             -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )             &
    1949                             +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )           &
    1950                             * adv_sca_5                 
    1951                ENDDO
    1952                
    1953             ENDIF
    1954            
    1955             DO  k = nzb_s_inner(j,i)+1, nzt
    1956                tend(k,j,i) = tend(k,j,i) - (                                  &
    1957                   ( flux_r(k) + diss_r(k)                                     &
    1958                 -   swap_flux_x_local(k,j) - swap_diss_x_local(k,j)   ) * ddx &
    1959                 + ( flux_n(k) + diss_n(k)                                     &
    1960                 -   swap_flux_y_local(k) - swap_diss_y_local(k)       ) * ddy)
    1961                    
     1943
     1944             flux_t(0) = 0.0
     1945             diss_t(0) = 0.0
     1946             flux_d    = 0.0
     1947             diss_d    = 0.0
     1948
     1949             DO  k = nzb+1, nzb_max
     1950
     1951                u_comp    = u(k,j,i+1) - u_gtrans
     1952                flux_r(k) = u_comp * (                                      &
     1953                     ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     1954                  +     7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     1955                  +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1    &
     1956                     ) *                                                    &
     1957                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
     1958              -      (  8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     1959                  +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     1960                     ) *                                                    &
     1961                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
     1962              +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     1963                     ) *     ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
     1964                                     )
     1965
     1966                diss_r(k) = -ABS( u_comp ) * (                              &
     1967                     ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     1968                  +     3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     1969                  +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1    &
     1970                     ) *                                                    &
     1971                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
     1972              -      (  5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     1973                  +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
     1974                     ) *                                                    &
     1975                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
     1976              +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
     1977                     ) *                                                    &
     1978                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
     1979                                             )
     1980
     1981                v_comp    = v(k,j+1,i) - v_gtrans
     1982                flux_n(k) = v_comp * (                                      &
     1983                     ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     1984                  +     7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     1985                  +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1    &
     1986                     ) *                                                    &
     1987                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
     1988              -      (  8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     1989                  +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     1990                     ) *                                                    &
     1991                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
     1992              +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     1993                     ) *                                                    &
     1994                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
     1995                                     )
     1996
     1997                diss_n(k) = -ABS( v_comp ) * (                              &
     1998                     ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     1999                  +     3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     2000                  +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1    &
     2001                     ) *                                                    &
     2002                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
     2003              -      (  5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     2004                  +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
     2005                     ) *                                                    &
     2006                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
     2007              +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
     2008                     ) *                                                    &
     2009                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
     2010                                             )
     2011!
     2012!--             k index has to be modified near bottom and top, else array
     2013!--             subscripts will be exceeded.
     2014                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
     2015                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1)  )
     2016                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
     2017
     2018
     2019                flux_t(k) = w(k,j,i) * (                                    &
     2020                     ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2021                  +     7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2022                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     2023                     ) *                                                    &
     2024                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
     2025              -      (  8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2026                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2027                     ) *                                                    &
     2028                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
     2029              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2030                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
     2031                                       )
     2032
     2033                diss_t(k) = -ABS( w(k,j,i) ) * (                            &
     2034                     ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2035                  +     3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2036                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     2037                     ) *                                                    &
     2038                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
     2039              -      (  5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2040                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2041                     ) *                                                    &
     2042                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
     2043              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2044                     ) *                                                    &
     2045                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
     2046                                               )
     2047!
     2048!--             Calculate the divergence of the velocity field. A respective
     2049!--             correction is needed to overcome numerical instabilities caused
     2050!--             by a not sufficient reduction of divergences near topography.
     2051                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
     2052                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
     2053                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     2054
     2055                tend(k,j,i) = tend(k,j,i) - (                                 &
     2056                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
     2057                          swap_diss_x_local(k,j)            ) * ddx           &
     2058                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
     2059                          swap_diss_y_local(k)              ) * ddy           &
     2060                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
     2061                                                               ) * ddzw(k)    &
     2062                                            ) + sk(k,j,i) * div
     2063
     2064                swap_flux_y_local(k)   = flux_n(k)
     2065                swap_diss_y_local(k)   = diss_n(k)
    19622066                swap_flux_x_local(k,j) = flux_r(k)
    19632067                swap_diss_x_local(k,j) = diss_r(k)
     2068                flux_d                 = flux_t(k)
     2069                diss_d                 = diss_t(k)
     2070
     2071             ENDDO
     2072
     2073             DO  k = nzb_max+1, nzt
     2074
     2075                u_comp    = u(k,j,i+1) - u_gtrans
     2076                flux_r(k) = u_comp * (                                      &
     2077                      37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                  &
     2078                    -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
     2079                    +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     2080                diss_r(k) = -ABS( u_comp ) * (                              &
     2081                      10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                  &
     2082                    -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
     2083                    +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     2084
     2085                v_comp    = v(k,j+1,i) - v_gtrans
     2086                flux_n(k) = v_comp * (                                      &
     2087                      37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                  &
     2088                    -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
     2089                    +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     2090                diss_n(k) = -ABS( v_comp ) * (                              &
     2091                      10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                  &
     2092                    -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                  &
     2093                    +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     2094!
     2095!--             k index has to be modified near bottom and top, else array
     2096!--             subscripts will be exceeded.
     2097                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
     2098                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1)  )
     2099                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
     2100
     2101                flux_t(k) = w(k,j,i) * (                                    &
     2102                     ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2103                  +     7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2104                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     2105                     ) *                                                    &
     2106                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
     2107              -      (  8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2108                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2109                     ) *                                                    &
     2110                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
     2111              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2112                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
     2113                                       )
     2114
     2115                diss_t(k) = -ABS( w(k,j,i) ) * (                            &
     2116                     ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2117                  +     3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2118                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
     2119                     ) *                                                    &
     2120                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
     2121              -      (  5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2122                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
     2123                     ) *                                                    &
     2124                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
     2125              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
     2126                     ) *                                                    &
     2127                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
     2128                                               )
     2129!
     2130!--             Calculate the divergence of the velocity field. A respective
     2131!--             correction is needed to overcome numerical instabilities introduced
     2132!--             by a not sufficient reduction of divergences near topography.
     2133                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
     2134                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
     2135                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     2136
     2137                tend(k,j,i) = tend(k,j,i) - (                                 &
     2138                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
     2139                          swap_diss_x_local(k,j)            ) * ddx           &
     2140                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
     2141                          swap_diss_y_local(k)              ) * ddy           &
     2142                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
     2143                                                               ) * ddzw(k)    &
     2144                                            ) + sk(k,j,i) * div
     2145
    19642146                swap_flux_y_local(k)   = flux_n(k)
    19652147                swap_diss_y_local(k)   = diss_n(k)
    1966             ENDDO
     2148                swap_flux_x_local(k,j) = flux_r(k)
     2149                swap_diss_x_local(k,j) = diss_r(k)
     2150                flux_d                 = flux_t(k)
     2151                diss_d                 = diss_t(k)
     2152
     2153             ENDDO
     2154!
     2155!--          evaluation of statistics
     2156             SELECT CASE ( sk_char )
     2157
     2158                 CASE ( 'pt' )
     2159                    DO  k = nzb, nzt
     2160                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
     2161                        + ( flux_t(k) + diss_t(k) )                          &
     2162                        *   weight_substep(intermediate_timestep_count)
     2163                    ENDDO
     2164                 CASE ( 'sa' )
     2165                    DO  k = nzb, nzt
     2166                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
     2167                        + ( flux_t(k) + diss_t(k) )                          &
     2168                        *   weight_substep(intermediate_timestep_count)
     2169                    ENDDO
     2170                 CASE ( 'q' )
     2171                    DO  k = nzb, nzt
     2172                       sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
     2173                        + ( flux_t(k) + diss_t(k) )                          &
     2174                        *   weight_substep(intermediate_timestep_count)
     2175                    ENDDO
     2176
     2177              END SELECT
     2178
    19672179         ENDDO
    19682180      ENDDO
    1969 
    1970 !
    1971 !--   Vertical advection, degradation of order near surface and top.
    1972 !--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    1973 !--   statistical evaluation the top flux at the surface should be 0
    1974       DO  i = nxl, nxr
    1975          DO  j = nys, nyn
    1976 !
    1977 !--         2nd order scheme (bottom)
    1978             k=nzb_s_inner(j,i)+1
    1979 !           
    1980 !--         The fluxes flux_d and diss_d at the surface are 0. Due to static
    1981 !--         reasons the top flux at the surface should be 0.
    1982             flux_t(nzb_s_inner(j,i)) = 0.0
    1983             diss_t(nzb_s_inner(j,i)) = 0.0
    1984             flux_d = flux_t(k-1)
    1985             diss_d = diss_t(k-1)
    1986             flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
    1987             diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i),        &
    1988                                   sk(k,j,i), sk(k,j,i), w(k,j,i),             &
    1989                                   0.5, ddzw(k) )
    1990             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    1991                                       -   flux_d    - diss_d       ) * ddzw(k)   
    1992 !
    1993 !--         WS3 as an intermediate step (bottom)
    1994             k         = nzb_s_inner(j,i)+2
    1995             flux_d    = flux_t(k-1)
    1996             diss_d    = diss_t(k-1)
    1997             flux_t(k) = w(k,j,i) * (                                          &
    1998                         7.0 * ( sk(k+1,j,i) + sk(k,j,i) )                     &
    1999                        - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    2000             diss_t(k) = - ABS( w(k,j,i) ) * (                                 &
    2001                         3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   &
    2002                       -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    2003             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2004                                       -   flux_d    - diss_d       ) * ddzw(k)
    2005 !
    2006 !--         WS5
    2007             DO  k = nzb_s_inner(j,i)+3, nzt-2
    2008                flux_d    = flux_t(k-1)
    2009                diss_d    = diss_t(k-1)
    2010                flux_t(k) = w(k,j,i) * (                                       &
    2011                            37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )               &
    2012                          -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
    2013                          +        ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
    2014                diss_t(k) = - ABS(w(k,j,i)) * (                                &
    2015                            10.0 * ( sk(k+1,j,i) -sk(k,j,i)    )               &
    2016                          -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
    2017                          +        ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
    2018 
    2019                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
    2020                                          -   flux_d    - diss_d     ) * ddzw(k)
    2021             ENDDO
    2022 !
    2023 !--         WS3 as an intermediate step (top)
    2024             k         = nzt - 1
    2025             flux_d    = flux_t(k-1)
    2026             diss_d    = diss_t(k-1)
    2027             flux_t(k) = w(k,j,i) * (                                          &
    2028                         7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )                   &
    2029                       -       ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    2030             diss_t(k) = - ABS(w(k,j,i)) * (                                   &
    2031                         3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   &
    2032                       -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    2033                        
    2034             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2035                                       -   flux_d    - diss_d       ) * ddzw(k)
    2036 !
    2037 !--         2nd order scheme (top)
    2038             k         = nzt
    2039             flux_d    = flux_t(k-1)
    2040             diss_d    = diss_t(k-1)
    2041             flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
    2042             diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i),        &
    2043                                   sk(k-1,j,i), sk(k-2,j,i), w(k,j,i),         &
    2044                                   0.5, ddzw(k) )
    2045 
    2046             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2047                                       -   flux_d    - diss_d       ) * ddzw(k)
    2048 !
    2049 !--         evaluation of statistics
    2050             SELECT CASE ( sk_char )
    2051 
    2052                CASE ( 'pt' )
    2053                  DO  k = nzb_s_inner(j,i), nzt
    2054                    sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)                    &
    2055                       + ( flux_t(k) + diss_t(k) )                             &
    2056                       *   weight_substep(intermediate_timestep_count)
    2057                  ENDDO
    2058                CASE ( 'sa' )
    2059                  DO  k = nzb_s_inner(j,i), nzt
    2060                    sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)                    &
    2061                       + ( flux_t(k) + diss_t(k) )                             &
    2062                       *   weight_substep(intermediate_timestep_count)
    2063                  ENDDO
    2064                CASE ( 'q' )
    2065                  DO  k = nzb_s_inner(j,i), nzt
    2066                    sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)                      &
    2067                       + ( flux_t(k) + diss_t(k) )                             &
    2068                       *   weight_substep(intermediate_timestep_count)
    2069                  ENDDO
    2070 
    2071             END SELECT
    2072          ENDDO
    2073       ENDDO
    2074 
    20752181
    20762182    END SUBROUTINE advec_s_ws
    … …  
    20912197       IMPLICIT NONE
    20922198
    2093        INTEGER ::  i, j, k, tn = 0
    2094        REAL    ::  gu, gv, flux_d, diss_d, v_comp, w_comp
     2199       INTEGER ::  i, j, k, tn = 0,  k_ppp, k_pp, k_mm
     2200       REAL    ::  gu, gv, flux_d, diss_d, v_comp, w_comp, div
    20952201       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u
    2096        REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u,             &
     2202       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u,           &
    20972203                                             swap_diss_x_local_u
    2098        REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n,  &
     2204       REAL, DIMENSION(nzb:nzt) :: flux_t, diss_t, flux_r, diss_r, flux_n,  &
    20992205                                     diss_n, u_comp
    21002206 
    … …  
    21062212       i = nxlu
    21072213       DO  j = nys, nyn
    2108           IF( boundary_flags(j,i) == 5 )  THEN
    2109          
    2110              DO  k = nzb_u_inner(j,i)+1, nzt
    2111                 u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
    2112                 swap_flux_x_local_u(k,j) = u_comp(k) *                         &
    2113                                            ( u(k,j,i) + u(k,j,i-1) ) * 0.25
    2114                 swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),   &
    2115                                                      u(k,j,i), u(k,j,i-1),     &
    2116                                                      u(k,j,i-1), u_comp(k),    &
    2117                                                      0.25, ddx )
    2118              ENDDO
    2119              
    2120           ELSE
    2121          
    2122              DO  k = nzb_u_inner(j,i)+1, nzt
    2123                 u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
    2124                 swap_flux_x_local_u(k,j) = u_comp(k) * (                       &
    2125                                            37.0 * ( u(k,j,i) + u(k,j,i-1)   )  &
    2126                                          -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )  &
    2127                                          +        (u(k,j,i+2)+u(k,j,i-3) )  )  &
    2128                                          * adv_mom_5
    2129                 swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                &
    2130                                            10.0 * ( u(k,j,i) - u(k,j,i-1)   )  &
    2131                                          -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )  &
    2132                                          +        ( u(k,j,i+2) - u(k,j,i-3) ) )&
    2133                                          * adv_mom_5
    2134              ENDDO
    2135              
    2136           ENDIF
    2137          
     2214          DO  k = nzb+1, nzb_max
     2215
     2216             u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
     2217             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
     2218                          ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     2219                       +     7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     2220                       +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1 &
     2221                          ) *                                                  &
     2222                                     ( u(k,j,i)   + u(k,j,i-1) )               &
     2223                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     2224                       +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     2225                          ) *                                                  &
     2226                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
     2227                   +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     2228                          ) *                                                  &
     2229                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
     2230                                                   )
     2231
     2232              swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * (                &
     2233                          ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     2234                       +     3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     2235                       +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1 &
     2236                          ) *                                                  &
     2237                                     ( u(k,j,i)   - u(k,j,i-1) )               &
     2238                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     2239                       +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
     2240                          ) *                                                  &
     2241                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
     2242                   +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
     2243                          ) *                                                  &
     2244                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
     2245                                                             )
     2246
     2247          ENDDO
     2248
     2249          DO  k = nzb_max+1, nzt
     2250
     2251             u_comp(k)         = u(k,j,i) + u(k,j,i-1) - gu
     2252             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
     2253                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )                &
     2254                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )                &
     2255                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
     2256             swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                   &
     2257                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )                &
     2258                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )                &
     2259                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
     2260
     2261          ENDDO
    21382262       ENDDO
    21392263
    21402264       DO i = nxlu, nxr
    21412265!       
    2142 !--      The following loop computes the fluxes for the south boundary points
    2143          j = nys
    2144          IF ( boundary_flags(j,i) == 8 )  THEN
    2145 !         
    2146 !--         Compute southside fluxes for the south boundary of PE domain
    2147             DO  k = nzb_u_inner(j,i)+1, nzt
    2148                v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
    2149                swap_flux_y_local_u(k) = v_comp *                              &
    2150                                         ( u(k,j,i) + u(k,j-1,i) ) * 0.25
    2151                swap_diss_y_local_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i),     &
    2152                                                   u(k,j,i), u(k,j-1,i),       &
    2153                                                   u(k,j-1,i), v_comp,         &
    2154                                                   0.25, ddy )
    2155             ENDDO
    2156            
    2157          ELSE
    2158          
    2159             DO  k = nzb_u_inner(j,i)+1, nzt
    2160                v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
    2161                swap_flux_y_local_u(k) = v_comp * (                             &
    2162                                         37.0 * ( u(k,j,i) + u(k,j-1,i)   )     &
    2163                                       -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )     &
    2164                                       +        ( u(k,j+2,i) + u(k,j-3,i) ) )   &
    2165                                       * adv_mom_5
    2166                swap_diss_y_local_u(k) = - ABS( v_comp ) * (                    &
    2167                                         10.0 * ( u(k,j,i) - u(k,j-1,i)   )     &
    2168                                       -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )     &
    2169                                       +        ( u(k,j+2,i) - u(k,j-3,i) ) )   &
    2170                                       * adv_mom_5
    2171             ENDDO
    2172            
    2173          ENDIF
    2174 !         
    2175 !--      Computation of interior fluxes and tendency terms
    2176          DO  j = nys, nyn
    2177             IF ( boundary_flags(j,i) /= 0 )  THEN
    2178 !           
    2179 !--            Degrade the order for Dirichlet bc. at the outflow boundary
    2180                SELECT CASE ( boundary_flags(j,i) )
    2181                
    2182                   CASE ( 1 )
    2183                      DO  k = nzb_u_inner(j,i)+1, nzt
    2184                         u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2185                         flux_r(k) = ( u_comp(k) - gu ) * (                     &
    2186                                     7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
    2187                                   -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
    2188                                   * adv_mom_3
    2189                         diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
    2190                                     3.0 * ( u(k,j,i+1) - u(k,j,i) )            &
    2191                                   -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
    2192                                   * adv_mom_3
    2193                      ENDDO
    2194                      
    2195                   CASE ( 2 )
    2196                      DO  k = nzb_u_inner(j,i)+1, nzt
    2197                         u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2198                         flux_r(k) = ( u_comp(k) - gu ) *                      &
    2199                                     ( u(k,j,i+1) + u(k,j,i) ) * 0.25
    2200                         diss_r(k) = diss_2nd( u(k,j,i+1), u(k,j,i+1),         &
    2201                                               u(k,j,i), u(k,j,i-1),           &
    2202                                               u(k,j,i-2), u_comp(k) ,0.25 ,ddx)
    2203                      ENDDO
    2204                      
    2205                   CASE ( 3 )
    2206                      DO  k = nzb_u_inner(j,i)+1, nzt
    2207                         v_comp   = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2208                         flux_n(k) = v_comp * (                                 &
    2209                                     7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
    2210                                   -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
    2211                                   * adv_mom_3
    2212                         diss_n(k) = - ABS( v_comp ) * (                        &
    2213                                     3.0 * ( u(k,j+1,i) - u(k,j,i)   )          &
    2214                                   -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
    2215                                   * adv_mom_3
    2216                      ENDDO
    2217                      
    2218                   CASE ( 4 )
    2219                      DO  k = nzb_u_inner(j,i)+1, nzt
    2220                         v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2221                         flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
    2222                         diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i),         &
    2223                                               u(k,j,i), u(k,j-1,i),           &
    2224                                               u(k,j-2,i), v_comp, 0.25, ddy )
    2225                      ENDDO
    2226                      
    2227                   CASE ( 5 )
    2228                      DO  k = nzb_u_inner(j,i)+1, nzt       
    2229                         u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2230                         flux_r(k) = ( u_comp(k) - gu ) * (                     &
    2231                                     7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
    2232                                   -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
    2233                                   * adv_mom_3
    2234                         diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
    2235                                     3.0 * ( u(k,j,i+1) - u(k,j,i)   )          &
    2236                                   -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
    2237                                   * adv_mom_3
    2238                      ENDDO
    2239                      
    2240                   CASE ( 7 )
    2241                      DO  k = nzb_u_inner(j,i)+1, nzt                           
    2242                         v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2243                         flux_n(k) = v_comp * (                                 &
    2244                                     7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
    2245                                   -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
    2246                                   * adv_mom_3
    2247                         diss_n(k) = - ABS( v_comp ) * (                        &
    2248                                     3.0 * ( u(k,j+1,i) - u(k,j,i)   )          &
    2249                                   -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
    2250                                   * adv_mom_3
    2251                      ENDDO
    2252                      
    2253                   CASE ( 8 )
    2254                      DO  k = nzb_u_inner(j,i)+1, nzt
    2255                         v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2256                         flux_n(k) = v_comp * (                                 &
    2257                                     7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
    2258                                   -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
    2259                                   * adv_mom_3
    2260                         diss_n(k) = - ABS( v_comp ) * (                        &
    2261                                     3.0 * ( u(k,j+1,i) - u(k,j,i)   )          &
    2262                                   -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
    2263                                   * adv_mom_3
    2264                      ENDDO 
    2265                      
    2266                   CASE DEFAULT
    2267                  
    2268                END SELECT
    2269 !               
    2270 !--            Compute the crosswise 5th order fluxes at the outflow
    2271                IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
    2272                     boundary_flags(j,i) == 5 )  THEN
    2273              
    2274                   DO  k = nzb_u_inner(j,i)+1, nzt
    2275                      v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2276                      flux_n(k) = v_comp * (                                    &
    2277                                  37.0 * ( u(k,j+1,i) + u(k,j,i)   )            &
    2278                                -  8.0 * ( u(k,j+2,i) +u(k,j-1,i)  )            &
    2279                                +        ( u(k,j+3,i) + u(k,j-2,i) ) )          &
    2280                                * adv_mom_5
    2281                      diss_n(k) = - ABS( v_comp ) * (                           &
    2282                                  10.0 * ( u(k,j+1,i) - u(k,j,i)   )            &
    2283                                -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )            &
    2284                                +        ( u(k,j+3,i) - u(k,j-2,i) ) )          &
    2285                                * adv_mom_5
    2286                   ENDDO
    2287                  
    2288                ELSE
    2289                
    2290                   DO  k = nzb_u_inner(j,i)+1, nzt
    2291                      u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2292                      flux_r(k) = ( u_comp(k) - gu ) * (                        &
    2293                                  37.0 * ( u(k,j,i+1) + u(k,j,i)   )            &
    2294                                -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )            &
    2295                                +        ( u(k,j,i+3) + u(k,j,i-2) ) )          &
    2296                                * adv_mom_5
    2297                      diss_r(k) = - ABS( u_comp(k) - gu ) * (                   &
    2298                                  10.0 * ( u(k,j,i+1) - u(k,j,i)   )            &
    2299                                -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )            &
    2300                                +        ( u(k,j,i+3) - u(k,j,i-2) ) )          &
    2301                                * adv_mom_5
    2302                   ENDDO
    2303                  
    2304                ENDIF
    2305                
    2306             ELSE
    2307            
    2308                DO  k = nzb_u_inner(j,i)+1, nzt
    2309                   u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2310                   flux_r(k) = ( u_comp(k) - gu ) * (                           &
    2311                               37.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
    2312                             -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )               &
    2313                             +        ( u(k,j,i+3) + u(k,j,i-2) ) )             &
    2314                             * adv_mom_5
    2315                   diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
    2316                               10.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
    2317                             -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )               &
    2318                             +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    2319 
    2320                   v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2321                   flux_n(k) = v_comp * (                                       &
    2322                               37.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
    2323                             -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )               &
    2324                             +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    2325                   diss_n(k) = - ABS( v_comp ) * (                              &
    2326                               10.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
    2327                             -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )               &
    2328                             +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    2329                                
    2330                ENDDO
    2331                
    2332             ENDIF
    2333            
    2334             DO  k = nzb_u_inner(j,i)+1, nzt
    2335            
    2336                tend(k,j,i) = tend(k,j,i) - (                                   &
    2337               ( flux_r(k) + diss_r(k)                                          &
    2338             -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j)   ) * ddx  &
    2339             + ( flux_n(k) + diss_n(k)                                          &
    2340             -   swap_flux_y_local_u(k) - swap_diss_y_local_u(k)       ) * ddy )
    2341                  
    2342                swap_flux_x_local_u(k,j) = flux_r(k)   
    2343                swap_diss_x_local_u(k,j) = diss_r(k)
    2344                swap_flux_y_local_u(k)   = flux_n(k)     
    2345                swap_diss_y_local_u(k)   = diss_n(k)
    2346                      
    2347                sums_us2_ws_l(k,tn)  = sums_us2_ws_l(k,tn)                            &
    2348                  + ( flux_r(k)                                                 &
    2349                  * ( u_comp(k) - 2.0 * hom(k,1,1,0) )                          &
    2350                  / ( u_comp(k) - gu + 1.0E-20 )                                &
    2351                  +   diss_r(k)                                                 &
    2352                  *   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )                     &
    2353                  / ( ABS( u_comp(k) - gu) + 1.0E-20) )                         &
    2354                  *   weight_substep(intermediate_timestep_count)
    2355             ENDDO
    2356             sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
    2357          ENDDO
     2266!--       The following loop computes the fluxes for the south boundary points
     2267          j = nys
     2268          DO  k = nzb+1, nzb_max
     2269
     2270             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
     2271             swap_flux_y_local_u(k) = v_comp * (                              &
     2272                         ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     2273                      +     7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     2274                      +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
     2275                         ) *                                                  &
     2276                                     ( u(k,j,i)   + u(k,j-1,i) )              &
     2277                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     2278                      +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     2279                         ) *                                                  &
     2280                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
     2281                  +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     2282                         ) *                                                  &
     2283                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
     2284                                               )
     2285
     2286             swap_diss_y_local_u(k) = - ABS ( v_comp ) * (                    &
     2287                         ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     2288                      +     3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     2289                      +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
     2290                         ) *                                                  &
     2291                                     ( u(k,j,i)   - u(k,j-1,i) )              &
     2292                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     2293                      +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
     2294                         ) *                                                  &
     2295                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
     2296                  +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
     2297                         ) *                                                  &
     2298                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
     2299                                                         )
     2300
     2301          ENDDO
     2302
     2303          DO  k = nzb_max+1, nzt
     2304
     2305             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
     2306             swap_flux_y_local_u(k) = v_comp * (                              &
     2307                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
     2308                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
     2309                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
     2310             swap_diss_y_local_u(k) = - ABS(v_comp) * (                       &
     2311                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
     2312                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
     2313                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
     2314
     2315          ENDDO
     2316!
     2317!--       Computation of interior fluxes and tendency terms
     2318          DO  j = nys, nyn
     2319
     2320             flux_t(0) = 0.0
     2321             diss_t(0) = 0.0
     2322             flux_d    = 0.0
     2323             diss_d    = 0.0
     2324
     2325             DO  k = nzb+1, nzb_max
     2326
     2327                u_comp(k) = u(k,j,i+1) + u(k,j,i)
     2328                flux_r(k) = ( u_comp(k) - gu ) * (                           &
     2329                     ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     2330                  +     7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     2331                  +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1    &
     2332                     ) *                                                     &
     2333                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
     2334              -      (  8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     2335                  +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     2336                     ) *                                                     &
     2337                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
     2338              +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     2339                     ) *                                                     &
     2340                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
     2341                                                 )
     2342
     2343                diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
     2344                     ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     2345                  +     3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     2346                  +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1    &
     2347                     ) *                                                     &
     2348                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
     2349              -      (  5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     2350                  +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
     2351                     ) *                                                     &
     2352                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
     2353              +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
     2354                     ) *                                                     &
     2355                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
     2356                                                     )
     2357
     2358                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
     2359                flux_n(k) = v_comp * (                                       &
     2360                     ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     2361                  +     7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     2362                  +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1    &
     2363                     ) *                                                     &
     2364                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
     2365              -      (  8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     2366                  +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     2367                     ) *                                                     &
     2368                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
     2369              +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     2370                     ) *                                                     &
     2371                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
     2372                                                 )
     2373
     2374                diss_n(k) = - ABS ( v_comp ) * (                             &
     2375                     ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     2376                  +     3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     2377                  +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1    &
     2378                     ) *                                                     &
     2379                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
     2380              -      (  5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     2381                  +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
     2382                     ) *                                                     &
     2383                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
     2384              +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
     2385                     ) *                                                     &
     2386                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
     2387                                                      )
     2388!
     2389!--             k index has to be modified near bottom and top, else array
     2390!--             subscripts will be exceeded.
     2391                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
     2392                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1)  )
     2393                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
     2394
     2395                w_comp    = w(k,j,i) + w(k,j,i-1)
     2396                flux_t(k) = w_comp  * (                                      &
     2397                     ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2398                  +     7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2399                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     2400                     ) *                                                     &
     2401                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
     2402              -      (  8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2403                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2404                     ) *                                                     &
     2405                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
     2406              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2407                     ) *                                                     &
     2408                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     2409                                      )
     2410
     2411                diss_t(k) = - ABS( w_comp ) * (                              &
     2412                     ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2413                  +     3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2414                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     2415                     ) *                                                     &
     2416                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
     2417              -      (  5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2418                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2419                     ) *                                                     &
     2420                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
     2421              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2422                     ) *                                                     &
     2423                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     2424                                              )
     2425!
     2426!--             Calculate the divergence of the velocity field. A respective
     2427!--             correction is needed to overcome numerical instabilities caused
     2428!--             by a not sufficient reduction of divergences near topography.
     2429                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
     2430                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
     2431                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
     2432                                                                    * ddzw(k)  &
     2433                      ) * 0.5
     2434
     2435                tend(k,j,i) = tend(k,j,i) - (                                  &
     2436                 ( flux_r(k) + diss_r(k)                                       &
     2437               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
     2438               + ( flux_n(k) + diss_n(k)                                       &
     2439               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
     2440               + ( flux_t(k) + diss_t(k)                                       &
     2441               -   flux_d    - diss_d                                          &
     2442                                                                  ) * ddzw(k)  &
     2443                                           ) + div * u(k,j,i)
     2444
     2445                swap_flux_x_local_u(k,j) = flux_r(k)
     2446                swap_diss_x_local_u(k,j) = diss_r(k)
     2447                swap_flux_y_local_u(k)   = flux_n(k)
     2448                swap_diss_y_local_u(k)   = diss_n(k)
     2449                flux_d                   = flux_t(k)
     2450                diss_d                   = diss_t(k)
     2451!
     2452!--             Statistical Evaluation of u'u'. The factor has to be applied
     2453!--             for right evaluation when gallilei_trans = .T. .
     2454                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                     &
     2455                              + ( flux_r(k) *                                 &
     2456                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
     2457                              / ( u_comp(k) - gu + 1.0E-20      )             &
     2458                              +   diss_r(k) *                                 &
     2459                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
     2460                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
     2461                              *   weight_substep(intermediate_timestep_count)
     2462!
     2463!--             Statistical Evaluation of w'u'.
     2464                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                   &
     2465                              + ( flux_t(k) + diss_t(k) )                     &
     2466                              *   weight_substep(intermediate_timestep_count)
     2467             ENDDO
     2468
     2469             DO  k = nzb_max+1, nzt
     2470
     2471                u_comp(k) = u(k,j,i+1) + u(k,j,i)
     2472                flux_r(k) = ( u_comp(k) - gu ) * (                            &
     2473                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
     2474                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
     2475                       +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
     2476                diss_r(k) = - ABS( u_comp(k) - gu ) * (                       &
     2477                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
     2478                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
     2479                       +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     2480
     2481                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
     2482                flux_n(k) = v_comp * (                                        &
     2483                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
     2484                       -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
     2485                       +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     2486                diss_n(k) = - ABS( v_comp ) * (                               &
     2487                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
     2488                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
     2489                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     2490!
     2491!--             k index has to be modified near bottom and top, else array
     2492!--             subscripts will be exceeded.
     2493                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
     2494                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1)  )
     2495                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
     2496
     2497                w_comp    = w(k,j,i) + w(k,j,i-1)
     2498                flux_t(k) = w_comp  * (                                      &
     2499                     ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2500                  +     7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2501                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     2502                     ) *                                                     &
     2503                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
     2504              -      (  8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2505                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2506                     ) *                                                     &
     2507                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
     2508              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2509                     ) *                                                     &
     2510                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     2511                                      )
     2512
     2513                 diss_t(k) = - ABS( w_comp ) * (                             &
     2514                     ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2515                  +     3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2516                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
     2517                     ) *                                                     &
     2518                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
     2519              -      (  5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2520                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
     2521                     ) *                                                     &
     2522                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
     2523              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
     2524                     ) *                                                     &
     2525                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     2526                                               )
     2527!
     2528!--             Calculate the divergence of the velocity field. A respective
     2529!--             correction is needed to overcome numerical instabilities caused
     2530!--             by a not sufficient reduction of divergences near topography.
     2531                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
     2532                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
     2533                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )       &
     2534                                                                    * ddzw(k) &
     2535                      ) * 0.5
     2536
     2537                 tend(k,j,i) = tend(k,j,i) - (                                 &
     2538                 ( flux_r(k) + diss_r(k)                                       &
     2539               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
     2540               + ( flux_n(k) + diss_n(k)                                       &
     2541               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
     2542               + ( flux_t(k) + diss_t(k)                                       &
     2543               -   flux_d    - diss_d                                          &
     2544                                                                  ) * ddzw(k)  &
     2545                                           ) + div * u(k,j,i)
     2546
     2547                 swap_flux_x_local_u(k,j) = flux_r(k)
     2548                 swap_diss_x_local_u(k,j) = diss_r(k)
     2549                 swap_flux_y_local_u(k)   = flux_n(k)
     2550                 swap_diss_y_local_u(k)   = diss_n(k)
     2551                 flux_d                   = flux_t(k)
     2552                 diss_d                   = diss_t(k)
     2553!
     2554!--              Statistical Evaluation of u'u'. The factor has to be applied
     2555!--              for right evaluation when gallilei_trans = .T. .
     2556                 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                    &
     2557                              + ( flux_r(k) *                                 &
     2558                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
     2559                              / ( u_comp(k) - gu + 1.0E-20      )             &
     2560                              +   diss_r(k) *                                 &
     2561                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
     2562                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
     2563                              *   weight_substep(intermediate_timestep_count)
     2564!
     2565!--              Statistical Evaluation of w'u'.
     2566                 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                  &
     2567                              + ( flux_t(k) + diss_t(k) )                     &
     2568                              *   weight_substep(intermediate_timestep_count)
    23582569       ENDDO
    2359 
    2360 !
    2361 !--   Vertical advection, degradation of order near surface and top.
    2362 !--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    2363 !--   statistical evaluation the top flux at the surface should be 0
    2364        DO i = nxlu, nxr
    2365           DO  j = nys, nyn
    2366              k = nzb_u_inner(j,i)+1
    2367 !             
    2368 !--         The fluxes flux_d and diss_d at the surface are 0. Due to static
    2369 !--         reasons the top flux at the surface should be 0.
    2370             flux_t(nzb_u_inner(j,i)) = 0.0
    2371             diss_t(nzb_u_inner(j,i)) = 0.0
    2372             flux_d = flux_t(k-1)
    2373             diss_d = diss_t(k-1)             
    2374 !
    2375 !--         2nd order scheme (bottom)
    2376              w_comp    = w(k,j,i) + w(k,j,i-1)
    2377              flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
    2378              diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),           &
    2379                                    0.0, 0.0, w_comp, 0.25, ddzw(k) )
    2380                                    
    2381              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2382                                        -   flux_d    -  diss_d       ) * ddzw(k)
    2383 !
    2384 !--         WS3 as an intermediate step (bottom)
    2385             k         = nzb_u_inner(j,i)+2
    2386             flux_d    = flux_t(k-1)
    2387             diss_d    = diss_t(k-1)
    2388             w_comp    = w(k,j,i) + w(k,j,i-1)
    2389             flux_t(k) = w_comp * (                                            &
    2390                         7.0 * ( u(k+1,j,i) + u(k,j,i)   )                     &
    2391                       -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    2392             diss_t(k) = - ABS( w_comp ) * (                                   &
    2393                         3.0 * ( u(k+1,j,i) - u(k,j,i)   )                     &
    2394                       -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    2395 
    2396             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
    2397                                       -   flux_d    -  diss_d      ) * ddzw(k)
    2398 !
    2399 !WS5
    2400              DO  k = nzb_u_inner(j,i)+3, nzt-2
    2401 
    2402                 flux_d    = flux_t(k-1)
    2403                 diss_d    = diss_t(k-1)
    2404                 w_comp    = w(k,j,i) + w(k,j,i-1)
    2405                 flux_t(k) = w_comp * (                                        &
    2406                             37.0 * ( u(k+1,j,i) + u(k,j,i)   )                &
    2407                           -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                &
    2408                           +        ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
    2409                 diss_t(k) = - ABS( w_comp ) * (                               &
    2410                             10.0 * ( u(k+1,j,i) - u(k,j,i)   )                &
    2411                           -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                &
    2412                           +        ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
    2413 
    2414                 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
    2415                                           -   flux_d    -  diss_d   ) * ddzw(k)
    2416 
    2417              ENDDO
    2418 !
    2419 !--          WS3 as an intermediate step (top)
    2420              k         = nzt-1
    2421              flux_d    = flux_t(k-1)
    2422              diss_d    = diss_t(k-1)
    2423              w_comp    = w(k,j,i) + w(k,j,i-1)
    2424              flux_t(k) = w_comp * (                                           &
    2425                          7.0 * ( u(k+1,j,i) + u(k,j,i) )                      &
    2426                        -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
    2427              diss_t(k) = - ABS( w_comp ) * (                                  &
    2428                          3.0 * ( u(k+1,j,i) - u(k,j,i) )                      &
    2429                        -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
    2430                            
    2431              tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
    2432                                        -   flux_d    -  diss_d      ) * ddzw(k)
    2433 !
    2434 !--         2nd order scheme (top)