Changeset 3551


Ignore:
Timestamp:
Nov 21, 2018 5:05:46 PM (3 years ago)
Author:
suehring
Message:

optimization of advec_ws

File:
1 edited

Legend:

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

    r3547 r3551  
    2525! -----------------
    2626! $Id$
     27! - Computation of vertical fluxes separated from computation of horizontal
     28!   fluxes. Loops are splitted in order to avoid indirect indexing and allow
     29!   for better vectorization.
     30! - Accelerate code by remove type-conversions of ibits
     31! - replace pointer definition in scalar routine by simple explicit definition
     32!   of the passed array
     33!
     34! 3547 2018-11-21 13:21:24Z suehring
    2735! variables documented
    2836!
     
    193201!
    194202! 888 2012-04-20 15:03:46Z suehring
    195 ! Number of IBITS() calls with identical arguments is reduced.
     203! Number of REAL( IBITS() calls with identical arguments is reduced.
    196204!
    197205! 862 2012-03-26 14:21:38Z suehring
     
    12001208       
    12011209       INTEGER(iwp) ::  i     !< grid index along x-direction
    1202        INTEGER(iwp) ::  ibit0 !< flag indicating 1st-order scheme along x-direction
    1203        INTEGER(iwp) ::  ibit1 !< flag indicating 3rd-order scheme along x-direction
    1204        INTEGER(iwp) ::  ibit2 !< flag indicating 5th-order scheme along x-direction
    1205        INTEGER(iwp) ::  ibit3 !< flag indicating 1st-order scheme along y-direction
    1206        INTEGER(iwp) ::  ibit4 !< flag indicating 3rd-order scheme along y-direction
    1207        INTEGER(iwp) ::  ibit5 !< flag indicating 5th-order scheme along y-direction
    1208        INTEGER(iwp) ::  ibit6 !< flag indicating 1st-order scheme along z-direction
    1209        INTEGER(iwp) ::  ibit7 !< flag indicating 3rd-order scheme along z-direction
    1210        INTEGER(iwp) ::  ibit8 !< flag indicating 5th-order scheme along z-direction
    12111210       INTEGER(iwp) ::  i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread
    12121211       INTEGER(iwp) ::  j     !< grid index along y-direction
     
    12171216       INTEGER(iwp) ::  tn    !< number of OpenMP thread
    12181217       
     1218       REAL(wp)     ::  ibit0  !< flag indicating 1st-order scheme along x-direction
     1219       REAL(wp)     ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
     1220       REAL(wp)     ::  ibit2  !< flag indicating 5th-order scheme along x-direction
     1221       REAL(wp)     ::  ibit3  !< flag indicating 1st-order scheme along y-direction
     1222       REAL(wp)     ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
     1223       REAL(wp)     ::  ibit5  !< flag indicating 5th-order scheme along y-direction
     1224       REAL(wp)     ::  ibit6  !< flag indicating 1st-order scheme along z-direction
     1225       REAL(wp)     ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
     1226       REAL(wp)     ::  ibit8  !< flag indicating 5th-order scheme along z-direction
    12191227       REAL(wp)     ::  diss_d !< artificial dissipation term at grid box bottom
    12201228       REAL(wp)     ::  div    !< diverence on scalar grid
     
    12221230       REAL(wp)     ::  u_comp !< advection velocity along x-direction
    12231231       REAL(wp)     ::  v_comp !< advection velocity along y-direction
    1224        
    1225 #if defined( __nopointer )
    1226        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< advected scalar
    1227 #else
    1228        REAL(wp), DIMENSION(:,:,:), POINTER    ::  sk     !< advected scalar
    1229 #endif
     1232!       
     1233!--    sk is an array from parameter list. It should not be a pointer, because
     1234!--    in that case the compiler can not assume a stride 1 and cannot perform
     1235!--    a strided one vector load. Adding the CONTIGUOUS keyword makes things
     1236!--    even worse, because the compiler cannot assume strided one in the
     1237!--    caller side.
     1238       REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
     1239
    12301240       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
    12311241       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
     
    12491259          DO  k = nzb+1, nzb_max
    12501260
    1251              ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
    1252              ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
    1253              ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
     1261             ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
     1262             ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
     1263             ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
    12541264
    12551265             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     
    13101320          DO  k = nzb+1, nzb_max
    13111321
    1312              ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
    1313              ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
    1314              ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
     1322             ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
     1323             ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
     1324             ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
    13151325
    13161326             u_comp                     = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    13651375           
    13661376       ENDIF
    1367 
    1368        flux_t(0) = 0.0_wp
    1369        diss_t(0) = 0.0_wp
    13701377!       
    13711378!--    Now compute the fluxes and tendency terms for the horizontal and
     
    13771384!--       flux at the end.
    13781385
    1379           ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
    1380           ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
    1381           ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
     1386          ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
     1387          ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
     1388          ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
    13821389
    13831390          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     
    14121419                                       )
    14131420
    1414           ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
    1415           ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
    1416           ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
     1421          ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
     1422          ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
     1423          ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
    14171424
    14181425          v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
     
    14461453                             ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
    14471454                                       )
    1448 !
    1449 !--       k index has to be modified near bottom and top, else array
    1450 !--       subscripts will be exceeded.
    1451           ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
    1452           ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
    1453           ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
    1454 
    1455           k_ppp = k + 3 * ibit8
    1456           k_pp  = k + 2 * ( 1 - ibit6  )
    1457           k_mm  = k - 2 * ibit8
    1458 
    1459 
    1460           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1461                      ( 37.0_wp * ibit8 * adv_sca_5                            &
    1462                   +     7.0_wp * ibit7 * adv_sca_3                            &
    1463                   +              ibit6 * adv_sca_1                            &
    1464                      ) *                                                      &
    1465                              ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
    1466               -      (  8.0_wp * ibit8 * adv_sca_5                            &
    1467                   +              ibit7 * adv_sca_3                            &
    1468                      ) *                                                      &
    1469                              ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
    1470               +      (           ibit8 * adv_sca_5                            &
    1471                      ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
    1472                                  )
    1473 
    1474           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1475                      ( 10.0_wp * ibit8 * adv_sca_5                            &
    1476                   +     3.0_wp * ibit7 * adv_sca_3                            &
    1477                   +              ibit6 * adv_sca_1                            &
    1478                      ) *                                                      &
    1479                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1480               -      (  5.0_wp * ibit8 * adv_sca_5                            &
    1481                   +              ibit7 * adv_sca_3                            &
    1482                      ) *                                                      &
    1483                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
    1484               +      (           ibit8 * adv_sca_5                            &
    1485                      ) *                                                      &
    1486                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    1487                                          )
    1488        ENDDO
    1489 
    1490        DO  k = nzb+1, nzb_max
    1491 
    1492           flux_d    = flux_t(k-1)
    1493           diss_d    = diss_t(k-1)
    1494 !
    1495 !--       Calculate the divergence of the velocity field. A respective
    1496 !--       correction is needed to overcome numerical instabilities caused
    1497 !--       by a not sufficient reduction of divergences near topography.
    1498           div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
    1499                           - u(k,j,i)   * ( IBITS(advc_flags_1(k,j,i-1),0,1)    &
    1500                                          + IBITS(advc_flags_1(k,j,i-1),1,1)    &
    1501                                          + IBITS(advc_flags_1(k,j,i-1),2,1)    &
    1502                                          )                                     &
    1503                           ) * ddx                                              &
    1504                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    1505                           - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
    1506                                          + IBITS(advc_flags_1(k,j-1,i),4,1)    &
    1507                                          + IBITS(advc_flags_1(k,j-1,i),5,1)    &
    1508                                          )                                     &
    1509                           ) * ddy                                              &
    1510                         + ( w(k,j,i) * rho_air_zw(k) *                         &
    1511                                          ( ibit6 + ibit7 + ibit8 )             &
    1512                           - w(k-1,j,i) * rho_air_zw(k-1) *                     &
    1513                                          ( IBITS(advc_flags_1(k-1,j,i),6,1)    &
    1514                                          + IBITS(advc_flags_1(k-1,j,i),7,1)    &
    1515                                          + IBITS(advc_flags_1(k-1,j,i),8,1)    &
    1516                                          )                                     &     
    1517                           ) * drho_air(k) * ddzw(k)
    1518 
    1519 
    1520           tend(k,j,i) = tend(k,j,i) - (                                       &
    1521                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
    1522                           swap_diss_x_local(k,j,tn)            ) * ddx        &
    1523                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    1524                           swap_diss_y_local(k,tn)              ) * ddy        &
    1525                       + ( ( flux_t(k) + diss_t(k) ) -                         &
    1526                           ( flux_d    + diss_d    )                           &
    1527                                                     ) * drho_air(k) * ddzw(k) &
    1528                                       ) + sk(k,j,i) * div
    1529 
    1530           swap_flux_y_local(k,tn)   = flux_n(k)
    1531           swap_diss_y_local(k,tn)   = diss_n(k)
    1532           swap_flux_x_local(k,j,tn) = flux_r(k)
    1533           swap_diss_x_local(k,j,tn) = diss_r(k)
    1534 
    15351455       ENDDO
    15361456!
     
    15591479                    -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
    15601480                    +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     1481
     1482       ENDDO
     1483!
     1484!--    Now, compute vertical fluxes. Split loop into a part treating the
     1485!--    lowest 2 grid points with indirect indexing, a main loop without
     1486!--    indirect indexing, and a loop for the uppermost 2 grip points with
     1487!--    indirect indexing. This allows better vectorization for the main loop.
     1488!--    First, compute the flux at model surface, which need has to be
     1489!--    calculated explicetely for the tendency at
     1490!--    the first w-level. For topography wall this is done implicitely by
     1491!--    advc_flags_1.
     1492       flux_t(nzb) = 0.0_wp
     1493       diss_t(nzb) = 0.0_wp
     1494       DO  k = nzb+1, nzb+2
     1495          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
     1496          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
     1497          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
    15611498!
    15621499!--       k index has to be modified near bottom and top, else array
    15631500!--       subscripts will be exceeded.
    1564           ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
    1565           ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
    1566           ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
    1567 
    15681501          k_ppp = k + 3 * ibit8
    15691502          k_pp  = k + 2 * ( 1 - ibit6  )
    15701503          k_mm  = k - 2 * ibit8
    15711504
     1505
    15721506          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1573                     ( 37.0_wp * ibit8 * adv_sca_5                             &
    1574                  +     7.0_wp * ibit7 * adv_sca_3                             &
    1575                  +              ibit6 * adv_sca_1                             &
    1576                     ) *                                                       &
    1577                              ( sk(k+1,j,i)  + sk(k,j,i)   )                   &
    1578               -     (  8.0_wp * ibit8 * adv_sca_5                             &
     1507                     ( 37.0_wp * ibit8 * adv_sca_5                            &
     1508                  +     7.0_wp * ibit7 * adv_sca_3                            &
     1509                  +              ibit6 * adv_sca_1                            &
     1510                     ) *                                                      &
     1511                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
     1512              -      (  8.0_wp * ibit8 * adv_sca_5                            &
    15791513                  +              ibit7 * adv_sca_3                            &
    1580                     ) *                                                       &
    1581                              ( sk(k_pp,j,i) + sk(k-1,j,i) )                   &
    1582               +     (           ibit8 * adv_sca_5                             &
    1583                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                   &
     1514                     ) *                                                      &
     1515                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
     1516              +      (           ibit8 * adv_sca_5                            &
     1517                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
    15841518                                 )
    15851519
    15861520          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1587                     ( 10.0_wp * ibit8 * adv_sca_5                             &
    1588                  +     3.0_wp * ibit7 * adv_sca_3                             &
    1589                  +              ibit6 * adv_sca_1                             &
    1590                     ) *                                                       &
     1521                     ( 10.0_wp * ibit8 * adv_sca_5                            &
     1522                  +     3.0_wp * ibit7 * adv_sca_3                            &
     1523                  +              ibit6 * adv_sca_1                            &
     1524                     ) *                                                      &
    15911525                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1592               -     (  5.0_wp * ibit8 * adv_sca_5                             &
    1593                  +              ibit7 * adv_sca_3                             &
    1594                     ) *                                                       &
     1526              -      (  5.0_wp * ibit8 * adv_sca_5                            &
     1527                  +              ibit7 * adv_sca_3                            &
     1528                     ) *                                                      &
    15951529                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
    1596               +     (           ibit8 * adv_sca_5                             &
    1597                     ) *                                                       &
     1530              +      (           ibit8 * adv_sca_5                            &
     1531                     ) *                                                      &
    15981532                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    15991533                                         )
    1600 
    1601        ENDDO
    1602 
    1603        DO  k = nzb_max+1, nzt
     1534       ENDDO
     1535       
     1536       DO  k = nzb+3, nzt-2
     1537          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
     1538          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
     1539          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     1540
     1541          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
     1542                     ( 37.0_wp * ibit8 * adv_sca_5                            &
     1543                  +     7.0_wp * ibit7 * adv_sca_3                            &
     1544                  +              ibit6 * adv_sca_1                            &
     1545                     ) *                                                      &
     1546                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
     1547              -      (  8.0_wp * ibit8 * adv_sca_5                            &
     1548                  +              ibit7 * adv_sca_3                            &
     1549                     ) *                                                      &
     1550                             ( sk(k+2,j,i) + sk(k-1,j,i)  )                   &
     1551              +      (           ibit8 * adv_sca_5                            &
     1552                     ) *     ( sk(k+3,j,i)+ sk(k-2,j,i) )                     &
     1553                                                 )
     1554
     1555          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
     1556                     ( 10.0_wp * ibit8 * adv_sca_5                            &
     1557                  +     3.0_wp * ibit7 * adv_sca_3                            &
     1558                  +              ibit6 * adv_sca_1                            &
     1559                     ) *                                                      &
     1560                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
     1561              -      (  5.0_wp * ibit8 * adv_sca_5                            &
     1562                  +              ibit7 * adv_sca_3                            &
     1563                     ) *                                                      &
     1564                             ( sk(k+2,j,i)  - sk(k-1,j,i)  )                  &
     1565              +      (           ibit8 * adv_sca_5                            &
     1566                     ) *                                                      &
     1567                             ( sk(k+3,j,i) - sk(k-2,j,i) )                    &
     1568                                                         )
     1569       ENDDO       
     1570       
     1571       DO  k = nzt-1, nzt
     1572          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
     1573          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
     1574          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     1575!
     1576!--       k index has to be modified near bottom and top, else array
     1577!--       subscripts will be exceeded.
     1578          k_ppp = k + 3 * ibit8
     1579          k_pp  = k + 2 * ( 1 - ibit6  )
     1580          k_mm  = k - 2 * ibit8
     1581
     1582
     1583          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
     1584                     ( 37.0_wp * ibit8 * adv_sca_5                            &
     1585                  +     7.0_wp * ibit7 * adv_sca_3                            &
     1586                  +              ibit6 * adv_sca_1                            &
     1587                     ) *                                                      &
     1588                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
     1589              -      (  8.0_wp * ibit8 * adv_sca_5                            &
     1590                  +              ibit7 * adv_sca_3                            &
     1591                     ) *                                                      &
     1592                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
     1593              +      (           ibit8 * adv_sca_5                            &
     1594                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
     1595                                                 )
     1596
     1597          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
     1598                     ( 10.0_wp * ibit8 * adv_sca_5                            &
     1599                  +     3.0_wp * ibit7 * adv_sca_3                            &
     1600                  +              ibit6 * adv_sca_1                            &
     1601                     ) *                                                      &
     1602                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
     1603              -      (  5.0_wp * ibit8 * adv_sca_5                            &
     1604                  +              ibit7 * adv_sca_3                            &
     1605                     ) *                                                      &
     1606                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
     1607              +      (           ibit8 * adv_sca_5                            &
     1608                     ) *                                                      &
     1609                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
     1610                                                         )
     1611       ENDDO
     1612
     1613       DO  k = nzb+1, nzt
    16041614
    16051615          flux_d    = flux_t(k-1)
     
    16101620!--       correction is needed to overcome numerical instabilities introduced
    16111621!--       by a not sufficient reduction of divergences near topography.
    1612           div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
    1613                         + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
    1614                         + ( w(k,j,i)   * rho_air_zw(k) -                      &
    1615                             w(k-1,j,i) * rho_air_zw(k-1)                      &
     1622          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )            &
     1623                          - u(k,j,i)   * (                                    &
     1624                        REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )   &
     1625                      + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )   &
     1626                      + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )   &
     1627                                         )                                    &
     1628                          ) * ddx                                             &
     1629                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )            &
     1630                          - v(k,j,i)   * (                                    &
     1631                        REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )   &
     1632                      + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )   &
     1633                      + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )   &
     1634                                         )                                    &
     1635                          ) * ddy                                             &
     1636                        + ( w(k,j,i) * rho_air_zw(k) *                        &
     1637                                         ( ibit6 + ibit7 + ibit8 )            &
     1638                          - w(k-1,j,i) * rho_air_zw(k-1) *                    &
     1639                                         (                                    &
     1640                        REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp )   &
     1641                      + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp )   &
     1642                      + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp )   &
     1643                                         )                                    &     
    16161644                          ) * drho_air(k) * ddzw(k)
    16171645
     
    17771805
    17781806       INTEGER(iwp) ::  i      !< grid index along x-direction
    1779        INTEGER(iwp) ::  ibit9  !< flag indicating 1st-order scheme along x-direction
    1780        INTEGER(iwp) ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
    1781        INTEGER(iwp) ::  ibit11 !< flag indicating 5th-order scheme along x-direction
    1782        INTEGER(iwp) ::  ibit12 !< flag indicating 1st-order scheme along y-direction
    1783        INTEGER(iwp) ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
    1784        INTEGER(iwp) ::  ibit14 !< flag indicating 5th-order scheme along y-direction
    1785        INTEGER(iwp) ::  ibit15 !< flag indicating 1st-order scheme along z-direction
    1786        INTEGER(iwp) ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
    1787        INTEGER(iwp) ::  ibit17 !< flag indicating 5th-order scheme along z-direction
    17881807       INTEGER(iwp) ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
    17891808       INTEGER(iwp) ::  j      !< grid index along y-direction
     
    17941813       INTEGER(iwp) ::  tn     !< number of OpenMP thread
    17951814       
     1815       REAL(wp)    ::  ibit9    !< flag indicating 1st-order scheme along x-direction
     1816       REAL(wp)    ::  ibit10   !< flag indicating 3rd-order scheme along x-direction
     1817       REAL(wp)    ::  ibit11   !< flag indicating 5th-order scheme along x-direction
     1818       REAL(wp)    ::  ibit12   !< flag indicating 1st-order scheme along y-direction
     1819       REAL(wp)    ::  ibit13   !< flag indicating 3rd-order scheme along y-direction
     1820       REAL(wp)    ::  ibit14   !< flag indicating 5th-order scheme along y-direction
     1821       REAL(wp)    ::  ibit15   !< flag indicating 1st-order scheme along z-direction
     1822       REAL(wp)    ::  ibit16   !< flag indicating 3rd-order scheme along z-direction
     1823       REAL(wp)    ::  ibit17   !< flag indicating 5th-order scheme along z-direction
    17961824       REAL(wp)    ::  diss_d   !< artificial dissipation term at grid box bottom
    17971825       REAL(wp)    ::  div      !< diverence on u-grid
     
    18191847          DO  k = nzb+1, nzb_max
    18201848
    1821              ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
    1822              ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
    1823              ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
     1849             ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
     1850             ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )
     1851             ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )
    18241852
    18251853             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
     
    18771905          DO  k = nzb+1, nzb_max
    18781906
    1879              ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
    1880              ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
    1881              ibit9  = IBITS(advc_flags_1(k,j,i-1),9,1)
     1907             ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
     1908             ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
     1909             ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
    18821910
    18831911             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     
    19291957         
    19301958       ENDIF
    1931 
    1932        flux_t(0) = 0.0_wp
    1933        diss_t(0) = 0.0_wp
    1934        w_comp(0) = 0.0_wp
    19351959!
    19361960!--    Now compute the fluxes tendency terms for the horizontal and
     
    19381962       DO  k = nzb+1, nzb_max
    19391963
    1940           ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
    1941           ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
    1942           ibit9  = IBITS(advc_flags_1(k,j,i),9,1)
     1964          ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
     1965          ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
     1966          ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1),  KIND = wp )
    19431967
    19441968          u_comp(k) = u(k,j,i+1) + u(k,j,i)
     
    19731997                                                )
    19741998
    1975           ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
    1976           ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
    1977           ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
     1999          ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
     2000          ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
     2001          ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
    19782002
    19792003          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
     
    20072031                                    ( u(k,j+3,i) - u(k,j-2,i) )               &
    20082032                                            )
     2033       ENDDO
     2034
     2035       DO  k = nzb_max+1, nzt
     2036
     2037          u_comp(k) = u(k,j,i+1) + u(k,j,i)
     2038          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
     2039                         37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                &
     2040                       -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
     2041                       +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
     2042          diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
     2043                         10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
     2044                       -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
     2045                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     2046
     2047          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
     2048          flux_n(k) = v_comp(k) * (                                           &
     2049                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
     2050                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
     2051                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     2052          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
     2053                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
     2054                       -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
     2055                       +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     2056
     2057       ENDDO
     2058!
     2059!--    Now, compute vertical fluxes. Split loop into a part treating the
     2060!--    lowest 2 grid points with indirect indexing, a main loop without
     2061!--    indirect indexing, and a loop for the uppermost 2 grip points with
     2062!--    indirect indexing. This allows better vectorization for the main loop.
     2063!--    First, compute the flux at model surface, which need has to be
     2064!--    calculated explicetely for the tendency at
     2065!--    the first w-level. For topography wall this is done implicitely by
     2066!--    advc_flags_1.
     2067       flux_t(nzb) = 0.0_wp
     2068       diss_t(nzb) = 0.0_wp
     2069       w_comp(nzb) = 0.0_wp
     2070       DO  k = nzb+1, nzb+2
    20092071!
    20102072!--       k index has to be modified near bottom and top, else array
    20112073!--       subscripts will be exceeded.
    2012           ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
    2013           ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
    2014           ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
     2074          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
     2075          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
     2076          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
    20152077
    20162078          k_ppp = k + 3 * ibit17
     
    20492111                                                           )
    20502112       ENDDO
    2051 
    2052        DO  k = nzb+1, nzb_max
     2113       
     2114       DO  k = nzb+3, nzt-2
     2115
     2116          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
     2117          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
     2118          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
     2119
     2120          w_comp(k) = w(k,j,i) + w(k,j,i-1)
     2121          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
     2122                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     2123                  +     7.0_wp * ibit16 * adv_mom_3                           &
     2124                  +              ibit15 * adv_mom_1                           &
     2125                     ) *                                                      &
     2126                                ( u(k+1,j,i)  + u(k,j,i)     )                &
     2127              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     2128                  +              ibit16 * adv_mom_3                           &
     2129                     ) *                                                      &
     2130                                ( u(k+2,j,i) + u(k-1,j,i)   )                 &
     2131              +      (           ibit17 * adv_mom_5                           &
     2132                     ) *                                                      &
     2133                                ( u(k+3,j,i) + u(k-2,j,i) )                   &
     2134                                                  )
     2135
     2136          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
     2137                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     2138                  +     3.0_wp * ibit16 * adv_mom_3                           &
     2139                  +              ibit15 * adv_mom_1                           &
     2140                     ) *                                                      &
     2141                                ( u(k+1,j,i)   - u(k,j,i)    )                &
     2142              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     2143                  +              ibit16 * adv_mom_3                           &
     2144                     ) *                                                      &
     2145                                ( u(k+2,j,i)  - u(k-1,j,i)  )                 &
     2146              +      (           ibit17 * adv_mom_5                           &
     2147                     ) *                                                      &
     2148                                ( u(k+3,j,i) - u(k-2,j,i) )                   &
     2149                                                           )
     2150       ENDDO
     2151       
     2152       DO  k = nzt-1, nzt
     2153!
     2154!--       k index has to be modified near bottom and top, else array
     2155!--       subscripts will be exceeded.
     2156          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
     2157          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
     2158          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
     2159
     2160          k_ppp = k + 3 * ibit17
     2161          k_pp  = k + 2 * ( 1 - ibit15 )
     2162          k_mm  = k - 2 * ibit17
     2163
     2164          w_comp(k) = w(k,j,i) + w(k,j,i-1)
     2165          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
     2166                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     2167                  +     7.0_wp * ibit16 * adv_mom_3                           &
     2168                  +              ibit15 * adv_mom_1                           &
     2169                     ) *                                                      &
     2170                                ( u(k+1,j,i)  + u(k,j,i)     )                &
     2171              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     2172                  +              ibit16 * adv_mom_3                           &
     2173                     ) *                                                      &
     2174                                ( u(k_pp,j,i) + u(k-1,j,i)   )                &
     2175              +      (           ibit17 * adv_mom_5                           &
     2176                     ) *                                                      &
     2177                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
     2178                                                  )
     2179
     2180          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
     2181                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     2182                  +     3.0_wp * ibit16 * adv_mom_3                           &
     2183                  +              ibit15 * adv_mom_1                           &
     2184                     ) *                                                      &
     2185                                ( u(k+1,j,i)   - u(k,j,i)    )                &
     2186              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     2187                  +              ibit16 * adv_mom_3                           &
     2188                     ) *                                                      &
     2189                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
     2190              +      (           ibit17 * adv_mom_5                           &
     2191                     ) *                                                      &
     2192                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
     2193                                                           )
     2194       ENDDO
     2195       
     2196       DO  k = nzb+1, nzt
    20532197
    20542198          flux_d    = flux_t(k-1)
     
    20572201!--       Calculate the divergence of the velocity field. A respective
    20582202!--       correction is needed to overcome numerical instabilities introduced
    2059 !--       by a not sufficient reduction of divergences near topography. 
     2203!--       by a not sufficient reduction of divergences near topography.
    20602204          div = ( ( u_comp(k)       * ( ibit9 + ibit10 + ibit11 )             &
    20612205                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
    2062                                     * ( IBITS(advc_flags_1(k,j,i-1),9,1)      &
    2063                                       + IBITS(advc_flags_1(k,j,i-1),10,1)     &
    2064                                       + IBITS(advc_flags_1(k,j,i-1),11,1)     &
     2206                                    * (                                       &
     2207                     REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )     &
     2208                   + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )     &
     2209                   + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )     &
    20652210                                      )                                       &
    20662211                  ) * ddx                                                     &
    20672212               +  ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 )         &
    20682213                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
    2069                                     * ( IBITS(advc_flags_1(k,j-1,i),12,1)     &
    2070                                       + IBITS(advc_flags_1(k,j-1,i),13,1)     &
    2071                                       + IBITS(advc_flags_1(k,j-1,i),14,1)     &
     2214                                    * (                                       &
     2215                     REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )     &
     2216                   + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )     &
     2217                   + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )     &
    20722218                                      )                                       &
    20732219                  ) * ddy                                                     &
    20742220               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
    20752221                -   w_comp(k-1) * rho_air_zw(k-1)                             &
    2076                                     * ( IBITS(advc_flags_1(k-1,j,i),15,1)     &
    2077                                       + IBITS(advc_flags_1(k-1,j,i),16,1)     &
    2078                                       + IBITS(advc_flags_1(k-1,j,i),17,1)     &
     2222                                    * (                                       &
     2223                     REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp )     &
     2224                   + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp )     &
     2225                   + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp )     &
    20792226                                      )                                       & 
    20802227                  ) * drho_air(k) * ddzw(k)                                   &
    20812228                ) * 0.5_wp
    20822229
    2083 
    2084           tend(k,j,i) = tend(k,j,i) - (                                       &
     2230           tend(k,j,i) = tend(k,j,i) - (                                      &
    20852231                            ( flux_r(k) + diss_r(k)                           &
    20862232                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
     
    20882234                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
    20892235                          + ( ( flux_t(k) + diss_t(k) )                       &
    2090                           -   ( flux_d    + diss_d )                          &
     2236                          -   ( flux_d    + diss_d    )                       &
    20912237                                                    ) * drho_air(k) * ddzw(k) &
    20922238                                       ) + div * u(k,j,i)
     
    21112257           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
    21122258                + ( flux_t(k)                                                  &
    2113                     * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                )     &
    2114                     / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        )     &
    2115                   + diss_t(k)                                                  &
    2116                     *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           )     &
    2117                     / ( ABS( w_comp(k) ) + 1.0E-20_wp                    )     &
    2118                   ) *   weight_substep(intermediate_timestep_count)
    2119        ENDDO
    2120 
    2121        DO  k = nzb_max+1, nzt
    2122 
    2123           u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2124           flux_r(k) = ( u_comp(k) - gu ) * (                                  &
    2125                          37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                &
    2126                        -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
    2127                        +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    2128           diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
    2129                          10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
    2130                        -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
    2131                        +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    2132 
    2133           v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2134           flux_n(k) = v_comp(k) * (                                           &
    2135                          37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
    2136                        -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
    2137                        +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    2138           diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    2139                          10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
    2140                        -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
    2141                        +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    2142 !
    2143 !--       k index has to be modified near bottom and top, else array
    2144 !--       subscripts will be exceeded.
    2145           ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
    2146           ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
    2147           ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
    2148 
    2149           k_ppp = k + 3 * ibit17
    2150           k_pp  = k + 2 * ( 1 - ibit15 )
    2151           k_mm  = k - 2 * ibit17
    2152 
    2153           w_comp(k) = w(k,j,i) + w(k,j,i-1)
    2154           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2155                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2156                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2157                   +              ibit15 * adv_mom_1                           &
    2158                      ) *                                                      &
    2159                                 ( u(k+1,j,i)  + u(k,j,i)     )                &
    2160               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2161                   +              ibit16 * adv_mom_3                           &
    2162                      ) *                                                      &
    2163                                 ( u(k_pp,j,i) + u(k-1,j,i)   )                &
    2164               +      (           ibit17 * adv_mom_5                           &
    2165                      ) *                                                      &
    2166                                 ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
    2167                                                   )
    2168 
    2169           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2170                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2171                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2172                   +              ibit15 * adv_mom_1                           &
    2173                      ) *                                                      &
    2174                                 ( u(k+1,j,i)   - u(k,j,i)    )                &
    2175               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2176                   +              ibit16 * adv_mom_3                           &
    2177                      ) *                                                      &
    2178                                 ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
    2179               +      (           ibit17 * adv_mom_5                           &
    2180                      ) *                                                      &
    2181                                 ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
    2182                                                            )
    2183 
    2184        ENDDO
    2185 
    2186        DO  k = nzb_max+1, nzt
    2187 
    2188           flux_d    = flux_t(k-1)
    2189           diss_d    = diss_t(k-1)
    2190 !
    2191 !--       Calculate the divergence of the velocity field. A respective
    2192 !--       correction is needed to overcome numerical instabilities introduced
    2193 !--       by a not sufficient reduction of divergences near topography.
    2194           div = ( ( u_comp(k)      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx    &
    2195                +  ( v_comp(k) + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy    &
    2196                +  ( w_comp(k)   * rho_air_zw(k)   -                           &
    2197                     w_comp(k-1) * rho_air_zw(k-1)                             &
    2198                   ) * drho_air(k) * ddzw(k)                                   &
    2199                 ) * 0.5_wp
    2200 
    2201           tend(k,j,i) = tend(k,j,i) - (                                       &
    2202                             ( flux_r(k) + diss_r(k)                           &
    2203                           -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
    2204                           + ( flux_n(k) + diss_n(k)                           &
    2205                           -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
    2206                           + ( ( flux_t(k) + diss_t(k) )                       &
    2207                           -   ( flux_d    + diss_d    )                       &
    2208                                                     ) * drho_air(k) * ddzw(k) &
    2209                                        ) + div * u(k,j,i)
    2210 
    2211            flux_l_u(k,j,tn) = flux_r(k)
    2212            diss_l_u(k,j,tn) = diss_r(k)
    2213            flux_s_u(k,tn)   = flux_n(k)
    2214            diss_s_u(k,tn)   = diss_n(k)
    2215 !
    2216 !--        Statistical Evaluation of u'u'. The factor has to be applied for
    2217 !--        right evaluation when gallilei_trans = .T. .
    2218            sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                           &
    2219                 + ( flux_r(k)                                                  &
    2220                     * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
    2221                     / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
    2222                   + diss_r(k)                                                  &
    2223                     *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
    2224                     / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
    2225                   ) *   weight_substep(intermediate_timestep_count)
    2226 !
    2227 !--        Statistical Evaluation of w'u'.
    2228            sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
    2229                 + ( flux_t(k)                                                  &
    22302259                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
    22312260                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
     
    22722301
    22732302       INTEGER(iwp)  ::  i      !< grid index along x-direction
    2274        INTEGER(iwp)  ::  ibit18 !< flag indicating 1st-order scheme along x-direction
    2275        INTEGER(iwp)  ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
    2276        INTEGER(iwp)  ::  ibit20 !< flag indicating 5th-order scheme along x-direction
    2277        INTEGER(iwp)  ::  ibit21 !< flag indicating 1st-order scheme along y-direction
    2278        INTEGER(iwp)  ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
    2279        INTEGER(iwp)  ::  ibit23 !< flag indicating 3rd-order scheme along y-direction
    2280        INTEGER(iwp)  ::  ibit24 !< flag indicating 1st-order scheme along z-direction
    2281        INTEGER(iwp)  ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
    2282        INTEGER(iwp)  ::  ibit26 !< flag indicating 3rd-order scheme along z-direction
    22832303       INTEGER(iwp)  ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
    22842304       INTEGER(iwp)  ::  j      !< grid index along y-direction
     
    22892309       INTEGER(iwp)  ::  tn     !< number of OpenMP thread
    22902310       
    2291        REAL(wp)     ::  diss_d   !< artificial dissipation term at grid box bottom
    2292        REAL(wp)     ::  div      !< divergence on v-grid
    2293        REAL(wp)     ::  flux_d   !< 6th-order flux at grid box bottom
    2294        REAL(wp)     ::  gu       !< Galilei-transformation velocity along x
    2295        REAL(wp)     ::  gv       !< Galilei-transformation velocity along y
    2296        REAL(wp)     ::  v_comp_l !< advection velocity along y on leftmost grid point on subdomain
     2311       REAL(wp)      ::  ibit18   !< flag indicating 1st-order scheme along x-direction
     2312       REAL(wp)      ::  ibit19   !< flag indicating 3rd-order scheme along x-direction
     2313       REAL(wp)      ::  ibit20   !< flag indicating 5th-order scheme along x-direction
     2314       REAL(wp)      ::  ibit21   !< flag indicating 1st-order scheme along y-direction
     2315       REAL(wp)      ::  ibit22   !< flag indicating 3rd-order scheme along y-direction
     2316       REAL(wp)      ::  ibit23   !< flag indicating 3rd-order scheme along y-direction
     2317       REAL(wp)      ::  ibit24   !< flag indicating 1st-order scheme along z-direction
     2318       REAL(wp)      ::  ibit25   !< flag indicating 3rd-order scheme along z-direction
     2319       REAL(wp)      ::  ibit26   !< flag indicating 3rd-order scheme along z-direction
     2320       REAL(wp)      ::  diss_d   !< artificial dissipation term at grid box bottom
     2321       REAL(wp)      ::  div      !< divergence on v-grid
     2322       REAL(wp)      ::  flux_d   !< 6th-order flux at grid box bottom
     2323       REAL(wp)      ::  gu       !< Galilei-transformation velocity along x
     2324       REAL(wp)      ::  gv       !< Galilei-transformation velocity along y
     2325       REAL(wp)      ::  v_comp_l !< advection velocity along y on leftmost grid point on subdomain
    22972326       
    22982327       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
     
    23152344          DO  k = nzb+1, nzb_max
    23162345
    2317              ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
    2318              ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
    2319              ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
     2346             ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
     2347             ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )
     2348             ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )
    23202349
    23212350             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
     
    23732402          DO  k = nzb+1, nzb_max
    23742403
    2375              ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
    2376              ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
    2377              ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
     2404             ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
     2405             ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )
     2406             ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
    23782407
    23792408             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     
    24252454         
    24262455       ENDIF
    2427 
    2428        flux_t(0) = 0.0_wp
    2429        diss_t(0) = 0.0_wp
    2430        w_comp(0) = 0.0_wp
    24312456!
    24322457!--    Now compute the fluxes and tendency terms for the horizontal and
     
    24342459       DO  k = nzb+1, nzb_max
    24352460
    2436           ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
    2437           ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
    2438           ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
     2461          ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
     2462          ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
     2463          ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
    24392464 
    24402465          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
     
    24692494                                           )
    24702495
    2471           ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
    2472           ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
    2473           ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
     2496          ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
     2497          ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
     2498          ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
    24742499
    24752500
     
    25042529                                    ( v(k,j+3,i) - v(k,j-2,i) )               &
    25052530                                                )
     2531       ENDDO
     2532
     2533       DO  k = nzb_max+1, nzt
     2534
     2535          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
     2536          flux_r(k) = u_comp(k) * (                                           &
     2537                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
     2538                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
     2539                    +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
     2540
     2541          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
     2542                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
     2543                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
     2544                    +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     2545
     2546
     2547          v_comp(k) = v(k,j+1,i) + v(k,j,i)
     2548          flux_n(k) = ( v_comp(k) - gv ) * (                                  &
     2549                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                   &
     2550                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                   &
     2551                      +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
     2552
     2553          diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
     2554                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                   &
     2555                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                   &
     2556                    +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
     2557       ENDDO
     2558!
     2559!--    Now, compute vertical fluxes. Split loop into a part treating the
     2560!--    lowest 2 grid points with indirect indexing, a main loop without
     2561!--    indirect indexing, and a loop for the uppermost 2 grip points with
     2562!--    indirect indexing. This allows better vectorization for the main loop.
     2563!--    First, compute the flux at model surface, which need has to be
     2564!--    calculated explicetely for the tendency at
     2565!--    the first w-level. For topography wall this is done implicitely by
     2566!--    advc_flags_1.
     2567       flux_t(nzb) = 0.0_wp
     2568       diss_t(nzb) = 0.0_wp
     2569       w_comp(nzb) = 0.0_wp
     2570       
     2571       DO  k = nzb+1, nzb+2
    25062572!
    25072573!--       k index has to be modified near bottom and top, else array
    25082574!--       subscripts will be exceeded.
    2509           ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
    2510           ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
    2511           ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
     2575          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
     2576          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
     2577          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
    25122578
    25132579          k_ppp = k + 3 * ibit26
     
    25462612                                                           )
    25472613       ENDDO
    2548 
    2549        DO  k = nzb+1, nzb_max
     2614       
     2615       DO  k = nzb+3, nzt-2
     2616
     2617          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
     2618          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
     2619          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
     2620
     2621          w_comp(k) = w(k,j-1,i) + w(k,j,i)
     2622          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
     2623                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     2624                  +     7.0_wp * ibit25 * adv_mom_3                           &
     2625                  +              ibit24 * adv_mom_1                           &
     2626                     ) *                                                      &
     2627                                ( v(k+1,j,i)   + v(k,j,i)    )                &
     2628              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     2629                  +              ibit25 * adv_mom_3                           &
     2630                     ) *                                                      &
     2631                                ( v(k+2,j,i)  + v(k-1,j,i)  )                 &
     2632              +      (           ibit26 * adv_mom_5                           &
     2633                     ) *                                                      &
     2634                                ( v(k+3,j,i) + v(k-2,j,i) )                   &
     2635                                                  )
     2636
     2637          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
     2638                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     2639                  +     3.0_wp * ibit25 * adv_mom_3                           &
     2640                  +              ibit24 * adv_mom_1                           &
     2641                     ) *                                                      &
     2642                                ( v(k+1,j,i)   - v(k,j,i)    )                &
     2643              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     2644                  +              ibit25 * adv_mom_3                           &
     2645                     ) *                                                      &
     2646                                ( v(k+2,j,i)  - v(k-1,j,i)  )                 &
     2647              +      (           ibit26 * adv_mom_5                           &
     2648                     ) *                                                      &
     2649                                ( v(k+3,j,i) - v(k-2,j,i) )                   &
     2650                                                           )
     2651       ENDDO
     2652       
     2653       DO  k = nzt-1, nzt
     2654!
     2655!--       k index has to be modified near bottom and top, else array
     2656!--       subscripts will be exceeded.
     2657          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
     2658          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
     2659          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
     2660
     2661          k_ppp = k + 3 * ibit26
     2662          k_pp  = k + 2 * ( 1 - ibit24  )
     2663          k_mm  = k - 2 * ibit26
     2664
     2665          w_comp(k) = w(k,j-1,i) + w(k,j,i)
     2666          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
     2667                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     2668                  +     7.0_wp * ibit25 * adv_mom_3                           &
     2669                  +              ibit24 * adv_mom_1                           &
     2670                     ) *                                                      &
     2671                                ( v(k+1,j,i)   + v(k,j,i)    )                &
     2672              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     2673                  +              ibit25 * adv_mom_3                           &
     2674                     ) *                                                      &
     2675                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
     2676              +      (           ibit26 * adv_mom_5                           &
     2677                     ) *                                                      &
     2678                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
     2679                                                  )
     2680
     2681          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
     2682                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     2683                  +     3.0_wp * ibit25 * adv_mom_3                           &
     2684                  +              ibit24 * adv_mom_1                           &
     2685                     ) *                                                      &
     2686                                ( v(k+1,j,i)   - v(k,j,i)    )                &
     2687              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     2688                  +              ibit25 * adv_mom_3                           &
     2689                     ) *                                                      &
     2690                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
     2691              +      (           ibit26 * adv_mom_5                           &
     2692                     ) *                                                      &
     2693                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
     2694                                                           )
     2695       ENDDO
     2696       
     2697       DO  k = nzb+1, nzt
    25502698
    25512699          flux_d    = flux_t(k-1)
     
    25542702!--       Calculate the divergence of the velocity field. A respective
    25552703!--       correction is needed to overcome numerical instabilities introduced
    2556 !--       by a not sufficient reduction of divergences near topography. 
     2704!--       by a not sufficient reduction of divergences near topography.
    25572705          div = ( ( ( u_comp(k)     + gu )                                    &
    25582706                                       * ( ibit18 + ibit19 + ibit20 )         &
    25592707                  - ( u(k,j-1,i) + u(k,j,i) )                                 &
    2560                                        * ( IBITS(advc_flags_1(k,j,i-1),18,1)  &
    2561                                          + IBITS(advc_flags_1(k,j,i-1),19,1)  &
    2562                                          + IBITS(advc_flags_1(k,j,i-1),20,1)  &
     2708                                       * (                                    &
     2709                        REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )  &
     2710                      + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )  &
     2711                      + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )  &
    25632712                                         )                                    &
    25642713                  ) * ddx                                                     &
     
    25662715                                       * ( ibit21 + ibit22 + ibit23 )         &
    25672716                - ( v(k,j,i)     + v(k,j-1,i) )                               &
    2568                                        * ( IBITS(advc_flags_1(k,j-1,i),21,1)  &
    2569                                          + IBITS(advc_flags_1(k,j-1,i),22,1)  &
    2570                                          + IBITS(advc_flags_1(k,j-1,i),23,1)  &
     2717                                       * (                                    &
     2718                        REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )  &
     2719                      + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )  &
     2720                      + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )  &
    25712721                                         )                                    &
    25722722                  ) * ddy                                                     &
    25732723               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&
    25742724                -   w_comp(k-1) * rho_air_zw(k-1)                             &
    2575                                        * ( IBITS(advc_flags_1(k-1,j,i),24,1)  &
    2576                                          + IBITS(advc_flags_1(k-1,j,i),25,1)  &
    2577                                          + IBITS(advc_flags_1(k-1,j,i),26,1)  &
     2725                                       * (                                    &
     2726                        REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp )  &
     2727                      + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp )  &
     2728                      + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp )  &
    25782729                                         )                                    &
    25792730                   ) * drho_air(k) * ddzw(k)                                  &
    25802731                ) * 0.5_wp
    2581 
    25822732
    25832733          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    25952745           flux_s_v(k,tn)   = flux_n(k)
    25962746           diss_s_v(k,tn)   = diss_n(k)
    2597 
    25982747!
    25992748!--        Statistical Evaluation of v'v'. The factor has to be applied for
     
    26192768
    26202769       ENDDO
    2621 
    2622        DO  k = nzb_max+1, nzt
    2623 
    2624           u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
    2625           flux_r(k) = u_comp(k) * (                                           &
    2626                       37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
    2627                     -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
    2628                     +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    2629 
    2630           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    2631                       10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
    2632                     -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
    2633                     +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    2634 
    2635 
    2636           v_comp(k) = v(k,j+1,i) + v(k,j,i)
    2637           flux_n(k) = ( v_comp(k) - gv ) * (                                  &
    2638                       37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                   &
    2639                     -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                   &
    2640                       +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    2641 
    2642           diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
    2643                       10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                   &
    2644                     -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                   &
    2645                     +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    2646 !
    2647 !--       k index has to be modified near bottom and top, else array
    2648 !--       subscripts will be exceeded.
    2649           ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
    2650           ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
    2651           ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
    2652 
    2653           k_ppp = k + 3 * ibit26
    2654           k_pp  = k + 2 * ( 1 - ibit24  )
    2655           k_mm  = k - 2 * ibit26
    2656 
    2657           w_comp(k) = w(k,j-1,i) + w(k,j,i)
    2658           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2659                      ( 37.0_wp * ibit26 * adv_mom_5                           &
    2660                   +     7.0_wp * ibit25 * adv_mom_3                           &
    2661                   +              ibit24 * adv_mom_1                           &
    2662                      ) *                                                      &
    2663                                 ( v(k+1,j,i)   + v(k,j,i)    )                &
    2664               -      (  8.0_wp * ibit26 * adv_mom_5                           &
    2665                   +              ibit25 * adv_mom_3                           &
    2666                      ) *                                                      &
    2667                                 ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
    2668               +      (           ibit26 * adv_mom_5                           &
    2669                      ) *                                                      &
    2670                                 ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
    2671                                                   )
    2672 
    2673           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2674                      ( 10.0_wp * ibit26 * adv_mom_5                           &
    2675                   +     3.0_wp * ibit25 * adv_mom_3                           &
    2676                   +              ibit24 * adv_mom_1                           &
    2677                      ) *                                                      &
    2678                                 ( v(k+1,j,i)   - v(k,j,i)    )                &
    2679               -      (  5.0_wp * ibit26 * adv_mom_5                           &
    2680                   +              ibit25 * adv_mom_3                           &
    2681                      ) *                                                      &
    2682                                 ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
    2683               +      (           ibit26 * adv_mom_5                           &
    2684                      ) *                                                      &
    2685                                 ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
    2686                                                            )
    2687        ENDDO
    2688 
    2689        DO  k = nzb_max+1, nzt
    2690 
    2691           flux_d    = flux_t(k-1)
    2692           diss_d    = diss_t(k-1)
    2693 !
    2694 !--       Calculate the divergence of the velocity field. A respective
    2695 !--       correction is needed to overcome numerical instabilities introduced
    2696 !--       by a not sufficient reduction of divergences near topography.
    2697           div = ( ( u_comp(k) + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx    &
    2698                +  ( v_comp(k)      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy    &
    2699                +  ( w_comp(k)   * rho_air_zw(k)   -                           &
    2700                     w_comp(k-1) * rho_air_zw(k-1)                             &
    2701                   ) * drho_air(k) * ddzw(k)                                   &
    2702                 ) * 0.5_wp
    2703 
    2704           tend(k,j,i) = tend(k,j,i) - (                                       &
    2705                          ( flux_r(k) + diss_r(k)                              &
    2706                        -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
    2707                        + ( flux_n(k) + diss_n(k)                              &
    2708                        -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
    2709                        + ( ( flux_t(k) + diss_t(k) )                          &
    2710                        -   ( flux_d    + diss_d    )                          &
    2711                                                    ) * drho_air(k) * ddzw(k)  &
    2712                                       ) + v(k,j,i) * div
    2713 
    2714            flux_l_v(k,j,tn) = flux_r(k)
    2715            diss_l_v(k,j,tn) = diss_r(k)
    2716            flux_s_v(k,tn)   = flux_n(k)
    2717            diss_s_v(k,tn)   = diss_n(k)
    2718 !
    2719 !--        Statistical Evaluation of v'v'. The factor has to be applied for
    2720 !--        right evaluation when gallilei_trans = .T. .
    2721            sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                           &
    2722                 + ( flux_n(k)                                                  &
    2723                     * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
    2724                     / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
    2725                   + diss_n(k)                                                  &
    2726                     *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
    2727                     / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
    2728                   ) *   weight_substep(intermediate_timestep_count)
    2729 !
    2730 !--        Statistical Evaluation of w'u'.
    2731            sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
    2732                 + ( flux_t(k)                                                  &
    2733                     * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
    2734                     / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
    2735                   + diss_t(k)                                                  &
    2736                     *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
    2737                     / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
    2738                   ) *   weight_substep(intermediate_timestep_count)
    2739 
    2740        ENDDO
    27412770       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
    27422771
     
    27742803
    27752804       INTEGER(iwp) ::  i      !< grid index along x-direction
    2776        INTEGER(iwp) ::  ibit27 !< flag indicating 1st-order scheme along x-direction
    2777        INTEGER(iwp) ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
    2778        INTEGER(iwp) ::  ibit29 !< flag indicating 5th-order scheme along x-direction
    2779        INTEGER(iwp) ::  ibit30 !< flag indicating 1st-order scheme along y-direction
    2780        INTEGER(iwp) ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
    2781        INTEGER(iwp) ::  ibit32 !< flag indicating 5th-order scheme along y-direction
    2782        INTEGER(iwp) ::  ibit33 !< flag indicating 1st-order scheme along z-direction
    2783        INTEGER(iwp) ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
    2784        INTEGER(iwp) ::  ibit35 !< flag indicating 5th-order scheme along z-direction
    27852805       INTEGER(iwp) ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
    27862806       INTEGER(iwp) ::  j      !< grid index along y-direction
     
    27912811       INTEGER(iwp) ::  tn     !< number of OpenMP thread
    27922812       
     2813       REAL(wp)    ::  ibit27  !< flag indicating 1st-order scheme along x-direction
     2814       REAL(wp)    ::  ibit28  !< flag indicating 3rd-order scheme along x-direction
     2815       REAL(wp)    ::  ibit29  !< flag indicating 5th-order scheme along x-direction
     2816       REAL(wp)    ::  ibit30  !< flag indicating 1st-order scheme along y-direction
     2817       REAL(wp)    ::  ibit31  !< flag indicating 3rd-order scheme along y-direction
     2818       REAL(wp)    ::  ibit32  !< flag indicating 5th-order scheme along y-direction
     2819       REAL(wp)    ::  ibit33  !< flag indicating 1st-order scheme along z-direction
     2820       REAL(wp)    ::  ibit34  !< flag indicating 3rd-order scheme along z-direction
     2821       REAL(wp)    ::  ibit35  !< flag indicating 5th-order scheme along z-direction
    27932822       REAL(wp)    ::  diss_d  !< discretized artificial dissipation at top of the grid box
    27942823       REAL(wp)    ::  div     !< divergence on w-grid
     
    28152844
    28162845          DO  k = nzb+1, nzb_max
    2817              ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
    2818              ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
    2819              ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
     2846             ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
     2847             ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
     2848             ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )
    28202849
    28212850             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
     
    28592888                         -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )              &
    28602889                         +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    2861              diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                         &
     2890             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
    28622891                           10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )              &
    28632892                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )              &
     
    28732902          DO  k = nzb+1, nzb_max
    28742903
    2875              ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
    2876              ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
    2877              ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
     2904             ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
     2905             ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )
     2906             ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )
    28782907
    28792908             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
     
    29262955       ENDIF
    29272956!
    2928 !--    The lower flux has to be calculated explicetely for the tendency at
     2957!--    Now compute the fluxes and tendency terms for the horizontal
     2958!--    and vertical parts.
     2959       DO  k = nzb+1, nzb_max
     2960
     2961          ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
     2962          ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
     2963          ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
     2964
     2965          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
     2966          flux_r(k) = u_comp(k) * (                                           &
     2967                     ( 37.0_wp * ibit29 * adv_mom_5                           &
     2968                  +     7.0_wp * ibit28 * adv_mom_3                           &
     2969                  +              ibit27 * adv_mom_1                           &
     2970                     ) *                                                      &
     2971                                    ( w(k,j,i+1) + w(k,j,i)   )               &
     2972              -      (  8.0_wp * ibit29 * adv_mom_5                           &
     2973                  +              ibit28 * adv_mom_3                           &
     2974                     ) *                                                      &
     2975                                    ( w(k,j,i+2) + w(k,j,i-1) )               &
     2976              +      (           ibit29 * adv_mom_5                           &
     2977                     ) *                                                      &
     2978                                    ( w(k,j,i+3) + w(k,j,i-2) )               &
     2979                                  )
     2980
     2981          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
     2982                     ( 10.0_wp * ibit29 * adv_mom_5                           &
     2983                  +     3.0_wp * ibit28 * adv_mom_3                           &
     2984                  +              ibit27 * adv_mom_1                           &
     2985                     ) *                                                      &
     2986                                    ( w(k,j,i+1) - w(k,j,i)   )               &
     2987              -      (  5.0_wp * ibit29 * adv_mom_5                           &
     2988                  +              ibit28 * adv_mom_3                           &
     2989                     ) *                                                      &
     2990                                    ( w(k,j,i+2) - w(k,j,i-1) )               &
     2991              +      (           ibit29 * adv_mom_5                           &
     2992                     ) *                                                      &
     2993                                    ( w(k,j,i+3) - w(k,j,i-2) )               &
     2994                                           )
     2995
     2996          ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),  KIND = wp )
     2997          ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
     2998          ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
     2999
     3000          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
     3001          flux_n(k) = v_comp(k) * (                                           &
     3002                     ( 37.0_wp * ibit32 * adv_mom_5                           &
     3003                  +     7.0_wp * ibit31 * adv_mom_3                           &
     3004                  +              ibit30 * adv_mom_1                           &
     3005                     ) *                                                      &
     3006                                    ( w(k,j+1,i) + w(k,j,i)   )               &
     3007              -      (  8.0_wp * ibit32 * adv_mom_5                           &
     3008                  +              ibit31 * adv_mom_3                           &
     3009                     ) *                                                      &
     3010                                    ( w(k,j+2,i) + w(k,j-1,i) )               &
     3011              +      (           ibit32 * adv_mom_5                           &
     3012                     ) *                                                      &
     3013                                    ( w(k,j+3,i) + w(k,j-2,i) )               &
     3014                                  )
     3015
     3016          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
     3017                     ( 10.0_wp * ibit32 * adv_mom_5                           &
     3018                  +     3.0_wp * ibit31 * adv_mom_3                           &
     3019                  +              ibit30 * adv_mom_1                           &
     3020                     ) *                                                      &
     3021                                    ( w(k,j+1,i) - w(k,j,i)  )                &
     3022              -      (  5.0_wp * ibit32 * adv_mom_5                           &
     3023                  +              ibit31 * adv_mom_3                           &
     3024                     ) *                                                      &
     3025                                   ( w(k,j+2,i) - w(k,j-1,i) )                &
     3026              +      (           ibit32 * adv_mom_5                           &
     3027                     ) *                                                      &
     3028                                   ( w(k,j+3,i) - w(k,j-2,i) )                &
     3029                                           )
     3030       ENDDO
     3031
     3032       DO  k = nzb_max+1, nzt
     3033
     3034          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
     3035          flux_r(k) = u_comp(k) * (                                           &
     3036                      37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
     3037                    -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
     3038                    +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
     3039
     3040          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
     3041                      10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
     3042                    -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
     3043                    +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     3044
     3045          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
     3046          flux_n(k) = v_comp(k) * (                                           &
     3047                      37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
     3048                    -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
     3049                    +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
     3050
     3051          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
     3052                      10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
     3053                    -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
     3054                    +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     3055       ENDDO
     3056
     3057!
     3058!--    Now, compute vertical fluxes. Split loop into a part treating the
     3059!--    lowest 2 grid points with indirect indexing, a main loop without
     3060!--    indirect indexing, and a loop for the uppermost 2 grip points with
     3061!--    indirect indexing. This allows better vectorization for the main loop.
     3062!--    First, compute the flux at model surface, which need has to be
     3063!--    calculated explicetely for the tendency at
    29293064!--    the first w-level. For topography wall this is done implicitely by
    29303065!--    advc_flags_1.
     
    29333068       flux_t(0) = w_comp(k)       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
    29343069       diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
    2935 !
    2936 !--    Now compute the fluxes and tendency terms for the horizontal
    2937 !--    and vertical parts.
    2938        DO  k = nzb+1, nzb_max
    2939 
    2940           ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
    2941           ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
    2942           ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
    2943 
    2944           u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
    2945           flux_r(k) = u_comp(k) * (                                           &
    2946                      ( 37.0_wp * ibit29 * adv_mom_5                           &
    2947                   +     7.0_wp * ibit28 * adv_mom_3                           &
    2948                   +              ibit27 * adv_mom_1                           &
    2949                      ) *                                                      &
    2950                                     ( w(k,j,i+1) + w(k,j,i)   )               &
    2951               -      (  8.0_wp * ibit29 * adv_mom_5                           &
    2952                   +              ibit28 * adv_mom_3                           &
    2953                      ) *                                                      &
    2954                                     ( w(k,j,i+2) + w(k,j,i-1) )               &
    2955               +      (           ibit29 * adv_mom_5                           &
    2956                      ) *                                                      &
    2957                                     ( w(k,j,i+3) + w(k,j,i-2) )               &
    2958                                   )
    2959 
    2960           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    2961                      ( 10.0_wp * ibit29 * adv_mom_5                           &
    2962                   +     3.0_wp * ibit28 * adv_mom_3                           &
    2963                   +              ibit27 * adv_mom_1                           &
    2964                      ) *                                                      &
    2965                                     ( w(k,j,i+1) - w(k,j,i)   )               &
    2966               -      (  5.0_wp * ibit29 * adv_mom_5                           &
    2967                   +              ibit28 * adv_mom_3                           &
    2968                      ) *                                                      &
    2969                                     ( w(k,j,i+2) - w(k,j,i-1) )               &
    2970               +      (           ibit29 * adv_mom_5                           &
    2971                      ) *                                                      &
    2972                                     ( w(k,j,i+3) - w(k,j,i-2) )               &
    2973                                            )
    2974 
    2975           ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
    2976           ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
    2977           ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
    2978 
    2979           v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
    2980           flux_n(k) = v_comp(k) * (                                           &
    2981                      ( 37.0_wp * ibit32 * adv_mom_5                           &
    2982                   +     7.0_wp * ibit31 * adv_mom_3                           &
    2983                   +              ibit30 * adv_mom_1                           &
    2984                      ) *                                                      &
    2985                                     ( w(k,j+1,i) + w(k,j,i)   )               &
    2986               -      (  8.0_wp * ibit32 * adv_mom_5                           &
    2987                   +              ibit31 * adv_mom_3                           &
    2988                      ) *                                                      &
    2989                                     ( w(k,j+2,i) + w(k,j-1,i) )               &
    2990               +      (           ibit32 * adv_mom_5                           &
    2991                      ) *                                                      &
    2992                                     ( w(k,j+3,i) + w(k,j-2,i) )               &
    2993                                   )
    2994 
    2995           diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    2996                      ( 10.0_wp * ibit32 * adv_mom_5                           &
    2997                   +     3.0_wp * ibit31 * adv_mom_3                           &
    2998                   +              ibit30 * adv_mom_1                           &
    2999                      ) *                                                      &
    3000                                     ( w(k,j+1,i) - w(k,j,i)  )                &
    3001               -      (  5.0_wp * ibit32 * adv_mom_5                           &
    3002                   +              ibit31 * adv_mom_3                           &
    3003                      ) *                                                      &
    3004                                    ( w(k,j+2,i) - w(k,j-1,i) )                &
    3005               +      (           ibit32 * adv_mom_5                           &
    3006                      ) *                                                      &
    3007                                    ( w(k,j+3,i) - w(k,j-2,i) )                &
    3008                                            )
     3070       
     3071       DO  k = nzb+1, nzb+2
    30093072!
    30103073!--       k index has to be modified near bottom and top, else array
    30113074!--       subscripts will be exceeded.
    3012           ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
    3013           ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
    3014           ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
     3075          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
     3076          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
     3077          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
    30153078
    30163079          k_ppp = k + 3 * ibit35
     
    30493112                                                          )
    30503113       ENDDO
    3051 
    3052        DO  k = nzb+1, nzb_max
     3114       
     3115       DO  k = nzb+3, nzt-2
     3116       
     3117          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
     3118          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
     3119          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
     3120
     3121          w_comp(k) = w(k+1,j,i) + w(k,j,i)
     3122          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
     3123                     ( 37.0_wp * ibit35 * adv_mom_5                           &
     3124                  +     7.0_wp * ibit34 * adv_mom_3                           &
     3125                  +              ibit33 * adv_mom_1                           &
     3126                     ) *                                                      &
     3127                                ( w(k+1,j,i)  + w(k,j,i)     )                &
     3128              -      (  8.0_wp * ibit35 * adv_mom_5                           &
     3129                  +              ibit34 * adv_mom_3                           &
     3130                     ) *                                                      &
     3131                                ( w(k+2,j,i)  + w(k-1,j,i)  )                 &
     3132              +      (           ibit35 * adv_mom_5                           &
     3133                     ) *                                                      &
     3134                                ( w(k+3,j,i) + w(k-2,j,i) )                   &
     3135                                                 )
     3136
     3137          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
     3138                     ( 10.0_wp * ibit35 * adv_mom_5                           &
     3139                  +     3.0_wp * ibit34 * adv_mom_3                           &
     3140                  +              ibit33 * adv_mom_1                           &
     3141                     ) *                                                      &
     3142                                ( w(k+1,j,i)   - w(k,j,i)    )                &
     3143              -      (  5.0_wp * ibit35 * adv_mom_5                           &
     3144                  +              ibit34 * adv_mom_3                           &
     3145                     ) *                                                      &
     3146                                ( w(k+2,j,i)  - w(k-1,j,i)  )                 &
     3147              +      (           ibit35 * adv_mom_5                           &
     3148                     ) *                                                      &
     3149                                ( w(k+3,j,i) - w(k-2,j,i) )                   &
     3150                                                          )
     3151       ENDDO
     3152       
     3153       DO  k = nzt-1, nzt
     3154!
     3155!--       k index has to be modified near bottom and top, else array
     3156!--       subscripts will be exceeded.
     3157          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
     3158          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
     3159          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
     3160
     3161          k_ppp = k + 3 * ibit35
     3162          k_pp  = k + 2 * ( 1 - ibit33  )
     3163          k_mm  = k - 2 * ibit35
     3164
     3165          w_comp(k) = w(k+1,j,i) + w(k,j,i)
     3166          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
     3167                     ( 37.0_wp * ibit35 * adv_mom_5                           &
     3168                  +     7.0_wp * ibit34 * adv_mom_3                           &
     3169                  +              ibit33 * adv_mom_1                           &
     3170                     ) *                                                      &
     3171                                ( w(k+1,j,i)  + w(k,j,i)     )                &
     3172              -      (  8.0_wp * ibit35 * adv_mom_5                           &
     3173                  +              ibit34 * adv_mom_3                           &
     3174                     ) *                                                      &
     3175                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
     3176              +      (           ibit35 * adv_mom_5                           &
     3177                     ) *                                                      &
     3178                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
     3179                                                 )
     3180
     3181          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
     3182                     ( 10.0_wp * ibit35 * adv_mom_5                           &
     3183                  +     3.0_wp * ibit34 * adv_mom_3                           &
     3184                  +              ibit33 * adv_mom_1                           &
     3185                     ) *                                                      &
     3186                                ( w(k+1,j,i)   - w(k,j,i)    )                &
     3187              -      (  5.0_wp * ibit35 * adv_mom_5                           &
     3188                  +              ibit34 * adv_mom_3                           &
     3189                     ) *                                                      &
     3190                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
     3191              +      (           ibit35 * adv_mom_5                           &
     3192                     ) *                                                      &
     3193                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
     3194                                                          )
     3195       ENDDO
     3196       
     3197       DO  k = nzb+1, nzt
    30533198
    30543199          flux_d    = flux_t(k-1)
     
    30573202!--       Calculate the divergence of the velocity field. A respective
    30583203!--       correction is needed to overcome numerical instabilities introduced
    3059 !--       by a not sufficient reduction of divergences near topography. 
     3204!--       by a not sufficient reduction of divergences near topography.
    30603205          div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 )         &
    30613206                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
    3062                                     * ( IBITS(advc_flags_1(k,j,i-1),27,1)     &
    3063                                       + IBITS(advc_flags_1(k,j,i-1),28,1)     &
    3064                                       + IBITS(advc_flags_1(k,j,i-1),29,1)     &
     3207                                    * (                                       &
     3208                     REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )     &
     3209                   + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )     &
     3210                   + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )     &
    30653211                                      )                                       &
    30663212                  ) * ddx                                                     &
    30673213              +   ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 )         &
    30683214                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
    3069                                     * ( IBITS(advc_flags_1(k,j-1,i),30,1)     &
    3070                                       + IBITS(advc_flags_1(k,j-1,i),31,1)     &
    3071                                       + IBITS(advc_flags_2(k,j-1,i),0,1)      &
     3215                                    * (                                       &
     3216                     REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )     &
     3217                   + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )     &
     3218                   + REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )     &
    30723219                                      )                                       &
    30733220                  ) * ddy                                                     &
     
    30753222                                            * ( ibit33 + ibit34 + ibit35 )    &
    30763223                - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
    3077                                     * ( IBITS(advc_flags_2(k-1,j,i),1,1)      &
    3078                                       + IBITS(advc_flags_2(k-1,j,i),2,1)      &
    3079                                       + IBITS(advc_flags_2(k-1,j,i),3,1)      &
     3224                                    * (                                       &
     3225                     REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp )      &
     3226                   + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp )      &
     3227                   + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp )      &
    30803228                                      )                                       &
    3081                   ) * drho_air_zw(k) * ddzu(k+1)                              &
    3082                 ) * 0.5_wp
    3083 
    3084 
    3085           tend(k,j,i) = tend(k,j,i) - (                                       &
    3086                       ( flux_r(k) + diss_r(k)                                 &
    3087                     -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
    3088                     + ( flux_n(k) + diss_n(k)                                 &
    3089                     -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
    3090                     + ( ( flux_t(k) + diss_t(k) )                             &
    3091                     -   ( flux_d    + diss_d    )                             &
    3092                                               ) * drho_air_zw(k) * ddzu(k+1)  &
    3093                                       ) + div * w(k,j,i)
    3094 
    3095           flux_l_w(k,j,tn) = flux_r(k)
    3096           diss_l_w(k,j,tn) = diss_r(k)
    3097           flux_s_w(k,tn)   = flux_n(k)
    3098           diss_s_w(k,tn)   = diss_n(k)
    3099 !
    3100 !--       Statistical Evaluation of w'w'.
    3101           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
    3102                       + ( flux_t(k)                                           &
    3103                        * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                ) &
    3104                        / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        ) &
    3105                         + diss_t(k)                                           &
    3106                        *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           ) &
    3107                        / ( ABS( w_comp(k) ) + 1.0E-20_wp                    ) &
    3108                         ) *   weight_substep(intermediate_timestep_count)
    3109 
    3110        ENDDO
    3111 
    3112        DO  k = nzb_max+1, nzt
    3113 
    3114           u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
    3115           flux_r(k) = u_comp(k) * (                                           &
    3116                       37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
    3117                     -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
    3118                     +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    3119 
    3120           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    3121                       10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
    3122                     -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
    3123                     +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    3124 
    3125           v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
    3126           flux_n(k) = v_comp(k) * (                                           &
    3127                       37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
    3128                     -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
    3129                     +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    3130 
    3131           diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    3132                       10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
    3133                     -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
    3134                     +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    3135 !
    3136 !--       k index has to be modified near bottom and top, else array
    3137 !--       subscripts will be exceeded.
    3138           ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
    3139           ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
    3140           ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
    3141 
    3142           k_ppp = k + 3 * ibit35
    3143           k_pp  = k + 2 * ( 1 - ibit33  )
    3144           k_mm  = k - 2 * ibit35
    3145 
    3146           w_comp(k) = w(k+1,j,i) + w(k,j,i)
    3147           flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    3148                      ( 37.0_wp * ibit35 * adv_mom_5                           &
    3149                   +     7.0_wp * ibit34 * adv_mom_3                           &
    3150                   +              ibit33 * adv_mom_1                           &
    3151                      ) *                                                      &
    3152                                 ( w(k+1,j,i)  + w(k,j,i)     )                &
    3153               -      (  8.0_wp * ibit35 * adv_mom_5                           &
    3154                   +              ibit34 * adv_mom_3                           &
    3155                      ) *                                                      &
    3156                                 ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
    3157               +      (           ibit35 * adv_mom_5                           &
    3158                      ) *                                                      &
    3159                                 ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
    3160                                                  )
    3161 
    3162           diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    3163                      ( 10.0_wp * ibit35 * adv_mom_5                           &
    3164                   +     3.0_wp * ibit34 * adv_mom_3                           &
    3165                   +              ibit33 * adv_mom_1                           &
    3166                      ) *                                                      &
    3167                                 ( w(k+1,j,i)   - w(k,j,i)    )                &
    3168               -      (  5.0_wp * ibit35 * adv_mom_5                           &
    3169                   +              ibit34 * adv_mom_3                           &
    3170                      ) *                                                      &
    3171                                 ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
    3172               +      (           ibit35 * adv_mom_5                           &
    3173                      ) *                                                      &
    3174                                 ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
    3175                                                           )
    3176        ENDDO
    3177 
    3178        DO  k = nzb_max+1, nzt
    3179 
    3180           flux_d    = flux_t(k-1)
    3181           diss_d    = diss_t(k-1)
    3182 !
    3183 !--       Calculate the divergence of the velocity field. A respective
    3184 !--       correction is needed to overcome numerical instabilities introduced
    3185 !--       by a not sufficient reduction of divergences near topography.
    3186           div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx      &
    3187               +   ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy      &
    3188               +   ( w_comp(k)               * rho_air(k+1) -                  &
    3189                   ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
    31903229                  ) * drho_air_zw(k) * ddzu(k+1)                              &
    31913230                ) * 0.5_wp
     
    32603299       
    32613300       INTEGER(iwp) ::  i      !< grid index along x-direction
    3262        INTEGER(iwp) ::  ibit0  !< flag indicating 1st-order scheme along x-direction
    3263        INTEGER(iwp) ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
    3264        INTEGER(iwp) ::  ibit2  !< flag indicating 5th-order scheme along x-direction
    3265        INTEGER(iwp) ::  ibit3  !< flag indicating 1st-order scheme along y-direction
    3266        INTEGER(iwp) ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
    3267        INTEGER(iwp) ::  ibit5  !< flag indicating 5th-order scheme along y-direction
    3268        INTEGER(iwp) ::  ibit6  !< flag indicating 1st-order scheme along z-direction
    3269        INTEGER(iwp) ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
    3270        INTEGER(iwp) ::  ibit8  !< flag indicating 5th-order scheme along z-direction
    32713301       INTEGER(iwp) ::  j      !< grid index along y-direction
    32723302       INTEGER(iwp) ::  k      !< grid index along z-direction
     
    32763306       INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
    32773307       
    3278 #if defined( __nopointer )
    3279        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< advected scalar
    3280 #else
    3281        REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !< advected scalar
    3282 #endif
    3283 
     3308!       
     3309!--    sk is an array from parameter list. It should not be a pointer, because
     3310!--    in that case the compiler can not assume a stride 1 and cannot perform
     3311!--    a strided one vector load. Adding the CONTIGUOUS keyword makes things
     3312!--    even worse, because the compiler cannot assume strided one in the
     3313!--    caller side.
     3314       REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
     3315
     3316       REAL(wp) ::  ibit0  !< flag indicating 1st-order scheme along x-direction
     3317       REAL(wp) ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
     3318       REAL(wp) ::  ibit2  !< flag indicating 5th-order scheme along x-direction
     3319       REAL(wp) ::  ibit3  !< flag indicating 1st-order scheme along y-direction
     3320       REAL(wp) ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
     3321       REAL(wp) ::  ibit5  !< flag indicating 5th-order scheme along y-direction
     3322       REAL(wp) ::  ibit6  !< flag indicating 1st-order scheme along z-direction
     3323       REAL(wp) ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
     3324       REAL(wp) ::  ibit8  !< flag indicating 5th-order scheme along z-direction
    32843325       REAL(wp) ::  diss_d !< artificial dissipation term at grid box bottom
    32853326       REAL(wp) ::  div    !< diverence on scalar grid
     
    33093350          DO  k = nzb+1, nzb_max
    33103351
    3311              ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
    3312              ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
    3313              ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
     3352             ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
     3353             ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
     3354             ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
    33143355
    33153356             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    33703411          DO  k = nzb+1, nzb_max
    33713412
    3372              ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
    3373              ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
    3374              ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
     3413             ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
     3414             ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
     3415             ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
    33753416
    33763417             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     
    34333474             DO  k = nzb+1, nzb_max
    34343475
    3435                 ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
    3436                 ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
    3437                 ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
     3476                ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
     3477                ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
     3478                ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
    34383479
    34393480                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     
    34683509                                             )
    34693510
    3470                 ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
    3471                 ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
    3472                 ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
     3511                ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
     3512                ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
     3513                ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
    34733514
    34743515                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
     
    35053546!--             k index has to be modified near bottom and top, else array
    35063547!--             subscripts will be exceeded.
    3507                 ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
    3508                 ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
    3509                 ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
     3548                ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
     3549                ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
     3550                ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
    35103551
    35113552                k_ppp = k + 3 * ibit8
     
    35473588!--             by a not sufficient reduction of divergences near topography.
    35483589                div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
    3549                           - u(k,j,i)   * ( IBITS(advc_flags_1(k,j,i-1),0,1)    &
    3550                                          + IBITS(advc_flags_1(k,j,i-1),1,1)    &
    3551                                          + IBITS(advc_flags_1(k,j,i-1),2,1)    &
     3590                          - u(k,j,i)   * (                                     &
     3591                        REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )    &
     3592                      + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )    &
     3593                      + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )    &
    35523594                                         )                                     &
    35533595                          ) * ddx                                              &
    35543596                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    3555                           - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
    3556                                          + IBITS(advc_flags_1(k,j-1,i),4,1)    &
    3557                                          + IBITS(advc_flags_1(k,j-1,i),5,1)    &
     3597                          - v(k,j,i)   * (                                     &
     3598                        REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )    &
     3599                      + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )    &
     3600                      + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )    &
    35583601                                         )                                     &
    35593602                          ) * ddy                                              &
     
    35613604                                         ( ibit6 + ibit7 + ibit8 )             &
    35623605                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
    3563                                          ( IBITS(advc_flags_1(k-1,j,i),6,1)    &
    3564                                          + IBITS(advc_flags_1(k-1,j,i),7,1)    &
    3565                                          + IBITS(advc_flags_1(k-1,j,i),8,1)    &
     3606                                         (                                     &
     3607                        REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp )    &
     3608                      + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp )    &
     3609                      + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp )    &
    35663610                                         )                                     &     
    35673611                          ) * drho_air(k) * ddzw(k)
     
    36113655!--             k index has to be modified near bottom and top, else array
    36123656!--             subscripts will be exceeded.
    3613                 ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
    3614                 ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
    3615                 ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
     3657                ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
     3658                ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
     3659                ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
    36163660
    36173661                k_ppp = k + 3 * ibit8
     
    38173861
    38183862       INTEGER(iwp) ::  i      !< grid index along x-direction
    3819        INTEGER(iwp) ::  ibit9  !< flag indicating 1st-order scheme along x-direction
    3820        INTEGER(iwp) ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
    3821        INTEGER(iwp) ::  ibit11 !< flag indicating 5th-order scheme along x-direction
    3822        INTEGER(iwp) ::  ibit12 !< flag indicating 1st-order scheme along y-direction
    3823        INTEGER(iwp) ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
    3824        INTEGER(iwp) ::  ibit14 !< flag indicating 5th-order scheme along y-direction
    3825        INTEGER(iwp) ::  ibit15 !< flag indicating 1st-order scheme along z-direction
    3826        INTEGER(iwp) ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
    3827        INTEGER(iwp) ::  ibit17 !< flag indicating 5th-order scheme along z-direction
    38283863       INTEGER(iwp) ::  j      !< grid index along y-direction
    38293864       INTEGER(iwp) ::  k      !< grid index along z-direction
     
    38333868       INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
    38343869       
     3870       REAL(wp)    ::  ibit9  !< flag indicating 1st-order scheme along x-direction
     3871       REAL(wp)    ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
     3872       REAL(wp)    ::  ibit11 !< flag indicating 5th-order scheme along x-direction
     3873       REAL(wp)    ::  ibit12 !< flag indicating 1st-order scheme along y-direction
     3874       REAL(wp)    ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
     3875       REAL(wp)    ::  ibit14 !< flag indicating 5th-order scheme along y-direction
     3876       REAL(wp)    ::  ibit15 !< flag indicating 1st-order scheme along z-direction
     3877       REAL(wp)    ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
     3878       REAL(wp)    ::  ibit17 !< flag indicating 5th-order scheme along z-direction
    38353879       REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
    38363880       REAL(wp)    ::  div    !< diverence on u-grid
     
    38643908          DO  k = nzb+1, nzb_max
    38653909
    3866              ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
    3867              ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
    3868              ibit9  = IBITS(advc_flags_1(k,j,i-1),9,1)
     3910             ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
     3911             ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
     3912             ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
    38693913
    38703914             u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
    38713915             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
    3872                                        ( 37.0_wp * ibit11 * adv_mom_5             &
    3873                                     +     7.0_wp * ibit10 * adv_mom_3             &
    3874                                     +              ibit9  * adv_mom_1             &
     3916                                       ( 37.0_wp * ibit11 * adv_mom_5          &
     3917                                    +     7.0_wp * ibit10 * adv_mom_3          &
     3918                                    +              ibit9  * adv_mom_1          &
    38753919                                       ) *                                     &
    38763920                                     ( u(k,j,i)   + u(k,j,i-1) )               &
    3877                                 -      (  8.0_wp * ibit11 * adv_mom_5             &
    3878                                     +              ibit10 * adv_mom_3             &
     3921                                -      (  8.0_wp * ibit11 * adv_mom_5          &
     3922                                    +              ibit10 * adv_mom_3          &
    38793923                                       ) *                                     &
    38803924                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
    3881                                 +      (           ibit11 * adv_mom_5             &
     3925                                +      (           ibit11 * adv_mom_5          &
    38823926                                       ) *                                     &
    38833927                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
     
    38853929
    38863930              swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * (                &
    3887                                        ( 10.0_wp * ibit11 * adv_mom_5             &
    3888                                     +     3.0_wp * ibit10 * adv_mom_3             &
    3889                                     +              ibit9  * adv_mom_1             &
     3931                                       ( 10.0_wp * ibit11 * adv_mom_5          &
     3932                                    +     3.0_wp * ibit10 * adv_mom_3          &
     3933                                    +              ibit9  * adv_mom_1          &
    38903934                                       ) *                                     &
    38913935                                     ( u(k,j,i)   - u(k,j,i-1) )               &
    3892                                 -      (  5.0_wp * ibit11 * adv_mom_5             &
    3893                                     +              ibit10 * adv_mom_3             &
     3936                                -      (  5.0_wp * ibit11 * adv_mom_5          &
     3937                                    +              ibit10 * adv_mom_3          &
    38943938                                       ) *                                     &
    38953939                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
    3896                                 +      (           ibit11 * adv_mom_5             &
     3940                                +      (           ibit11 * adv_mom_5          &
    38973941                                       ) *                                     &
    38983942                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
     
    39053949             u_comp(k)         = u(k,j,i) + u(k,j,i-1) - gu
    39063950             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
    3907                              37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )                &
    3908                            -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )                &
     3951                             37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )             &
     3952                           -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )             &
    39093953                           +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
    39103954             swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                   &
    3911                              10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )                &
    3912                            -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )                &
     3955                             10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )             &
     3956                           -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )             &
    39133957                           +           ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
    39143958
     
    39223966          DO  k = nzb+1, nzb_max
    39233967
    3924              ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
    3925              ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
    3926              ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
     3968             ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
     3969             ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )
     3970             ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )
    39273971
    39283972             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
    39293973             swap_flux_y_local_u(k) = v_comp * (                              &
    3930                                    ( 37.0_wp * ibit14 * adv_mom_5                &
    3931                                 +     7.0_wp * ibit13 * adv_mom_3                &
    3932                                 +              ibit12 * adv_mom_1                &
     3974                                   ( 37.0_wp * ibit14 * adv_mom_5             &
     3975                                +     7.0_wp * ibit13 * adv_mom_3             &
     3976                                +              ibit12 * adv_mom_1             &
    39333977                                   ) *                                        &
    39343978                                     ( u(k,j,i)   + u(k,j-1,i) )              &
    3935                             -      (  8.0_wp * ibit14 * adv_mom_5                &
    3936                             +                  ibit13 * adv_mom_3                    &
     3979                            -      (  8.0_wp * ibit14 * adv_mom_5             &
     3980                            +                  ibit13 * adv_mom_3             &
    39373981                                   ) *                                        &
    39383982                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
    3939                         +      (               ibit14 * adv_mom_5                    &
     3983                        +      (               ibit14 * adv_mom_5             &
    39403984                               ) *                                            &
    39413985                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
     
    39433987
    39443988             swap_diss_y_local_u(k) = - ABS ( v_comp ) * (                    &
    3945                                    ( 10.0_wp * ibit14 * adv_mom_5                &
    3946                                 +     3.0_wp * ibit13 * adv_mom_3                &
    3947                                 +              ibit12 * adv_mom_1                &
     3989                                   ( 10.0_wp * ibit14 * adv_mom_5             &
     3990                                +     3.0_wp * ibit13 * adv_mom_3             &
     3991                                +              ibit12 * adv_mom_1             &
    39483992                                   ) *                                        &
    39493993                                     ( u(k,j,i)   - u(k,j-1,i) )              &
    3950                             -      (  5.0_wp * ibit14 * adv_mom_5                &
    3951                                 +              ibit13 * adv_mom_3                &
     3994                            -      (  5.0_wp * ibit14 * adv_mom_5             &
     3995                                +              ibit13 * adv_mom_3             &
    39523996                                   ) *                                        &
    39533997                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
    3954                             +      (           ibit14 * adv_mom_5                &
     3998                            +      (           ibit14 * adv_mom_5             &
    39553999                                   ) *                                        &
    39564000                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
     
    39834027             DO  k = nzb+1, nzb_max
    39844028
    3985                 ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
    3986                 ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
    3987                 ibit9  = IBITS(advc_flags_1(k,j,i),9,1)
     4029                ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
     4030                ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
     4031                ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1),  KIND = wp )
    39884032
    39894033                u_comp(k) = u(k,j,i+1) + u(k,j,i)
    39904034                flux_r(k) = ( u_comp(k) - gu ) * (                           &
    3991                           ( 37.0_wp * ibit11 * adv_mom_5                        &
    3992                        +     7.0_wp * ibit10 * adv_mom_3                        &
    3993                        +              ibit9  * adv_mom_1                        &
     4035                          ( 37.0_wp * ibit11 * adv_mom_5                     &
     4036                       +     7.0_wp * ibit10 * adv_mom_3                     &
     4037                       +              ibit9  * adv_mom_1                     &
    39944038                          ) *                                                &
    39954039                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
    3996                    -      (  8.0_wp * ibit11 * adv_mom_5                        &
    3997                        +              ibit10 * adv_mom_3                        &
     4040                   -      (  8.0_wp * ibit11 * adv_mom_5                     &
     4041                       +              ibit10 * adv_mom_3                     &
    39984042                          ) *                                                &
    39994043                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
    4000                    +      (           ibit11 * adv_mom_5                        &
     4044                   +      (           ibit11 * adv_mom_5                     &
    40014045                          ) *                                                &
    40024046                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
     
    40044048
    40054049                diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
    4006                           ( 10.0_wp * ibit11 * adv_mom_5                        &
    4007                        +     3.0_wp * ibit10 * adv_mom_3                        &
    4008                        +              ibit9  * adv_mom_1                        &
     4050                          ( 10.0_wp * ibit11 * adv_mom_5                     &
     4051                       +     3.0_wp * ibit10 * adv_mom_3                     &
     4052                       +              ibit9  * adv_mom_1                     &
    40094053                          ) *                                                &
    40104054                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
    4011                    -      (  5.0_wp * ibit11 * adv_mom_5                        &
    4012                        +              ibit10 * adv_mom_3                        &
     4055                   -      (  5.0_wp * ibit11 * adv_mom_5                     &
     4056                       +              ibit10 * adv_mom_3                     &
    40134057                          ) *                                                &
    40144058                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
    4015                    +      (           ibit11 * adv_mom_5                        &
     4059                   +      (           ibit11 * adv_mom_5                     &
    40164060                          ) *                                                &
    40174061                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
    40184062                                                     )
    40194063
    4020                 ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
    4021                 ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
    4022                 ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
     4064                ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
     4065                ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
     4066                ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
    40234067
    40244068                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    40254069                flux_n(k) = v_comp * (                                       &
    4026                           ( 37.0_wp * ibit14 * adv_mom_5                        &
    4027                        +     7.0_wp * ibit13 * adv_mom_3                        &
    4028                        +              ibit12 * adv_mom_1                        &
     4070                          ( 37.0_wp * ibit14 * adv_mom_5                     &
     4071                       +     7.0_wp * ibit13 * adv_mom_3                     &
     4072                       +              ibit12 * adv_mom_1                     &
    40294073                          ) *                                                &
    40304074                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
    4031                    -      (  8.0_wp * ibit14 * adv_mom_5                        &
    4032                        +              ibit13 * adv_mom_3                        &
     4075                   -      (  8.0_wp * ibit14 * adv_mom_5                     &
     4076                       +              ibit13 * adv_mom_3                     &
    40334077                          ) *                                                &
    40344078                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
    4035                    +      (           ibit14 * adv_mom_5                        &
     4079                   +      (           ibit14 * adv_mom_5                     &
    40364080                          ) *                                                &
    40374081                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
     
    40394083
    40404084                diss_n(k) = - ABS ( v_comp ) * (                             &
    4041                           ( 10.0_wp * ibit14 * adv_mom_5                        &
    4042                        +     3.0_wp * ibit13 * adv_mom_3                        &
    4043                        +              ibit12 * adv_mom_1                        &
     4085                          ( 10.0_wp * ibit14 * adv_mom_5                     &
     4086                       +     3.0_wp * ibit13 * adv_mom_3                     &
     4087                       +              ibit12 * adv_mom_1                     &
    40444088                          ) *                                                &
    40454089                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
    4046                    -      (  5.0_wp * ibit14 * adv_mom_5                        &
    4047                        +              ibit13 * adv_mom_3                        &
     4090                   -      (  5.0_wp * ibit14 * adv_mom_5                     &
     4091                       +              ibit13 * adv_mom_3                     &
    40484092                          ) *                                                &
    40494093                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
    4050                    +      (           ibit14 * adv_mom_5                        &
     4094                   +      (           ibit14 * adv_mom_5                     &
    40514095                          ) *                                                &
    40524096                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
     
    40554099!--             k index has to be modified near bottom and top, else array
    40564100!--             subscripts will be exceeded.
    4057                 ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
    4058                 ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
    4059                 ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
     4101                ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
     4102                ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
     4103                ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
    40604104
    40614105                k_ppp = k + 3 * ibit17
     
    40654109                w_comp    = w(k,j,i) + w(k,j,i-1)
    40664110                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
    4067                           ( 37.0_wp * ibit17 * adv_mom_5                        &
    4068                        +     7.0_wp * ibit16 * adv_mom_3                        &
    4069                        +              ibit15 * adv_mom_1                        &
     4111                          ( 37.0_wp * ibit17 * adv_mom_5                     &
     4112                       +     7.0_wp * ibit16 * adv_mom_3                     &
     4113                       +              ibit15 * adv_mom_1                     &
    40704114                          ) *                                                &
    40714115                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
    4072                    -      (  8.0_wp * ibit17 * adv_mom_5                        &
    4073                        +              ibit16 * adv_mom_3                        &
     4116                   -      (  8.0_wp * ibit17 * adv_mom_5                     &
     4117                       +              ibit16 * adv_mom_3                     &
    40744118                          ) *                                                &
    40754119                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
    4076                    +      (           ibit17 * adv_mom_5                        &
     4120                   +      (           ibit17 * adv_mom_5                     &
    40774121                          ) *                                                &
    40784122                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     
    40804124
    40814125                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
    4082                           ( 10.0_wp * ibit17 * adv_mom_5                        &
    4083                        +     3.0_wp * ibit16 * adv_mom_3                        &
    4084                        +              ibit15 * adv_mom_1                        &
     4126                          ( 10.0_wp * ibit17 * adv_mom_5                     &
     4127                       +     3.0_wp * ibit16 * adv_mom_3                     &
     4128                       +              ibit15 * adv_mom_1                     &
    40854129                          ) *                                                &
    40864130                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
    4087                    -      (  5.0_wp * ibit17 * adv_mom_5                        &
    4088                        +              ibit16 * adv_mom_3                        &
     4131                   -      (  5.0_wp * ibit17 * adv_mom_5                     &
     4132                       +              ibit16 * adv_mom_3                     &
    40894133                          ) *                                                &
    40904134                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
    4091                    +      (           ibit17 * adv_mom_5                        &
     4135                   +      (           ibit17 * adv_mom_5                     &
    40924136                           ) *                                               &
    40934137                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     
    40994143                div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 )             &
    41004144                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
    4101                                     * ( IBITS(advc_flags_1(k,j,i-1),9,1)      &
    4102                                       + IBITS(advc_flags_1(k,j,i-1),10,1)     &
    4103                                       + IBITS(advc_flags_1(k,j,i-1),11,1)     &
     4145                                    * (                                       &
     4146                    REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )      &
     4147                  + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )      &
     4148                  + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )      &
    41044149                                      )                                       &
    41054150                  ) * ddx                                                     &
    41064151               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    41074152                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
    4108                                     * ( IBITS(advc_flags_1(k,j-1,i),12,1)     &
    4109                                       + IBITS(advc_flags_1(k,j-1,i),13,1)     &
    4110                                       + IBITS(advc_flags_1(k,j-1,i),14,1)     &
     4153                                    * (                                       &
     4154                     REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )     &
     4155                   + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )     &
     4156                   + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )     &
    41114157                                      )                                       &
    41124158                  ) * ddy                                                     &
    41134159               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
    41144160                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
    4115                                     * ( IBITS(advc_flags_1(k-1,j,i),15,1)     &
    4116                                       + IBITS(advc_flags_1(k-1,j,i),16,1)     &
    4117                                       + IBITS(advc_flags_1(k-1,j,i),17,1)     &
     4161                                    * (                                       &
     4162                     REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp )     &
     4163                   + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp )     &
     4164                   + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp )     &
    41184165                                      )                                       & 
    41194166                  ) * drho_air(k) * ddzw(k)                                   &
     
    41864233!--             k index has to be modified near bottom and top, else array
    41874234!--             subscripts will be exceeded.
    4188                 ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
    4189                 ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
    4190                 ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
     4235                ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
     4236                ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
     4237                ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
    41914238
    41924239                k_ppp = k + 3 * ibit17
     
    43094356
    43104357       INTEGER(iwp) ::  i      !< grid index along x-direction
    4311        INTEGER(iwp) ::  ibit18 !< flag indicating 1st-order scheme along x-direction
    4312        INTEGER(iwp) ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
    4313        INTEGER(iwp) ::  ibit20 !< flag indicating 5th-order scheme along x-direction
    4314        INTEGER(iwp) ::  ibit21 !< flag indicating 1st-order scheme along y-direction
    4315        INTEGER(iwp) ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
    4316        INTEGER(iwp) ::  ibit23 !< flag indicating 5th-order scheme along y-direction
    4317        INTEGER(iwp) ::  ibit24 !< flag indicating 1st-order scheme along z-direction
    4318        INTEGER(iwp) ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
    4319        INTEGER(iwp) ::  ibit26 !< flag indicating 5th-order scheme along z-direction
    43204358       INTEGER(iwp) ::  j      !< grid index along y-direction
    43214359       INTEGER(iwp) ::  k      !< grid index along z-direction
     
    43254363       INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
    43264364       
     4365       REAL(wp)    ::  ibit18 !< flag indicating 1st-order scheme along x-direction
     4366       REAL(wp)    ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
     4367       REAL(wp)    ::  ibit20 !< flag indicating 5th-order scheme along x-direction
     4368       REAL(wp)    ::  ibit21 !< flag indicating 1st-order scheme along y-direction
     4369       REAL(wp)    ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
     4370       REAL(wp)    ::  ibit23 !< flag indicating 5th-order scheme along y-direction
     4371       REAL(wp)    ::  ibit24 !< flag indicating 1st-order scheme along z-direction
     4372       REAL(wp)    ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
     4373       REAL(wp)    ::  ibit26 !< flag indicating 5th-order scheme along z-direction
    43274374       REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
    43284375       REAL(wp)    ::  div    !< diverence on v-grid
     
    43554402          DO  k = nzb+1, nzb_max
    43564403
    4357              ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
    4358              ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
    4359              ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
     4404             ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
     4405             ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )
     4406             ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )
    43604407
    43614408             u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
     
    44134460          DO  k = nzb+1, nzb_max
    44144461
    4415              ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
    4416              ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
    4417              ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
     4462             ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
     4463             ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )
     4464             ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
    44184465
    44194466             v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
     
    44734520             DO  k = nzb+1, nzb_max
    44744521
    4475                 ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
    4476                 ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
    4477                 ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
     4522                ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
     4523                ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
     4524                ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
    44784525
    44794526                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
     
    45084555                                              )
    45094556
    4510                 ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
    4511                 ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
    4512                 ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
     4557                ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
     4558                ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
     4559                ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
    45134560
    45144561                v_comp(k) = v(k,j+1,i) + v(k,j,i)
     
    45454592!--             k index has to be modified near bottom and top, else array
    45464593!--             subscripts will be exceeded.
    4547                 ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
    4548                 ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
    4549                 ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
     4594                ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
     4595                ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
     4596                ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
    45504597
    45514598                k_ppp = k + 3 * ibit26
     
    45904637                                       * ( ibit18 + ibit19 + ibit20 )         &
    45914638                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
    4592                                        * ( IBITS(advc_flags_1(k,j,i-1),18,1)  &
    4593                                          + IBITS(advc_flags_1(k,j,i-1),19,1)  &
    4594                                          + IBITS(advc_flags_1(k,j,i-1),20,1)  &
     4639                                       * (                                    &
     4640                        REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )  &
     4641                      + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )  &
     4642                      + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )  &
    45954643                                         )                                    &
    45964644                  ) * ddx                                                     &
     
    45984646                                       * ( ibit21 + ibit22 + ibit23 )         &
    45994647                - ( v(k,j,i)     + v(k,j-1,i) )                               &
    4600                                        * ( IBITS(advc_flags_1(k,j-1,i),21,1)  &
    4601                                          + IBITS(advc_flags_1(k,j-1,i),22,1)  &
    4602                                          + IBITS(advc_flags_1(k,j-1,i),23,1)  &
     4648                                       * (                                    &
     4649                        REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )  &
     4650                      + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )  &
     4651                      + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )  &
    46034652                                         )                                    &
    46044653                  ) * ddy                                                     &
     
    46064655                                       * ( ibit24 + ibit25 + ibit26 )         &
    46074656                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
    4608                                        * ( IBITS(advc_flags_1(k-1,j,i),24,1)  &
    4609                                          + IBITS(advc_flags_1(k-1,j,i),25,1)  &
    4610                                          + IBITS(advc_flags_1(k-1,j,i),26,1)  &
     4657                                       * (                                    &
     4658                        REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp )  &
     4659                      + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp )  &
     4660                      + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp )  &
    46114661                                         )                                    &
    46124662                   ) * drho_air(k) * ddzw(k)                                  &
     
    46844734!--             k index has to be modified near bottom and top, else array
    46854735!--             subscripts will be exceeded.
    4686                 ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
    4687                 ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
    4688                 ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
     4736                ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
     4737                ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
     4738                ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
    46894739
    46904740                k_ppp = k + 3 * ibit26
     
    48114861
    48124862       INTEGER(iwp) ::  i      !< grid index along x-direction
    4813        INTEGER(iwp) ::  ibit27 !< flag indicating 1st-order scheme along x-direction
    4814        INTEGER(iwp) ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
    4815        INTEGER(iwp) ::  ibit29 !< flag indicating 5th-order scheme along x-direction
    4816        INTEGER(iwp) ::  ibit30 !< flag indicating 1st-order scheme along y-direction
    4817        INTEGER(iwp) ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
    4818        INTEGER(iwp) ::  ibit32 !< flag indicating 5th-order scheme along y-direction
    4819        INTEGER(iwp) ::  ibit33 !< flag indicating 1st-order scheme along z-direction
    4820        INTEGER(iwp) ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
    4821        INTEGER(iwp) ::  ibit35 !< flag indicating 5th-order scheme along z-direction
    48224863       INTEGER(iwp) ::  j      !< grid index along y-direction
    48234864       INTEGER(iwp) ::  k      !< grid index along z-direction
     
    48274868       INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
    48284869       
     4870       REAL(wp)    ::  ibit27 !< flag indicating 1st-order scheme along x-direction
     4871       REAL(wp)    ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
     4872       REAL(wp)    ::  ibit29 !< flag indicating 5th-order scheme along x-direction
     4873       REAL(wp)    ::  ibit30 !< flag indicating 1st-order scheme along y-direction
     4874       REAL(wp)    ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
     4875       REAL(wp)    ::  ibit32 !< flag indicating 5th-order scheme along y-direction
     4876       REAL(wp)    ::  ibit33 !< flag indicating 1st-order scheme along z-direction
     4877       REAL(wp)    ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
     4878       REAL(wp)    ::  ibit35 !< flag indicating 5th-order scheme along z-direction
    48294879       REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
    48304880       REAL(wp)    ::  div    !< divergence on w-grid
     
    48574907          DO  k = nzb+1, nzb_max
    48584908
    4859              ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
    4860              ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
    4861              ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
     4909             ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
     4910             ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )
     4911             ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )
    48624912
    48634913             u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
     
    49154965          DO  k = nzb+1, nzb_max
    49164966
    4917              ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
    4918              ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
    4919              ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
     4967             ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
     4968             ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
     4969             ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )
    49204970
    49214971             v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
     
    49815031             DO  k = nzb+1, nzb_max
    49825032
    4983                 ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
    4984                 ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
    4985                 ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
     5033                ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
     5034                ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
     5035                ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
    49865036
    49875037                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
     
    50165066                                              )
    50175067
    5018                 ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
    5019                 ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
    5020                 ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
     5068                ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),  KIND = wp )
     5069                ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
     5070                ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
    50215071
    50225072                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
     
    50535103!--             k index has to be modified near bottom and top, else array
    50545104!--             subscripts will be exceeded.
    5055                 ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
    5056                 ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
    5057                 ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
     5105                ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
     5106                ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
     5107                ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
    50585108
    50595109                k_ppp = k + 3 * ibit35
     
    50975147                div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
    50985148                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
    5099                                     * ( IBITS(advc_flags_1(k,j,i-1),27,1)     &
    5100                                       + IBITS(advc_flags_1(k,j,i-1),28,1)     &
    5101                                       + IBITS(advc_flags_1(k,j,i-1),29,1)     &
     5149                                    * (                                       &
     5150                     REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )     &
     5151                   + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )     &
     5152                   + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )     &
    51025153                                      )                                       &
    51035154                  ) * ddx                                                     &
    51045155              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    51055156                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
    5106                                     * ( IBITS(advc_flags_1(k,j-1,i),30,1)     &
    5107                                       + IBITS(advc_flags_1(k,j-1,i),31,1)     &
    5108                                       + IBITS(advc_flags_2(k,j-1,i),0,1)      &
     5157                                    * (                                       &
     5158                     REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )     &
     5159                   + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )     &
     5160                   + REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )     &
    51095161                                      )                                       &
    51105162                  ) * ddy                                                     &
    51115163              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
    51125164                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
    5113                                     * ( IBITS(advc_flags_2(k-1,j,i),1,1)      &
    5114                                       + IBITS(advc_flags_2(k-1,j,i),2,1)      &
    5115                                       + IBITS(advc_flags_2(k-1,j,i),3,1)      &
     5165                                    * (                                       &
     5166                     REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp )      &
     5167                   + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp )      &
     5168                   + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp )      &
    51165169                                      )                                       &
    51175170                  ) * drho_air_zw(k) * ddzu(k+1)                              &
     
    51765229!--             k index has to be modified near bottom and top, else array
    51775230!--             subscripts will be exceeded.
    5178                 ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
    5179                 ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
    5180                 ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
     5231                ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
     5232                ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
     5233                ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
    51815234
    51825235                k_ppp = k + 3 * ibit35
Note: See TracChangeset for help on using the changeset viewer.