Home


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

Changeset 1353


Ignore:
Timestamp:
Apr 8, 2014 3:21:23 PM (11 years ago)
Author:
heinze
Message:

REAL constants provided with KIND-attribute

Location:
palm/trunk/SOURCE
Files:
76 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/advec_s_bc.f90

    r1347 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    170170!-- Array sk_p requires 2 extra elements for each dimension
    171171    ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
    172     sk_p = 0.0
     172    sk_p = 0.0_wp
    173173
    174174!
    175175!-- Assign reciprocal values in order to avoid divisions later
    176     f2    = 0.5
    177     f4    = 0.25
    178     f8    = 0.125
    179     f12   = 0.8333333333333333E-01
    180     f24   = 0.4166666666666666E-01
    181     f48   = 0.2083333333333333E-01
    182     f1920 = 0.5208333333333333E-03
     176    f2    = 0.5_wp
     177    f4    = 0.25_wp
     178    f8    = 0.125_wp
     179    f12   = 0.8333333333333333E-01_wp
     180    f24   = 0.4166666666666666E-01_wp
     181    f48   = 0.2083333333333333E-01_wp
     182    f1920 = 0.5208333333333333E-03_wp
    183183
    184184!
    … …  
    236236!
    237237!-- Initialise control density
    238     d = 0.0
     238    d = 0.0_wp
    239239
    240240!
    241241!-- Determine maxima of the first and second derivative in x-direction
    242     fmax_l = 0.0
     242    fmax_l = 0.0_wp
    243243    DO  i = nxl, nxr
    244244       DO  j = nys, nyn
    245245          DO  k = nzb+1, nzt
    246              numera = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
     246             numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
    247247             denomi  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    248248             fmax_l(1) = MAX( fmax_l(1) , numera )
    … …  
    258258#endif
    259259
    260     fmax = 0.04 * fmax
     260    fmax = 0.04_wp * fmax
    261261
    262262!
    … …  
    271271              sw(nzb+1:nzt,nxl-1:nxr+1)                                        &
    272272            )
    273     imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
     273    imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    274274
    275275!
    … …  
    287287       DO  i = nxl-1, nxr+1
    288288          DO  k = nzb+1, nzt
    289              a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    290              a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i)              &
    291                                               + sk_p(k,j,i-1) )
    292              a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1)           &
    293                          + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1)        &
    294                          + 9.0 * sk_p(k,j,i-2) ) * f1920
    295              a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1)           &
    296                          - 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2)          &
     289             a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
     290             a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i)        &
     291                                                 + sk_p(k,j,i-1) )
     292             a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)  &
     293                         + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)  &
     294                         + 9.0_wp * sk_p(k,j,i-2) ) * f1920
     295             a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)   &
     296                         - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)    &
    297297                       ) * f48
    298              a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1)           &
    299                          - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1)           &
    300                          - 3.0 * sk_p(k,j,i-2) ) * f48
     298             a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)     &
     299                         - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)     &
     300                         - 3.0_wp * sk_p(k,j,i-2) ) * f48
    301301          ENDDO
    302302       ENDDO
    … …  
    309309             cip  =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    310310             cim  = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    311              cipf = 1.0 - 2.0 * cip
    312              cimf = 1.0 - 2.0 * cim
    313              ip   =   a0(k,i)   * f2  * ( 1.0 - cipf )                         &
    314                     + a1(k,i)   * f8  * ( 1.0 - cipf*cipf )                    &
    315                     + a2(k,i)   * f24 * ( 1.0 - cipf*cipf*cipf )
    316              im   =   a0(k,i+1) * f2  * ( 1.0 - cimf )                         &
    317                     - a1(k,i+1) * f8  * ( 1.0 - cimf*cimf )                    &
    318                     + a2(k,i+1) * f24 * ( 1.0 - cimf*cimf*cimf )
     311             cipf = 1.0_wp - 2.0_wp * cip
     312             cimf = 1.0_wp - 2.0_wp * cim
     313             ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                      &
     314                    + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     315                    + a2(k,i)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
     316             im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                      &
     317                    - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
     318                    + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    319319             ip   = MAX( ip, 0.0_wp )
    320320             im   = MAX( im, 0.0_wp )
    321              ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
    322              impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15) )
     321             ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     322             impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) )
    323323
    324324             cip  =  MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    325325             cim  = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    326              cipf = 1.0 - 2.0 * cip
    327              cimf = 1.0 - 2.0 * cim
    328              ip   =   a0(k,i-1) * f2  * ( 1.0 - cipf )                         &
    329                     + a1(k,i-1) * f8  * ( 1.0 - cipf*cipf )                    &
    330                     + a2(k,i-1) * f24 * ( 1.0 - cipf*cipf*cipf )
    331              im   =   a0(k,i)   * f2  * ( 1.0 - cimf )                         &
    332                     - a1(k,i)   * f8  * ( 1.0 - cimf*cimf )                    &
    333                     + a2(k,i)   * f24 * ( 1.0 - cimf*cimf*cimf )
     326             cipf = 1.0_wp - 2.0_wp * cip
     327             cimf = 1.0_wp - 2.0_wp * cim
     328             ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                      &
     329                    + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
     330                    + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
     331             im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                      &
     332                    - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     333                    + a2(k,i)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    334334             ip   = MAX( ip, 0.0_wp )
    335335             im   = MAX( im, 0.0_wp )
    336              ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15) )
    337              immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
     336             ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) )
     337             immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    338338          ENDDO
    339339       ENDDO
    … …  
    343343       DO  i = nxl-2, nxr+2
    344344          DO  k = nzb+1, nzt
    345              m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
     345             m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
    346346             m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    347              IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
     347             IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
    348348                m1(k,i) = m1z / m1n
    349                 IF ( m1(k,i) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0
     349                IF ( m1(k,i) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0_wp
    350350             ELSEIF ( m1n < m1z )  THEN
    351                 m1(k,i) = -1.0
     351                m1(k,i) = -1.0_wp
    352352             ELSE
    353                 m1(k,i) = 0.0
     353                m1(k,i) = 0.0_wp
    354354             ENDIF
    355355          ENDDO
    … …  
    358358!
    359359!--    Compute switch sw
    360        sw = 0.0
     360       sw = 0.0_wp
    361361       DO  i = nxl-1, nxr+1
    362362          DO  k = nzb+1, nzt
    363              m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) /                            &
     363             m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) /                         &
    364364                  MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
    365              IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0
    366 
    367              m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) /                            &
     365             IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0_wp
     366
     367             m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) /                         &
    368368                  MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
    369              IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0
    370 
    371              t1 = 0.35
    372              t2 = 0.35
    373              IF ( m1(k,i) == -1.0 )  t2 = 0.12
     369             IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0_wp
     370
     371             t1 = 0.35_wp
     372             t2 = 0.35_wp
     373             IF ( m1(k,i) == -1.0_wp )  t2 = 0.12_wp
    374374
    375375!--          *VOCL STMT,IF(10)
    376              IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0   &
    377                   .OR.  m2 > t2  .OR.  m3 > t2  .OR.                           &
    378                   ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0  .AND.              &
    379                     m1(k,i) /= -1.0  .AND.  m1(k,i+1) /= -1.0 )                &
    380                 )  sw(k,i) = 1.0
     376             IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp                   &
     377                  .OR. m1(k,i+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
     378                  ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0_wp  .AND.           &
     379                    m1(k,i) /= -1.0_wp  .AND.  m1(k,i+1) /= -1.0_wp )          &
     380                )  sw(k,i) = 1.0_wp
    381381          ENDDO
    382382       ENDDO
    … …  
    389389
    390390!--          *VOCL STMT,IF(10)
    391              IF ( sw(k,i) == 1.0 )  THEN
     391             IF ( sw(k,i) == 1.0_wp )  THEN
    392392                snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
    393                 IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
     393                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    394394                sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
    395395                sterm = MIN( sterm, 0.9999_wp )
    … …  
    402402                ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                    &
    403403                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
    404                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
     404                            eex(ix) -                                          &
     405                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    405406                                                                )              &
    406407                                                          )
    407                 IF ( sterm == 0.0001 )  ippe(k,i) = sk_p(k,j,i) * cip
    408                 IF ( sterm == 0.9999 )  ippe(k,i) = sk_p(k,j,i) * cip
     408                IF ( sterm == 0.0001_wp )  ippe(k,i) = sk_p(k,j,i) * cip
     409                IF ( sterm == 0.9999_wp )  ippe(k,i) = sk_p(k,j,i) * cip
    409410
    410411                snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
    411                 IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
     412                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    412413                sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
    413414                sterm = MIN( sterm, 0.9999_wp )
    … …  
    420421                imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                    &
    421422                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
    422                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
     423                            eex(ix) -                                          &
     424                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    423425                                                                )              &
    424426                                                          )
    425                 IF ( sterm == 0.0001 )  imme(k,i) = sk_p(k,j,i) * cim
    426                 IF ( sterm == 0.9999 )  imme(k,i) = sk_p(k,j,i) * cim
     427                IF ( sterm == 0.0001_wp )  imme(k,i) = sk_p(k,j,i) * cim
     428                IF ( sterm == 0.9999_wp )  imme(k,i) = sk_p(k,j,i) * cim
    427429             ENDIF
    428430
    429431!--          *VOCL STMT,IF(10)
    430              IF ( sw(k,i+1) == 1.0 )  THEN
     432             IF ( sw(k,i+1) == 1.0_wp )  THEN
    431433                snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
    432                 IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
     434                IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
    433435                sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
    434436                sterm = MIN( sterm, 0.9999_wp )
    … …  
    441443                impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                    &
    442444                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
    443                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
     445                            eex(ix) -                                          &
     446                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    444447                                                                )              &
    445448                                                          )
    … …  
    449452
    450453!--          *VOCL STMT,IF(10)
    451              IF ( sw(k,i-1) == 1.0 )  THEN
     454             IF ( sw(k,i-1) == 1.0_wp )  THEN
    452455                snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
    453                 IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
     456                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    454457                sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
    455458                sterm = MIN( sterm, 0.9999_wp )
    … …  
    462465                ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                    &
    463466                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
    464                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0_wp - 2.0 * cip ) ) &
     467                            eex(ix) -                                          &
     468                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    465469                                                                )              &
    466470                                                          )
    … …  
    538542!
    539543!-- Determine the maxima of the first and second derivative in y-direction
    540     fmax_l = 0.0
     544    fmax_l = 0.0_wp
    541545    DO  i = nxl, nxr
    542546       DO  j = nys, nyn
    543547          DO  k = nzb+1, nzt
    544              numera = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
     548             numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
    545549             denomi  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    546550             fmax_l(1) = MAX( fmax_l(1) , numera )
    … …  
    556560#endif
    557561
    558     fmax = 0.04 * fmax
     562    fmax = 0.04_wp * fmax
    559563
    560564!
    … …  
    569573              sw(nzb+1:nzt,nys-1:nyn+1)                                        &
    570574            )
    571     imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
     575    imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    572576
    573577!
    … …  
    579583       DO  j = nys-1, nyn+1
    580584          DO  k = nzb+1, nzt
    581              a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    582              a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i)              &
    583                                               + sk_p(k,j-1,i) )
    584              a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i)           &
    585                          + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i)        &
    586                          + 9.0 * sk_p(k,j-2,i) ) * f1920
    587              a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i)           &
    588                          - 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i)          &
     585             a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
     586             a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i)        &
     587                                                 + sk_p(k,j-1,i) )
     588             a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)  &
     589                         + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)  &
     590                         + 9.0_wp * sk_p(k,j-2,i) ) * f1920
     591             a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)   &
     592                         - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)   &
    589593                       ) * f48
    590              a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i)           &
    591                          - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i)           &
    592                          - 3.0 * sk_p(k,j-2,i) ) * f48
     594             a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)     &
     595                         - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)     &
     596                         - 3.0_wp * sk_p(k,j-2,i) ) * f48
    593597          ENDDO
    594598       ENDDO
    … …  
    601605             cip  =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    602606             cim  = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    603              cipf = 1.0 - 2.0 * cip
    604              cimf = 1.0 - 2.0 * cim
    605              ip   =   a0(k,j)   * f2  * ( 1.0 - cipf )                         &
    606                     + a1(k,j)   * f8  * ( 1.0 - cipf*cipf )                    &
    607                     + a2(k,j)   * f24 * ( 1.0 - cipf*cipf*cipf )
    608              im   =   a0(k,j+1) * f2  * ( 1.0 - cimf )                         &
    609                     - a1(k,j+1) * f8  * ( 1.0 - cimf*cimf )                    &
    610                     + a2(k,j+1) * f24 * ( 1.0 - cimf*cimf*cimf )
     607             cipf = 1.0_wp - 2.0_wp * cip
     608             cimf = 1.0_wp - 2.0_wp * cim
     609             ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
     610                    + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     611                    + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
     612             im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                      &
     613                    - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
     614                    + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    611615             ip   = MAX( ip, 0.0_wp )
    612616             im   = MAX( im, 0.0_wp )
    613              ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
    614              impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15) )
     617             ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     618             impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) )
    615619
    616620             cip  =  MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    617621             cim  = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    618              cipf = 1.0 - 2.0 * cip
    619              cimf = 1.0 - 2.0 * cim
    620              ip   =   a0(k,j-1) * f2  * ( 1.0 - cipf )                         &
    621                     + a1(k,j-1) * f8  * ( 1.0 - cipf*cipf )                    &
    622                     + a2(k,j-1) * f24 * ( 1.0 - cipf*cipf*cipf )
    623              im   =   a0(k,j)   * f2  * ( 1.0 - cimf )                         &
    624                     - a1(k,j)   * f8  * ( 1.0 - cimf*cimf )                    &
    625                     + a2(k,j)   * f24 * ( 1.0 - cimf*cimf*cimf )
     622             cipf = 1.0_wp - 2.0_wp * cip
     623             cimf = 1.0_wp - 2.0_wp * cim
     624             ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                      &
     625                    + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
     626                    + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
     627             im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
     628                    - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     629                    + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    626630             ip   = MAX( ip, 0.0_wp )
    627631             im   = MAX( im, 0.0_wp )
    628              ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15) )
    629              immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
     632             ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) )
     633             immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    630634          ENDDO
    631635       ENDDO
    … …  
    635639       DO  j = nys-2, nyn+2
    636640          DO  k = nzb+1, nzt
    637              m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
     641             m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
    638642             m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    639              IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
     643             IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
    640644                m1(k,j) = m1z / m1n
    641                 IF ( m1(k,j) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0
     645                IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
    642646             ELSEIF ( m1n < m1z )  THEN
    643                 m1(k,j) = -1.0
     647                m1(k,j) = -1.0_wp
    644648             ELSE
    645                 m1(k,j) = 0.0
     649                m1(k,j) = 0.0_wp
    646650             ENDIF
    647651          ENDDO
    … …  
    650654!
    651655!--    Compute switch sw
    652        sw = 0.0
     656       sw = 0.0_wp
    653657       DO  j = nys-1, nyn+1
    654658          DO  k = nzb+1, nzt
    655              m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) /                            &
     659             m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                            &
    656660                  MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
    657              IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0
    658 
    659              m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) /                            &
     661             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
     662
     663             m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                            &
    660664                  MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
    661              IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0
    662 
    663              t1 = 0.35
    664              t2 = 0.35
    665              IF ( m1(k,j) == -1.0 )  t2 = 0.12
     665             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
     666
     667             t1 = 0.35_wp
     668             t2 = 0.35_wp
     669             IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
    666670
    667671!--          *VOCL STMT,IF(10)
    668              IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0   &
    669                   .OR.  m2 > t2  .OR.  m3 > t2  .OR.                           &
    670                   ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0  .AND.              &
    671                     m1(k,j) /= -1.0  .AND.  m1(k,j+1) /= -1.0 )                &
    672                 )  sw(k,j) = 1.0
     672             IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
     673                  .OR. m1(k,j+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
     674                  ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0_wp  .AND.           &
     675                    m1(k,j) /= -1.0_wp  .AND.  m1(k,j+1) /= -1.0_wp )          &
     676                )  sw(k,j) = 1.0_wp
    673677          ENDDO
    674678       ENDDO
    … …  
    681685
    682686!--          *VOCL STMT,IF(10)
    683              IF ( sw(k,j) == 1.0 )  THEN
     687             IF ( sw(k,j) == 1.0_wp )  THEN
    684688                snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
    685                 IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
     689                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    686690                sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
    687691                sterm = MIN( sterm, 0.9999_wp )
    … …  
    694698                ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                    &
    695699                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
    696                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
     700                            eex(ix) -                                          &
     701                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    697702                                                                )              &
    698703                                                          )
    … …  
    701706
    702707                snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
    703                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
     708                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    704709                sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
    705710                sterm = MIN( sterm, 0.9999_wp )
    … …  
    712717                imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                    &
    713718                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
    714                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
     719                            eex(ix) -                                          &
     720                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    715721                                                                )              &
    716722                                                          )
    … …  
    720726
    721727!--          *VOCL STMT,IF(10)
    722              IF ( sw(k,j+1) == 1.0 )  THEN
     728             IF ( sw(k,j+1) == 1.0_wp )  THEN
    723729                snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
    724                 IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
     730                IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
    725731                sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
    726732                sterm = MIN( sterm, 0.9999_wp )
    … …  
    733739                impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                    &
    734740                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
    735                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
     741                            eex(ix) -                                          &
     742                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    736743                                                                )              &
    737744                                                          )
    … …  
    741748
    742749!--          *VOCL STMT,IF(10)
    743              IF ( sw(k,j-1) == 1.0 )  THEN
     750             IF ( sw(k,j-1) == 1.0_wp )  THEN
    744751                snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
    745                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
     752                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    746753                sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
    747754                sterm = MIN( sterm, 0.9999_wp )
    … …  
    754761                ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                    &
    755762                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
    756                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
     763                            eex(ix) -                                          &
     764                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    757765                                                                )              &
    758766                                                          )
    … …  
    769777       DO  j = nys, nyn
    770778          DO  k = nzb+1, nzt
    771              fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    772                     - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
    773              fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
    774                     - ( 1.0 - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
     779             fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
     780                    - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
     781             fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
     782                    - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
    775783             tendcy = fplus - fminus
    776784!
    … …  
    781789!--          Density correction because of possible remaining divergences
    782790             d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
    783              sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
    784                            ( 1.0 + d_new )
     791             sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
     792                           ( 1.0_wp + d_new )
    785793             d(k,j,i)  = d_new
    786794          ENDDO
    … …  
    800808!-- Initialise for the computation of heat fluxes (see below; required in
    801809!-- UP flow_statistics)
    802     IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0
     810    IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0_wp
    803811
    804812!
    … …  
    939947!
    940948!-- Determine the maxima of the first and second derivative in z-direction
    941     fmax_l = 0.0
     949    fmax_l = 0.0_wp
    942950    DO  i = nxl, nxr
    943951       DO  j = nys, nyn
    944952          DO  k = nzb, nzt+1
    945              numera = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
     953             numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
    946954             denomi  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
    947955             fmax_l(1) = MAX( fmax_l(1) , numera )
    … …  
    957965#endif
    958966
    959     fmax = 0.04 * fmax
     967    fmax = 0.04_wp * fmax
    960968
    961969!
    … …  
    970978              sw(nzb:nzt+1,nys:nyn)                                            &
    971979            )
    972     imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
     980    imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    973981
    974982!
    … …  
    980988       DO  j = nys, nyn
    981989          DO  k = nzb, nzt+1
    982              a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
    983              a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i)              &
    984                                               + sk_p(k-1,j,i) )
    985              a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i)           &
    986                          + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i)        &
    987                          + 9.0 * sk_p(k-2,j,i) ) * f1920
    988              a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i)           &
    989                          - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i)          &
     990             a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
     991             a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i)        &
     992                                                 + sk_p(k-1,j,i) )
     993             a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)  &
     994                         + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)  &
     995                         + 9.0_wp * sk_p(k-2,j,i) ) * f1920
     996             a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)   &
     997                         - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)   &
    990998                       ) * f48
    991              a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i)           &
    992                          - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i)           &
    993                          - 3.0 * sk_p(k-2,j,i) ) * f48
     999             a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)     &
     1000                         - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)     &
     1001                         - 3.0_wp * sk_p(k-2,j,i) ) * f48
    9941002          ENDDO
    9951003       ENDDO
    … …  
    10021010             cip  =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    10031011             cim  = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    1004              cipf = 1.0 - 2.0 * cip
    1005              cimf = 1.0 - 2.0 * cim
    1006              ip   =   a0(k,j)   * f2  * ( 1.0 - cipf )                         &
    1007                     + a1(k,j)   * f8  * ( 1.0 - cipf*cipf )                    &
    1008                     + a2(k,j)   * f24 * ( 1.0 - cipf*cipf*cipf )
    1009              im   =   a0(k+1,j) * f2  * ( 1.0 - cimf )                         &
    1010                     - a1(k+1,j) * f8  * ( 1.0 - cimf*cimf )                    &
    1011                     + a2(k+1,j) * f24 * ( 1.0 - cimf*cimf*cimf )
     1012             cipf = 1.0_wp - 2.0_wp * cip
     1013             cimf = 1.0_wp - 2.0_wp * cim
     1014             ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
     1015                    + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     1016                    + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
     1017             im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                      &
     1018                    - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                 &
     1019                    + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    10121020             ip   = MAX( ip, 0.0_wp )
    10131021             im   = MAX( im, 0.0_wp )
    1014              ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
    1015              impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15) )
     1022             ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
     1023             impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) )
    10161024
    10171025             cip  =  MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    10181026             cim  = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    1019              cipf = 1.0 - 2.0 * cip
    1020              cimf = 1.0 - 2.0 * cim
    1021              ip   =   a0(k-1,j) * f2  * ( 1.0 - cipf )                         &
    1022                     + a1(k-1,j) * f8  * ( 1.0 - cipf*cipf )                    &
    1023                     + a2(k-1,j) * f24 * ( 1.0 - cipf*cipf*cipf )
    1024              im   =   a0(k,j)   * f2  * ( 1.0 - cimf )                         &
    1025                     - a1(k,j)   * f8  * ( 1.0 - cimf*cimf )                    &
    1026                     + a2(k,j)   * f24 * ( 1.0 - cimf*cimf*cimf )
     1027             cipf = 1.0_wp - 2.0_wp * cip
     1028             cimf = 1.0_wp - 2.0_wp * cim
     1029             ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                      &
     1030                    + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                 &
     1031                    + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf )
     1032             im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
     1033                    - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     1034                    + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    10271035             ip   = MAX( ip, 0.0_wp )
    10281036             im   = MAX( im, 0.0_wp )
    1029              ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15) )
    1030              immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
     1037             ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) )
     1038             immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
    10311039          ENDDO
    10321040       ENDDO
    … …  
    10361044       DO  j = nys, nyn
    10371045          DO  k = nzb-1, nzt+2
    1038              m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
     1046             m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
    10391047             m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
    1040              IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
     1048             IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
    10411049                m1(k,j) = m1z / m1n
    1042                 IF ( m1(k,j) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0
     1050                IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
    10431051             ELSEIF ( m1n < m1z )  THEN
    1044                 m1(k,j) = -1.0
     1052                m1(k,j) = -1.0_wp
    10451053             ELSE
    1046                 m1(k,j) = 0.0
     1054                m1(k,j) = 0.0_wp
    10471055             ENDIF
    10481056          ENDDO
    … …  
    10511059!
    10521060!--    Compute switch sw
    1053        sw = 0.0
     1061       sw = 0.0_wp
    10541062       DO  j = nys, nyn
    10551063          DO  k = nzb, nzt+1
    1056              m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) /                            &
     1064             m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                         &
    10571065                  MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
    1058              IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0
    1059 
    1060              m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) /                            &
     1066             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
     1067
     1068             m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                         &
    10611069                  MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
    1062              IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0
    1063 
    1064              t1 = 0.35
    1065              t2 = 0.35
    1066              IF ( m1(k,j) == -1.0 )  t2 = 0.12
     1070             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
     1071
     1072             t1 = 0.35_wp
     1073             t2 = 0.35_wp
     1074             IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
    10671075
    10681076!--          *VOCL STMT,IF(10)
    1069              IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0   &
    1070                   .OR.  m2 > t2  .OR.  m3 > t2  .OR.                           &
    1071                   ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0  .AND.              &
    1072                     m1(k,j) /= -1.0  .AND.  m1(k+1,j) /= -1.0 )                &
    1073                 )  sw(k,j) = 1.0
     1077             IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
     1078                  .OR. m1(k+1,j) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
     1079                  ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0_wp  .AND.           &
     1080                    m1(k,j) /= -1.0_wp  .AND.  m1(k+1,j) /= -1.0_wp )          &
     1081                )  sw(k,j) = 1.0_wp
    10741082          ENDDO
    10751083       ENDDO
    … …  
    10821090
    10831091!--          *VOCL STMT,IF(10)
    1084              IF ( sw(k,j) == 1.0 )  THEN
     1092             IF ( sw(k,j) == 1.0_wp )  THEN
    10851093                snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
    1086                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
     1094                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    10871095                sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
    10881096                sterm = MIN( sterm, 0.9999_wp )
    … …  
    10951103                ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
    10961104                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
    1097                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
     1105                            eex(ix) -                                          &
     1106                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    10981107                                                                )              &
    10991108                                                          )
    1100                 IF ( sterm == 0.0001 )  ippe(k,j) = sk_p(k,j,i) * cip
    1101                 IF ( sterm == 0.9999 )  ippe(k,j) = sk_p(k,j,i) * cip
     1109                IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     1110                IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
    11021111
    11031112                snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
    1104                 IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
     1113                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    11051114                sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
    11061115                sterm = MIN( sterm, 0.9999_wp )
    … …  
    11131122                imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
    11141123                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
    1115                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
     1124                            eex(ix) -                                          &
     1125                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) &
    11161126                                                                )              &
    11171127                                                          )
    … …  
    11211131
    11221132!--          *VOCL STMT,IF(10)
    1123              IF ( sw(k+1,j) == 1.0 )  THEN
     1133             IF ( sw(k+1,j) == 1.0_wp )  THEN
    11241134                snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
    1125                 IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
     1135                IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
    11261136                sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
    11271137                sterm = MIN( sterm, 0.9999_wp )
    … …  
    11341144                impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
    11351145                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
    1136                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
     1146                            eex(ix) -                                          &
     1147                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    11371148                                                                )              &
    11381149                                                          )
    … …  
    11421153
    11431154!--          *VOCL STMT,IF(10)
    1144              IF ( sw(k-1,j) == 1.0 )  THEN
     1155             IF ( sw(k-1,j) == 1.0_wp )  THEN
    11451156                snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
    1146                 IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
     1157                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
    11471158                sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
    11481159                sterm = MIN( sterm, 0.9999_wp )
    … …  
    11551166                ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
    11561167                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
    1157                             eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
     1168                            eex(ix) -                                          &
     1169                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    11581170                                                                )              &
    11591171                                                          )
    … …  
    11701182       DO  j = nys, nyn
    11711183          DO  k = nzb+1, nzt
    1172              fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    1173                     - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
    1174              fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
    1175                     - ( 1.0 - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
     1184             fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
     1185                    - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
     1186             fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
     1187                    - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
    11761188             tendcy = fplus - fminus
    11771189!
    … …  
    11821194!--          Density correction because of possible remaining divergences
    11831195             d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
    1184              sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
    1185                            ( 1.0 + d_new )
     1196             sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
     1197                           ( 1.0_wp + d_new )
    11861198!
    11871199!--          Store heat flux for subsequent statistics output.

  • palm/trunk/SOURCE/advec_s_pw.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    102102             DO  k = nzb_s_inner(j,i)+1, nzt
    103103                tend(k,j,i) = tend(k,j,i)                                      &
    104               -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
    105                      - ( u(k,j,i)   - u_gtrans ) * ( sk(k,j,i-1) - sk(k,j,i) ) &
    106                      ) * ddx                                                   &
    107               -0.5 * ( ( v(k,j+1,i) - v_gtrans ) * ( sk(k,j+1,i) - sk(k,j,i) ) &
    108                      - ( v(k,j,i)   - v_gtrans ) * ( sk(k,j-1,i) - sk(k,j,i) ) &
    109                      ) * ddy                                                   &
    110               -      (   w(k,j,i)                * ( sk(k+1,j,i) - sk(k,j,i) ) &
    111                      -   w(k-1,j,i)              * ( sk(k-1,j,i) - sk(k,j,i) ) &
    112                      ) * dd2zu(k)
     104              -0.5_wp * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
     105                        - ( u(k,j,i)   - u_gtrans ) * ( sk(k,j,i-1) - sk(k,j,i) ) &
     106                        ) * ddx                                                   &
     107              -0.5_wp * ( ( v(k,j+1,i) - v_gtrans ) * ( sk(k,j+1,i) - sk(k,j,i) ) &
     108                        - ( v(k,j,i)   - v_gtrans ) * ( sk(k,j-1,i) - sk(k,j,i) ) &
     109                        ) * ddy                                                   &
     110              -         (   w(k,j,i)                * ( sk(k+1,j,i) - sk(k,j,i) ) &
     111                        -   w(k-1,j,i)              * ( sk(k-1,j,i) - sk(k,j,i) ) &
     112                        ) * dd2zu(k)
    113113             ENDDO
    114114          ENDDO
    … …  
    153153       DO  k = nzb_s_inner(j,i)+1, nzt
    154154          tend(k,j,i) = tend(k,j,i)                                            &
    155               -0.5 * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
    156                      - ( u(k,j,i)   - u_gtrans ) * ( sk(k,j,i-1) - sk(k,j,i) ) &
    157                      ) * ddx                                                   &
    158               -0.5 * ( ( v(k,j+1,i) - v_gtrans ) * ( sk(k,j+1,i) - sk(k,j,i) ) &
    159                      - ( v(k,j,i)   - v_gtrans ) * ( sk(k,j-1,i) - sk(k,j,i) ) &
    160                      ) * ddy                                                   &
    161               -      (   w(k,j,i)                * ( sk(k+1,j,i) - sk(k,j,i) ) &
    162                      -   w(k-1,j,i)              * ( sk(k-1,j,i) - sk(k,j,i) ) &
    163                      ) * dd2zu(k)
     155              -0.5_wp * ( ( u(k,j,i+1) - u_gtrans ) * ( sk(k,j,i+1) - sk(k,j,i) ) &
     156                        - ( u(k,j,i)   - u_gtrans ) * ( sk(k,j,i-1) - sk(k,j,i) ) &
     157                        ) * ddx                                                   &
     158              -0.5_wp * ( ( v(k,j+1,i) - v_gtrans ) * ( sk(k,j+1,i) - sk(k,j,i) ) &
     159                        - ( v(k,j,i)   - v_gtrans ) * ( sk(k,j-1,i) - sk(k,j,i) ) &
     160                        ) * ddy                                                   &
     161              -         (   w(k,j,i)                * ( sk(k+1,j,i) - sk(k,j,i) ) &
     162                        -   w(k-1,j,i)              * ( sk(k-1,j,i) - sk(k,j,i) ) &
     163                        ) * dd2zu(k)
    164164       ENDDO
    165165

  • palm/trunk/SOURCE/advec_s_up.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    106106!
    107107!--             x-direction
    108                 ukomp = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
    109                 IF ( ukomp > 0.0 )  THEN
     108                ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
     109                IF ( ukomp > 0.0_wp )  THEN
    110110                   tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    111111                                         ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
    … …  
    116116!
    117117!--             y-direction
    118                 vkomp = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
    119                 IF ( vkomp > 0.0 )  THEN
     118                vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
     119                IF ( vkomp > 0.0_wp )  THEN
    120120                   tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    121121                                         ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
    … …  
    126126!
    127127!--             z-direction
    128                 wkomp = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
    129                 IF ( wkomp > 0.0 )  THEN
     128                wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     129                IF ( wkomp > 0.0_wp )  THEN
    130130                   tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    131131                                         ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
    … …  
    182182!
    183183!--       x-direction
    184           ukomp = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
    185           IF ( ukomp > 0.0 )  THEN
     184          ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
     185          IF ( ukomp > 0.0_wp )  THEN
    186186             tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    187187                                         ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
    … …  
    192192!
    193193!--       y-direction
    194           vkomp = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
    195           IF ( vkomp > 0.0 )  THEN
     194          vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
     195          IF ( vkomp > 0.0_wp )  THEN
    196196             tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    197197                                         ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
    … …  
    202202!
    203203!--       z-direction
    204           wkomp = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
    205           IF ( wkomp > 0.0 )  THEN
     204          wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     205          IF ( wkomp > 0.0_wp )  THEN
    206206             tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    207207                                         ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)

  • palm/trunk/SOURCE/advec_u_pw.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    8989       REAL(wp)    ::  gv !:
    9090 
    91        gu = 2.0 * u_gtrans
    92        gv = 2.0 * v_gtrans
     91       gu = 2.0_wp * u_gtrans
     92       gv = 2.0_wp * v_gtrans
    9393       DO  i = nxlu, nxr
    9494          DO  j = nys, nyn
    9595             DO  k = nzb_u_inner(j,i)+1, nzt
    96                 tend(k,j,i) = tend(k,j,i) - 0.25 * (                           &
     96                tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                        &
    9797                         ( u(k,j,i+1) * ( u(k,j,i+1) + u(k,j,i) - gu )         &
    9898                         - u(k,j,i-1) * ( u(k,j,i) + u(k,j,i-1) - gu ) ) * ddx &
    … …  
    102102                         - u(k-1,j,i) * ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
    103103                                                                  * ddzw(k)    &
    104                                                    )
     104                                                      )
    105105             ENDDO
    106106          ENDDO
    … …  
    139139       REAL(wp)    ::  gv !:
    140140
    141        gu = 2.0 * u_gtrans
    142        gv = 2.0 * v_gtrans
     141       gu = 2.0_wp * u_gtrans
     142       gv = 2.0_wp * v_gtrans
    143143       DO  k = nzb_u_inner(j,i)+1, nzt
    144           tend(k,j,i) = tend(k,j,i) - 0.25 * (                                 &
     144          tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                              &
    145145                         ( u(k,j,i+1) * ( u(k,j,i+1) + u(k,j,i) - gu )         &
    146146                         - u(k,j,i-1) * ( u(k,j,i) + u(k,j,i-1) - gu ) ) * ddx &
    … …  
    150150                         - u(k-1,j,i) * ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
    151151                                                                  * ddzw(k)    &
    152                                              )
     152                                                )
    153153       ENDDO
    154154

  • palm/trunk/SOURCE/advec_u_up.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    9696!--             x-direction
    9797                ukomp = u(k,j,i) - u_gtrans
    98                 IF ( ukomp > 0.0 )  THEN
     98                IF ( ukomp > 0.0_wp )  THEN
    9999                   tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    100100                                         ( u(k,j,i) - u(k,j,i-1) ) * ddx
    … …  
    105105!
    106106!--             y-direction
    107                 vkomp = 0.25 * ( v(k,j,i)   + v(k,j+1,i) +                     &
     107                vkomp = 0.25_wp * ( v(k,j,i)   + v(k,j+1,i) +                  &
    108108                                 v(k,j,i-1) + v(k,j+1,i-1) ) - v_gtrans
    109                 IF ( vkomp > 0.0 )  THEN
     109                IF ( vkomp > 0.0_wp )  THEN
    110110                   tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    111111                                         ( u(k,j,i) - u(k,j-1,i) ) * ddy
    … …  
    116116!
    117117!--             z-direction
    118                 wkomp = 0.25 * ( w(k,j,i)   + w(k-1,j,i) +                     &
     118                wkomp = 0.25_wp * ( w(k,j,i)   + w(k-1,j,i) +                  &
    119119                                 w(k,j,i-1) + w(k-1,j,i-1) )
    120                 IF ( wkomp > 0.0 )  THEN
     120                IF ( wkomp > 0.0_wp )  THEN
    121121                   tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    122122                                         ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)
    … …  
    168168!--       x-direction
    169169          ukomp = u(k,j,i) - u_gtrans
    170           IF ( ukomp > 0.0 )  THEN
     170          IF ( ukomp > 0.0_wp )  THEN
    171171             tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    172172                                         ( u(k,j,i) - u(k,j,i-1) ) * ddx
    … …  
    177177!
    178178!--       y-direction
    179           vkomp = 0.25 * ( v(k,j,i) + v(k,j+1,i) + v(k,j,i-1) + v(k,j+1,i-1)  &
     179          vkomp = 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k,j,i-1) + v(k,j+1,i-1) &
    180180                         ) - v_gtrans
    181           IF ( vkomp > 0.0 )  THEN
     181          IF ( vkomp > 0.0_wp )  THEN
    182182             tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    183183                                         ( u(k,j,i) - u(k,j-1,i) ) * ddy
    … …  
    188188!
    189189!--       z-direction
    190           wkomp = 0.25 * ( w(k,j,i) + w(k-1,j,i) + w(k,j,i-1) + w(k-1,j,i-1) )
    191           IF ( wkomp > 0.0 )  THEN
     190          wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j,i-1) + w(k-1,j,i-1) )
     191          IF ( wkomp > 0.0_wp )  THEN
    192192             tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    193193                                         ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)

  • palm/trunk/SOURCE/advec_v_pw.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    9090 
    9191
    92        gu = 2.0 * u_gtrans
    93        gv = 2.0 * v_gtrans
     92       gu = 2.0_wp * u_gtrans
     93       gv = 2.0_wp * v_gtrans
    9494       DO  i = nxl, nxr
    9595          DO  j = nysv, nyn
    9696             DO  k = nzb_v_inner(j,i)+1, nzt
    97                 tend(k,j,i) = tend(k,j,i) - 0.25 * (                           &
     97                tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                        &
    9898                         ( v(k,j,i+1) * ( u(k,j-1,i+1) + u(k,j,i+1) - gu )     &
    9999                         - v(k,j,i-1) * ( u(k,j-1,i) + u(k,j,i) - gu ) ) * ddx &
    … …  
    103103                         - v(k-1,j,i) * ( w(k-1,j-1,i) + w(k-1,j,i) ) )        &
    104104                                                                  * ddzw(k)    &
    105                                                    )
     105                                                      )
    106106             ENDDO
    107107          ENDDO
    … …  
    141141
    142142
    143        gu = 2.0 * u_gtrans
    144        gv = 2.0 * v_gtrans
     143       gu = 2.0_wp * u_gtrans
     144       gv = 2.0_wp * v_gtrans
    145145       DO  k = nzb_v_inner(j,i)+1, nzt
    146           tend(k,j,i) = tend(k,j,i) - 0.25 * (                                 &
     146          tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                              &
    147147                         ( v(k,j,i+1) * ( u(k,j-1,i+1) + u(k,j,i+1) - gu )     &
    148148                         - v(k,j,i-1) * ( u(k,j-1,i) + u(k,j,i) - gu ) ) * ddx &
    … …  
    152152                         - v(k-1,j,i) * ( w(k-1,j-1,i) + w(k-1,j,i) ) )        &
    153153                                                                  * ddzw(k)    &
    154                                              )
     154                                                )
    155155       ENDDO
    156156 

  • palm/trunk/SOURCE/advec_v_up.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    9595!
    9696!--             x-direction
    97                 ukomp = 0.25 * ( u(k,j,i)   + u(k,j-1,i) +                     &
     97                ukomp = 0.25_wp * ( u(k,j,i)   + u(k,j-1,i) +                  &
    9898                                 u(k,j,i+1) + u(k,j-1,i+1) ) - u_gtrans
    99                 IF ( ukomp > 0.0 )  THEN
     99                IF ( ukomp > 0.0_wp )  THEN
    100100                   tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    101101                                         ( v(k,j,i) - v(k,j,i-1) ) * ddx
    … …  
    107107!--             y-direction
    108108                vkomp = v(k,j,i) - v_gtrans
    109                 IF ( vkomp > 0.0 )  THEN
     109                IF ( vkomp > 0.0_wp )  THEN
    110110                   tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    111111                                         ( v(k,j,i) - v(k,j-1,i) ) * ddy
    … …  
    116116!
    117117!--             z-direction
    118                 wkomp = 0.25 * ( w(k,j,i)  + w(k-1,j,i) +                      &
     118                wkomp = 0.25_wp * ( w(k,j,i)  + w(k-1,j,i) +                   &
    119119                                 w(k,j-1,i) + w(k-1,j-1,i) )
    120                 IF ( wkomp > 0.0 )  THEN
     120                IF ( wkomp > 0.0_wp )  THEN
    121121                   tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    122122                                         ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)
    … …  
    167167!
    168168!--       x-direction
    169           ukomp = 0.25 * ( u(k,j,i) + u(k,j-1,i) + u(k,j,i+1) + u(k,j-1,i+1)  &
     169          ukomp = 0.25_wp * ( u(k,j,i) + u(k,j-1,i) + u(k,j,i+1) + u(k,j-1,i+1) &
    170170                         ) - u_gtrans
    171           IF ( ukomp > 0.0 )  THEN
     171          IF ( ukomp > 0.0_wp )  THEN
    172172             tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    173173                                         ( v(k,j,i) - v(k,j,i-1) ) * ddx
    … …  
    179179!--       y-direction
    180180          vkomp = v(k,j,i) - v_gtrans
    181           IF ( vkomp > 0.0 )  THEN
     181          IF ( vkomp > 0.0_wp )  THEN
    182182             tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    183183                                         ( v(k,j,i) - v(k,j-1,i) ) * ddy
    … …  
    188188!
    189189!--       z-direction
    190           wkomp = 0.25 * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )
    191           IF ( wkomp > 0.0 )  THEN
     190          wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )
     191          IF ( wkomp > 0.0_wp )  THEN
    192192             tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    193193                                         ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)

  • palm/trunk/SOURCE/advec_w_pw.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    9090
    9191 
    92        gu = 2.0 * u_gtrans
    93        gv = 2.0 * v_gtrans
     92       gu = 2.0_wp * u_gtrans
     93       gv = 2.0_wp * v_gtrans
    9494       DO  i = nxl, nxr
    9595          DO  j = nys, nyn
    9696             DO  k = nzb_w_inner(j,i)+1, nzt
    97                 tend(k,j,i) = tend(k,j,i) - 0.25 * (                           &
     97                tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                        &
    9898                         ( w(k,j,i+1) * ( u(k+1,j,i+1) + u(k,j,i+1) - gu )     &
    9999                         - w(k,j,i-1) * ( u(k+1,j,i) + u(k,j,i) - gu ) ) * ddx &
    … …  
    103103                         - w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i) ) )            &
    104104                                                                  * ddzu(k+1)  &
    105                                                    )
     105                                                      )
    106106             ENDDO
    107107          ENDDO
    … …  
    140140       REAL(wp)    ::  gv !:
    141141
    142        gu = 2.0 * u_gtrans
    143        gv = 2.0 * v_gtrans
     142       gu = 2.0_wp * u_gtrans
     143       gv = 2.0_wp * v_gtrans
    144144       DO  k = nzb_w_inner(j,i)+1, nzt
    145           tend(k,j,i) = tend(k,j,i) - 0.25 * (                                 &
     145          tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                              &
    146146                         ( w(k,j,i+1) * ( u(k+1,j,i+1) + u(k,j,i+1) - gu )     &
    147147                         - w(k,j,i-1) * ( u(k+1,j,i) + u(k,j,i) - gu ) ) * ddx &

  • palm/trunk/SOURCE/advec_w_up.f90

    r1321 r1353  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
    … …  
    9494!
    9595!--             x-direction
    96                 ukomp = 0.25 * ( u(k,j,i)   + u(k,j,i+1) +                     &
     96                ukomp = 0.25_wp * ( u(k,j,i)   + u(k,j,i+1) +                  &
    9797                                 u(k+1,j,i) + u(k+1,j,i+1) ) - u_gtrans
    98                 IF ( ukomp > 0.0 )  THEN
     98                IF ( ukomp > 0.0_wp )  THEN
    9999                   tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    100100                                         ( w(k,j,i) - w(k,j,i-1) ) * ddx
    … …  
    105105!
    106106!--             y-direction
    107                 vkomp = 0.25 * ( v(k,j,i) + v(k,j+1,i) +                       &
     107                vkomp = 0.25_wp * ( v(k,j,i) + v(k,j+1,i) +                    &
    108108                                 v(k+1,j,i) + v(k+1,j+1,i) ) - v_gtrans
    109                 IF ( vkomp > 0.0 )  THEN
     109                IF ( vkomp > 0.0_wp )  THEN
    110110                   tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    111111                                         ( w(k,j,i) - w(k,j-1,i) ) * ddy
    … …  
    116116!
    117117!--             z-direction
    118                 IF ( w(k,j,i) > 0.0 )  THEN
     118                IF ( w(k,j,i) > 0.0_wp )  THEN
    119119                   tend(k,j,i) = tend(k,j,i) - w(k,j,i) *                      &
    120120                                         ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    … …  
    164164!
    165165!--       x-direction
    166           ukomp = 0.25 * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1)  &
     166          ukomp = 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) &
    167167                         ) - u_gtrans
    168           IF ( ukomp > 0.0 )  THEN
     168          IF ( ukomp > 0.0_wp )  THEN
    169169             tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    170170                                         ( w(k,j,i) - w(k,j,i-1) ) * ddx
    … …  
    175175!
    176176!--       y-direction
    177           vkomp = 0.25 * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i)  &
     177          vkomp = 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) &
    178178                         ) - v_gtrans
    179           IF ( vkomp > 0.0 )  THEN
     179          IF ( vkomp > 0.0_wp )  THEN
    180180             tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    181181                                         ( w(k,j,i) - w(k,j-1,i) ) * ddy
    … …  
    186186!
    187187!--       z-direction
    188           IF ( w(k,j,i) > 0.0 )  THEN
     188          IF ( w(k,j,i) > 0.0_wp )  THEN
    189189             tend(k,j,i) = tend(k,j,i) - w(k,j,i) *                            &
    190190                                         ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)

  • palm/trunk/SOURCE/advec_ws.f90

    r1323 r1353  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! REAL constants provided with KIND-attribute,
     23! module kinds added
     24! some formatting adjustments
    2325!
    2426! Former revisions:
    … …  
    116118
    117119    PRIVATE
    118     PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, &
    119              advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, &
     120    PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc,          &
     121             advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc,          &
    120122             ws_init, ws_statistics
    121123
    … …  
    187189
    188190       USE indices,                                                            &
    189            ONLY:  nyn, nys, nzb, nzt 
     191           ONLY:  nyn, nys, nzb, nzt
     192
     193       USE kinds
    190194       
    191195       USE pegrid
    … …  
    214218                    sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    215219
    216           sums_wsus_ws_l = 0.0
    217           sums_wsvs_ws_l = 0.0
    218           sums_us2_ws_l  = 0.0
    219           sums_vs2_ws_l  = 0.0
    220           sums_ws2_ws_l  = 0.0
     220          sums_wsus_ws_l = 0.0_wp
     221          sums_wsvs_ws_l = 0.0_wp
     222          sums_us2_ws_l  = 0.0_wp
     223          sums_vs2_ws_l  = 0.0_wp
     224          sums_ws2_ws_l  = 0.0_wp
    221225
    222226       ENDIF
    … …  
    225229
    226230          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    227           sums_wspts_ws_l = 0.0
     231          sums_wspts_ws_l = 0.0_wp
    228232
    229233          IF ( humidity  .OR.  passive_scalar )  THEN
    230234             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    231              sums_wsqs_ws_l = 0.0
     235             sums_wsqs_ws_l = 0.0_wp
    232236          ENDIF
    233237
    … …  
    236240             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    237241             ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    238              sums_wsqrs_ws_l = 0.0
    239              sums_wsnrs_ws_l = 0.0
     242             sums_wsqrs_ws_l = 0.0_wp
     243             sums_wsnrs_ws_l = 0.0_wp
    240244          ENDIF
    241245
    242246          IF ( ocean )  THEN
    243247             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    244              sums_wssas_ws_l = 0.0
     248             sums_wssas_ws_l = 0.0_wp
    245249          ENDIF
    246250
    … …  
    324328                  precipitation, ocean, ws_scheme_mom, ws_scheme_sca
    325329
     330       USE kinds
     331
    326332       USE statistics,                                                         &
    327333           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
    … …  
    335341!--    beginning of prognostic_equations.
    336342       IF ( ws_scheme_mom )  THEN
    337           sums_wsus_ws_l = 0.0
    338           sums_wsvs_ws_l = 0.0
    339           sums_us2_ws_l  = 0.0
    340           sums_vs2_ws_l  = 0.0
    341           sums_ws2_ws_l  = 0.0
     343          sums_wsus_ws_l = 0.0_wp
     344          sums_wsvs_ws_l = 0.0_wp
     345          sums_us2_ws_l  = 0.0_wp
     346          sums_vs2_ws_l  = 0.0_wp
     347          sums_ws2_ws_l  = 0.0_wp
    342348       ENDIF
    343349
    344350       IF ( ws_scheme_sca )  THEN
    345           sums_wspts_ws_l = 0.0
    346           IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0
     351          sums_wspts_ws_l = 0.0_wp
     352          IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0_wp
    347353          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
    348354               precipitation )  THEN
    349              sums_wsqrs_ws_l = 0.0
    350              sums_wsnrs_ws_l = 0.0
     355             sums_wsqrs_ws_l = 0.0_wp
     356             sums_wsnrs_ws_l = 0.0_wp
    351357          ENDIF
    352           IF ( ocean )  sums_wssas_ws_l = 0.0
     358          IF ( ocean )  sums_wssas_ws_l = 0.0_wp
    353359
    354360       ENDIF
    … …  
    446452
    447453             v_comp                  = v(k,j,i) - v_gtrans
    448              swap_flux_y_local(k,tn) = v_comp * (                             &
    449                                                   ( 37.0 * ibit5 * adv_sca_5  &
    450                                                +     7.0 * ibit4 * adv_sca_3  &
    451                                                +           ibit3 * adv_sca_1  &
    452                                                   ) *                         &
     454             swap_flux_y_local(k,tn) = v_comp *         (                     &
     455                                               ( 37.0_wp * ibit5 * adv_sca_5  &
     456                                            +     7.0_wp * ibit4 * adv_sca_3  &
     457                                            +              ibit3 * adv_sca_1  &
     458                                               ) *                            &
    453459                                           ( sk(k,j,i)  + sk(k,j-1,i)     )   &
    454                                             -     (  8.0 * ibit5 * adv_sca_5  &
    455                                                +           ibit4 * adv_sca_3  &
    456                                                    ) *                        &
     460                                         -     (  8.0_wp * ibit5 * adv_sca_5  &
     461                                            +              ibit4 * adv_sca_3  &
     462                                                ) *                           &
    457463                                           ( sk(k,j+1,i) + sk(k,j-2,i)    )   &
    458                                            +      (        ibit5 * adv_sca_5  &
    459                                                   ) *                         &
     464                                         +     (           ibit5 * adv_sca_5  &
     465                                               ) *                            &
    460466                                           ( sk(k,j+2,i) + sk(k,j-3,i)    )   &
    461                                                 )
     467                                                        )
    462468
    463469             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
    464                                                   ( 10.0 * ibit5 * adv_sca_5  &
    465                                                +     3.0 * ibit4 * adv_sca_3  &
    466                                                +           ibit3 * adv_sca_1  &
    467                                                   ) *                         &
    468                                             ( sk(k,j,i)   - sk(k,j-1,i)    )  &
    469                                            -      (  5.0 * ibit5 * adv_sca_5  &
    470                                                +           ibit4 * adv_sca_3  &
     470                                               ( 10.0_wp * ibit5 * adv_sca_5  &
     471                                            +     3.0_wp * ibit4 * adv_sca_3  &
     472                                            +              ibit3 * adv_sca_1  &
     473                                               ) *                            &
     474                                            ( sk(k,j,i)   - sk(k,j-1,i)  )    &
     475                                        -      (  5.0_wp * ibit5 * adv_sca_5  &
     476                                            +              ibit4 * adv_sca_3  &
    471477                                            ) *                               &
    472478                                            ( sk(k,j+1,i) - sk(k,j-2,i)  )    &
    473                                            +      (        ibit5 * adv_sca_5  &
    474                                                   ) *                         &
    475                                             ( sk(k,j+2,i) - sk(k,j-3,i) )     &
     479                                        +      (           ibit5 * adv_sca_5  &
     480                                               ) *                            &
     481                                            ( sk(k,j+2,i) - sk(k,j-3,i)  )    &
    476482                                                        )
    477483
    … …  
    482488
    483489             v_comp                  = v(k,j,i) - v_gtrans
    484              swap_flux_y_local(k,tn) = v_comp * (                            &
    485                                     37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
    486                                   -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
    487                                   +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
    488                                              ) * adv_sca_5
    489               swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                   &
    490                                     10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
    491                                   -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
    492                                   +          sk(k,j+2,i) - sk(k,j-3,i)       &
    493                                                         ) * adv_sca_5
     490             swap_flux_y_local(k,tn) = v_comp * (                             &
     491                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )   &
     492                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )   &
     493                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )   &
     494                                                ) * adv_sca_5
     495              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                    &
     496                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )   &
     497                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )   &
     498                                  +             sk(k,j+2,i) - sk(k,j-3,i)     &
     499                                                         ) * adv_sca_5
    494500
    495501          ENDDO
    … …  
    508514             u_comp                     = u(k,j,i) - u_gtrans
    509515             swap_flux_x_local(k,j,tn) = u_comp * (                           &
    510                                                   ( 37.0 * ibit2 * adv_sca_5  &
    511                                                +     7.0 * ibit1 * adv_sca_3  &
    512                                                +           ibit0 * adv_sca_1  &
    513                                                   ) *                         &
     516                                               ( 37.0_wp * ibit2 * adv_sca_5  &
     517                                            +     7.0_wp * ibit1 * adv_sca_3  &
     518                                            +              ibit0 * adv_sca_1  &
     519                                               ) *                            &
    514520                                            ( sk(k,j,i)   + sk(k,j,i-1)    )  &
    515                                            -      (  8.0 * ibit2 * adv_sca_5  &
    516                                                +           ibit1 * adv_sca_3  &
    517                                                   ) *                         &
     521                                        -      (  8.0_wp * ibit2 * adv_sca_5  &
     522                                            +              ibit1 * adv_sca_3  &
     523                                               ) *                            &
    518524                                            ( sk(k,j,i+1) + sk(k,j,i-2)    )  &
    519                                            +      (        ibit2 * adv_sca_5  &
    520                                                   ) *                         &
     525                                        +      (           ibit2 * adv_sca_5  &
     526                                               ) *                            &
    521527                                            ( sk(k,j,i+2) + sk(k,j,i-3)    )  &
    522528                                                  )
    523529
    524530              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
    525                                                   ( 10.0 * ibit2 * adv_sca_5  &
    526                                                +     3.0 * ibit1 * adv_sca_3  &
    527                                                +           ibit0 * adv_sca_1  &
    528                                                   ) *                         &
     531                                               ( 10.0_wp * ibit2 * adv_sca_5  &
     532                                            +     3.0_wp * ibit1 * adv_sca_3  &
     533                                            +              ibit0 * adv_sca_1  &
     534                                               ) *                            &
    529535                                            ( sk(k,j,i)   - sk(k,j,i-1)    )  &
    530                                            -      (  5.0 * ibit2 * adv_sca_5  &
    531                                                +           ibit1 * adv_sca_3  &
    532                                                   ) *                         &
    533                                             ( sk(k,j,i+1) - sk(k,j,i-2)  )    &
    534                                            +      (        ibit2 * adv_sca_5  &
    535                                                   ) *                         &
    536                                             ( sk(k,j,i+2) - sk(k,j,i-3) )     &
     536                                        -      (  5.0_wp * ibit2 * adv_sca_5  &
     537                                            +              ibit1 * adv_sca_3  &
     538                                               ) *                            &
     539                                            ( sk(k,j,i+1) - sk(k,j,i-2)    )  &
     540                                        +      (           ibit2 * adv_sca_5  &
     541                                               ) *                            &
     542                                            ( sk(k,j,i+2) - sk(k,j,i-3)    )  &
    537543                                                           )
    538544
    … …  
    542548
    543549             u_comp                 = u(k,j,i) - u_gtrans
    544              swap_flux_x_local(k,j,tn) = u_comp * (                          &
    545                                       37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )  &
    546                                     -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
    547                                     +        ( sk(k,j,i+2) + sk(k,j,i-3) )  &
     550             swap_flux_x_local(k,j,tn) = u_comp * (                           &
     551                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
     552                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
     553                                    +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
    548554                                                  ) * adv_sca_5
    549555
    550              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
    551                                       10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )  &
    552                                     -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
    553                                     +        ( sk(k,j,i+2) - sk(k,j,i-3) )  &
     556             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                   &
     557                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
     558                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
     559                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
    554560                                                          ) * adv_sca_5
    555561
    … …  
    558564       ENDIF
    559565
    560        flux_t(0) = 0.0
    561        diss_t(0) = 0.0
    562        flux_d    = 0.0
    563        diss_d    = 0.0
     566       flux_t(0) = 0.0_wp
     567       diss_t(0) = 0.0_wp
     568       flux_d    = 0.0_wp
     569       diss_d    = 0.0_wp
    564570!       
    565571!--    Now compute the fluxes and tendency terms for the horizontal and
    … …  
    576582
    577583          u_comp    = u(k,j,i+1) - u_gtrans
    578           flux_r(k) = u_comp * (                                            &
    579                      ( 37.0 * ibit2 * adv_sca_5                             &
    580                   +     7.0 * ibit1 * adv_sca_3                             &
    581                   +           ibit0 * adv_sca_1                             &
    582                      ) *                                                    &
    583                              ( sk(k,j,i+1) + sk(k,j,i)   )                  &
    584               -      (  8.0 * ibit2 * adv_sca_5                             &
    585                   +           ibit1 * adv_sca_3                             &
    586                      ) *                                                    &
    587                              ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
    588               +      (        ibit2 * adv_sca_5                             &
    589                      ) *                                                    &
    590                              ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
     584          flux_r(k) = u_comp * (                                              &
     585                     ( 37.0_wp * ibit2 * adv_sca_5                            &
     586                  +     7.0_wp * ibit1 * adv_sca_3                            &
     587                  +              ibit0 * adv_sca_1                            &
     588                     ) *                                                      &
     589                             ( sk(k,j,i+1) + sk(k,j,i)   )                    &
     590              -      (  8.0_wp * ibit2 * adv_sca_5                            &
     591                  +              ibit1 * adv_sca_3                            &
     592                     ) *                                                      &
     593                             ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
     594              +      (           ibit2 * adv_sca_5                            &
     595                     ) *                                                      &
     596                             ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
    591597                               )
    592598
    593           diss_r(k) = -ABS( u_comp ) * (                                    &
    594                      ( 10.0 * ibit2 * adv_sca_5                             &
    595                   +     3.0 * ibit1 * adv_sca_3                             &
    596                   +           ibit0 * adv_sca_1                             &
    597                      ) *                                                    &
    598                              ( sk(k,j,i+1) - sk(k,j,i)  )                   &
    599               -      (  5.0 * ibit2 * adv_sca_5                             &
    600                   +           ibit1 * adv_sca_3                             &
    601                      ) *                                                    &
    602                              ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
    603               +      (        ibit2 * adv_sca_5                             &
    604                      ) *                                                    &
    605                              ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
     599          diss_r(k) = -ABS( u_comp ) * (                                      &
     600                     ( 10.0_wp * ibit2 * adv_sca_5                            &
     601                  +     3.0_wp * ibit1 * adv_sca_3                            &
     602                  +              ibit0 * adv_sca_1                            &
     603                     ) *                                                      &
     604                             ( sk(k,j,i+1) - sk(k,j,i)  )                     &
     605              -      (  5.0_wp * ibit2 * adv_sca_5                            &
     606                  +              ibit1 * adv_sca_3                            &
     607                     ) *                                                      &
     608                             ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
     609              +      (           ibit2 * adv_sca_5                            &
     610                     ) *                                                      &
     611                             ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
    606612                                       )
    607613
    … …  
    611617
    612618          v_comp    = v(k,j+1,i) - v_gtrans
    613           flux_n(k) = v_comp * (                                            &
    614                      ( 37.0 * ibit5 * adv_sca_5                             &
    615                   +     7.0 * ibit4 * adv_sca_3                             &
    616                   +           ibit3 * adv_sca_1                             &
    617                      ) *                                                    &
    618                              ( sk(k,j+1,i) + sk(k,j,i)   )                  &
    619               -      (  8.0 * ibit5 * adv_sca_5                             &
    620                   +           ibit4 * adv_sca_3                             &
    621                      ) *                                                    &
    622                              ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
    623               +      (        ibit5 * adv_sca_5                             &
    624                      ) *                                                    &
    625                              ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
     619          flux_n(k) = v_comp * (                                              &
     620                     ( 37.0_wp * ibit5 * adv_sca_5                            &
     621                  +     7.0_wp * ibit4 * adv_sca_3                            &
     622                  +              ibit3 * adv_sca_1                            &
     623                     ) *                                                      &
     624                             ( sk(k,j+1,i) + sk(k,j,i)   )                    &
     625              -      (  8.0_wp * ibit5 * adv_sca_5                            &
     626                  +              ibit4 * adv_sca_3                            &
     627                     ) *                                                      &
     628                             ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
     629              +      (           ibit5 * adv_sca_5                            &
     630                     ) *                                                      &
     631                             ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
    626632                               )
    627633
    628           diss_n(k) = -ABS( v_comp ) * (                                    &
    629                      ( 10.0 * ibit5 * adv_sca_5                             &
    630                   +     3.0 * ibit4 * adv_sca_3                             &
    631                   +           ibit3 * adv_sca_1                             &
    632                      ) *                                                    &
    633                              ( sk(k,j+1,i) - sk(k,j,i)    )                 &
    634               -      (  5.0 * ibit5 * adv_sca_5                             &
    635                   +           ibit4 * adv_sca_3                             &
    636                      ) *                                                    &
    637                              ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
    638               +      (        ibit5 * adv_sca_5                             &
    639                      ) *                                                    &
    640                              ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
     634          diss_n(k) = -ABS( v_comp ) * (                                      &
     635                     ( 10.0_wp * ibit5 * adv_sca_5                            &
     636                  +     3.0_wp * ibit4 * adv_sca_3                            &
     637                  +              ibit3 * adv_sca_1                            &
     638                     ) *                                                      &
     639                             ( sk(k,j+1,i) - sk(k,j,i)   )                    &
     640              -      (  5.0_wp * ibit5 * adv_sca_5                            &
     641                  +              ibit4 * adv_sca_3                            &
     642                     ) *                                                      &
     643                             ( sk(k,j+2,i) - sk(k,j-1,i) )                    &
     644              +      (           ibit5 * adv_sca_5                            &
     645                     ) *                                                      &
     646                             ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
    641647                                       )
    642648!
    … …  
    652658
    653659
    654           flux_t(k) = w(k,j,i) * (                                          &
    655                      ( 37.0 * ibit8 * adv_sca_5                             &
    656                   +     7.0 * ibit7 * adv_sca_3                             &
    657                   +           ibit6 * adv_sca_1                             &
    658                      ) *                                                    &
    659                              ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
    660               -      (  8.0 * ibit8 * adv_sca_5                             &
    661                   +           ibit7 * adv_sca_3                             &
    662                      ) *                                                    &
    663                              ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
    664               +      (        ibit8 * adv_sca_5                             &
    665                      ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
     660          flux_t(k) = w(k,j,i) * (                                            &
     661                     ( 37.0_wp * ibit8 * adv_sca_5                            &
     662                  +     7.0_wp * ibit7 * adv_sca_3                            &
     663                  +              ibit6 * adv_sca_1                            &
     664                     ) *                                                      &
     665                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
     666              -      (  8.0_wp * ibit8 * adv_sca_5                            &
     667                  +              ibit7 * adv_sca_3                            &
     668                     ) *                                                      &
     669                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
     670              +      (           ibit8 * adv_sca_5                            &
     671                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
    666672                                 )
    667673
    668           diss_t(k) = -ABS( w(k,j,i) ) * (                                  &
    669                      ( 10.0 * ibit8 * adv_sca_5                             &
    670                   +     3.0 * ibit7 * adv_sca_3                             &
    671                   +           ibit6 * adv_sca_1                             &
    672                      ) *                                                    &
    673                              ( sk(k+1,j,i)   - sk(k,j,i)    )               &
    674               -      (  5.0 * ibit8 * adv_sca_5                             &
    675                   +           ibit7 * adv_sca_3                             &
    676                      ) *                                                    &
    677                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
    678               +      (        ibit8 * adv_sca_5                             &
    679                      ) *                                                    &
    680                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
     674          diss_t(k) = -ABS( w(k,j,i) ) * (                                    &
     675                     ( 10.0_wp * ibit8 * adv_sca_5                            &
     676                  +     3.0_wp * ibit7 * adv_sca_3                            &
     677                  +              ibit6 * adv_sca_1                            &
     678                     ) *                                                      &
     679                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
     680              -      (  5.0_wp * ibit8 * adv_sca_5                            &
     681                  +              ibit7 * adv_sca_3                            &
     682                     ) *                                                      &
     683                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
     684              +      (           ibit8 * adv_sca_5                            &
     685                     ) *                                                      &
     686                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    681687                                         )
    682688!
    … …  
    712718
    713719          u_comp    = u(k,j,i+1) - u_gtrans
    714           flux_r(k) = u_comp * (                                            &
    715                       37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                  &
    716                     -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
    717                     +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    718           diss_r(k) = -ABS( u_comp ) * (                                    &
    719                       10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                  &
    720                     -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
    721                     +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     720          flux_r(k) = u_comp * (                                              &
     721                      37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
     722                    -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
     723                    +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     724          diss_r(k) = -ABS( u_comp ) * (                                      &
     725                      10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
     726                    -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
     727                    +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    722728
    723729          v_comp    = v(k,j+1,i) - v_gtrans
    724           flux_n(k) = v_comp * (                                            &
    725                       37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                  &
    726                     -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
    727                     +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    728           diss_n(k) = -ABS( v_comp ) * (                                    &
    729                       10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                  &
    730                     -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                  &
    731                     +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     730          flux_n(k) = v_comp * (                                              &
     731                      37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
     732                    -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
     733                    +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     734          diss_n(k) = -ABS( v_comp ) * (                                      &
     735                      10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
     736                    -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
     737                    +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    732738!
    733739!--       k index has to be modified near bottom and top, else array
    … …  
    742748
    743749
    744           flux_t(k) = w(k,j,i) * (                                          &
    745                      ( 37.0 * ibit8 * adv_sca_5                             &
    746                   +     7.0 * ibit7 * adv_sca_3                             &
    747                   +           ibit6 * adv_sca_1                             &
    748                      ) *                                                    &
    749                              ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
    750               -      (  8.0 * ibit8 * adv_sca_5                             &
    751                   +           ibit7 * adv_sca_3                             &
    752                      ) *                                                    &
    753                              ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
    754               +      (        ibit8 * adv_sca_5                             &
    755                      ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
     750          flux_t(k) = w(k,j,i) * (                                            &
     751                    ( 37.0_wp * ibit8 * adv_sca_5                             &
     752                 +     7.0_wp * ibit7 * adv_sca_3                             &
     753                 +              ibit6 * adv_sca_1                             &
     754                    ) *                                                       &
     755                             ( sk(k+1,j,i)  + sk(k,j,i)   )                   &
     756              -     (  8.0_wp * ibit8 * adv_sca_5                             &
     757                  +              ibit7 * adv_sca_3                            &
     758                    ) *                                                       &
     759                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                   &
     760              +     (           ibit8 * adv_sca_5                             &
     761                    ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                   &
    756762                                 )
    757763
    758           diss_t(k) = -ABS( w(k,j,i) ) * (                                  &
    759                      ( 10.0 * ibit8 * adv_sca_5                             &
    760                   +     3.0 * ibit7 * adv_sca_3                             &
    761                   +           ibit6 * adv_sca_1                             &
    762                      ) *                                                    &
    763                              ( sk(k+1,j,i)   - sk(k,j,i)    )               &
    764               -      (  5.0 * ibit8 * adv_sca_5                             &
    765                   +           ibit7 * adv_sca_3                             &
    766                      ) *                                                    &
    767                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
    768               +      (        ibit8 * adv_sca_5                             &
    769                      ) *                                                    &
    770                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
     764          diss_t(k) = -ABS( w(k,j,i) ) * (                                    &
     765                    ( 10.0_wp * ibit8 * adv_sca_5                             &
     766                 +     3.0_wp * ibit7 * adv_sca_3                             &
     767                 +              ibit6 * adv_sca_1                             &
     768                    ) *                                                       &
     769                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
     770              -     (  5.0_wp * ibit8 * adv_sca_5                             &
     771                 +              ibit7 * adv_sca_3                             &
     772                    ) *                                                       &
     773                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
     774              +     (           ibit8 * adv_sca_5                             &
     775                    ) *                                                       &
     776                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    771777                                         )
    772778!
    … …  
    803809
    804810             DO  k = nzb, nzt
    805                 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                &
    806                                        ( flux_t(k) + diss_t(k) )               &
     811                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +               &
     812                                       ( flux_t(k) + diss_t(k) )              &
    807813                                 * weight_substep(intermediate_timestep_count)
    808814             ENDDO
    … …  
    811817
    812818             DO  k = nzb, nzt
    813                 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                &
    814                                        ( flux_t(k) + diss_t(k) )               &
     819                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +               &
     820                                       ( flux_t(k) + diss_t(k) )              &
    815821                                 * weight_substep(intermediate_timestep_count)
    816822             ENDDO
    … …  
    819825
    820826             DO  k = nzb, nzt
    821                 sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                 &
    822                                       ( flux_t(k) + diss_t(k) )                &
     827                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                &
     828                                      ( flux_t(k) + diss_t(k) )               &
    823829                                 * weight_substep(intermediate_timestep_count)
    824830             ENDDO
    … …  
    827833
    828834             DO  k = nzb, nzt
    829                 sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +               &
    830                                       ( flux_t(k) + diss_t(k) )                &
     835                sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +              &
     836                                      ( flux_t(k) + diss_t(k) )               &
    831837                                 * weight_substep(intermediate_timestep_count)
    832838             ENDDO
    … …  
    835841
    836842             DO  k = nzb, nzt
    837                 sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +               &
    838                                       ( flux_t(k) + diss_t(k) )                &
     843                sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +              &
     844                                      ( flux_t(k) + diss_t(k) )               &
    839845                                 * weight_substep(intermediate_timestep_count)
    840846             ENDDO
    … …  
    852858    SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
    853859
    854        USE arrays_3d,                                                          &
     860       USE arrays_3d,                                                         &
    855861           ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w
    856862
    857        USE constants,                                                          &
     863       USE constants,                                                         &
    858864           ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
    859865
    860        USE control_parameters,                                                 &
     866       USE control_parameters,                                                &
    861867           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    862868
    863        USE grid_variables,                                                     &
     869       USE grid_variables,                                                    &
    864870           ONLY:  ddx, ddy
    865871
    866        USE indices,                                                            &
     872       USE indices,                                                           &
    867873           ONLY:  nxl, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0
    868874
    869875       USE kinds
    870876
    871        USE statistics,                                                         &
     877       USE statistics,                                                        &
    872878           ONLY:  hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
    873879
    … …  
    909915       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_comp !:
    910916
    911        gu = 2.0 * u_gtrans
    912        gv = 2.0 * v_gtrans
     917       gu = 2.0_wp * u_gtrans
     918       gv = 2.0_wp * v_gtrans
    913919!
    914920!--    Compute southside fluxes for the respective boundary of PE
    … …  
    923929             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    924930             flux_s_u(k,tn) = v_comp * (                                      &
    925                             ( 37.0 * ibit14 * adv_mom_5                       &
    926                          +     7.0 * ibit13 * adv_mom_3                       &
    927                          +           ibit12 * adv_mom_1                       &
     931                            ( 37.0_wp * ibit14 * adv_mom_5                    &
     932                         +     7.0_wp * ibit13 * adv_mom_3                    &
     933                         +              ibit12 * adv_mom_1                    &
    928934                            ) *                                               &
    929                                      ( u(k,j,i)   + u(k,j-1,i) )              &
    930                      -      (  8.0 * ibit14 * adv_mom_5                       &
    931                          +           ibit13 * adv_mom_3                       &
     935                                        ( u(k,j,i)   + u(k,j-1,i) )           &
     936                     -      (  8.0_wp * ibit14 * adv_mom_5                    &
     937                         +              ibit13 * adv_mom_3                    &
    932938                            ) *                                               &
    933                                      ( u(k,j+1,i) + u(k,j-2,i) )              &
    934                      +      (        ibit14 * adv_mom_5                       &
     939                                        ( u(k,j+1,i) + u(k,j-2,i) )           &
     940                     +      (           ibit14 * adv_mom_5                    &
    935941                            ) *                                               &
    936                                      ( u(k,j+2,i) + u(k,j-3,i) )              &
     942                                        ( u(k,j+2,i) + u(k,j-3,i) )           &
    937943                                       )
    938944
    939945             diss_s_u(k,tn) = - ABS ( v_comp ) * (                            &
    940                             ( 10.0 * ibit14 * adv_mom_5                       &
    941                          +     3.0 * ibit13 * adv_mom_3                       &
    942                          +           ibit12 * adv_mom_1                       &
     946                            ( 10.0_wp * ibit14 * adv_mom_5                    &
     947                         +     3.0_wp * ibit13 * adv_mom_3                    &
     948                         +              ibit12 * adv_mom_1                    &
    943949                            ) *                                               &
    944                                      ( u(k,j,i)   - u(k,j-1,i) )              &
    945                      -      (  5.0 * ibit14 * adv_mom_5                       &
    946                          +           ibit13 * adv_mom_3                       &
     950                                        ( u(k,j,i)   - u(k,j-1,i) )           &
     951                     -      (  5.0_wp * ibit14 * adv_mom_5                    &
     952                         +              ibit13 * adv_mom_3                    &
    947953                            ) *                                               &
    948                                      ( u(k,j+1,i) - u(k,j-2,i) )              &
    949                      +      (        ibit14 * adv_mom_5                       &
     954                                        ( u(k,j+1,i) - u(k,j-2,i) )           &
     955                     +      (           ibit14 * adv_mom_5                    &
    950956                            ) *                                               &
    951                                      ( u(k,j+2,i) - u(k,j-3,i) )              &
     957                                        ( u(k,j+2,i) - u(k,j-3,i) )           &
    952958                                                 )
    953959
    … …  
    958964             v_comp         = v(k,j,i) + v(k,j,i-1) - gv
    959965             flux_s_u(k,tn) = v_comp * (                                      &
    960                            37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
    961                          -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
    962                          +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
     966                           37.0_wp * ( u(k,j,i) + u(k,j-1,i)   )              &
     967                         -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )              &
     968                         +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    963969             diss_s_u(k,tn) = - ABS(v_comp) * (                               &
    964                            10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
    965                          -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
    966                          +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
     970                           10.0_wp * ( u(k,j,i) - u(k,j-1,i)   )              &
     971                         -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )              &
     972                         +           ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
    967973
    968974          ENDDO
    … …  
    980986
    981987             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
    982              flux_l_u(k,j,tn) = u_comp_l * (                                   &
    983                               ( 37.0 * ibit11 * adv_mom_5                      &
    984                            +     7.0 * ibit10 * adv_mom_3                      &
    985                            +           ibit9  * adv_mom_1                      &
    986                               ) *                                              &
    987                                      ( u(k,j,i)   + u(k,j,i-1) )               &
    988                        -      (  8.0 * ibit11 * adv_mom_5                      &
    989                            +           ibit10 * adv_mom_3                      &
    990                               ) *                                              &
    991                                      ( u(k,j,i+1) + u(k,j,i-2) )               &
    992                        +      (        ibit11 * adv_mom_5                      &
    993                               ) *                                              &
    994                                      ( u(k,j,i+2) + u(k,j,i-3) )               &
     988             flux_l_u(k,j,tn) = u_comp_l * (                                  &
     989                              ( 37.0_wp * ibit11 * adv_mom_5                  &
     990                           +     7.0_wp * ibit10 * adv_mom_3                  &
     991                           +              ibit9  * adv_mom_1                  &
     992                              ) *                                             &
     993                                          ( u(k,j,i)   + u(k,j,i-1) )         &
     994                       -      (  8.0_wp * ibit11 * adv_mom_5                  &
     995                           +              ibit10 * adv_mom_3                  &
     996                              ) *                                             &
     997                                          ( u(k,j,i+1) + u(k,j,i-2) )         &
     998                       +      (           ibit11 * adv_mom_5                  &
     999                              ) *                                             &
     1000                                          ( u(k,j,i+2) + u(k,j,i-3) )         &
    9951001                                           )
    9961002
    997               diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                         &
    998                               ( 10.0 * ibit11 * adv_mom_5                      &
    999                            +     3.0 * ibit10 * adv_mom_3                      &
    1000                            +           ibit9  * adv_mom_1                      &
    1001                               ) *                                              &
    1002                                      ( u(k,j,i)   - u(k,j,i-1) )               &
    1003                        -      (  5.0 * ibit11 * adv_mom_5                      &
    1004                            +           ibit10 * adv_mom_3                      &
    1005                               ) *                                              &
    1006                                      ( u(k,j,i+1) - u(k,j,i-2) )               &
    1007                        +      (        ibit11 * adv_mom_5                      &
    1008                               ) *                                              &
    1009                                      ( u(k,j,i+2) - u(k,j,i-3) )               &
     1003              diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                        &
     1004                              ( 10.0_wp * ibit11 * adv_mom_5                  &
     1005                           +     3.0_wp * ibit10 * adv_mom_3                  &
     1006                           +              ibit9  * adv_mom_1                  &
     1007                              ) *                                             &
     1008                                        ( u(k,j,i)   - u(k,j,i-1) )           &
     1009                       -      (  5.0_wp * ibit11 * adv_mom_5                  &
     1010                           +              ibit10 * adv_mom_3                  &
     1011                              ) *                                             &
     1012                                        ( u(k,j,i+1) - u(k,j,i-2) )           &
     1013                       +      (           ibit11 * adv_mom_5                  &
     1014                              ) *                                             &
     1015                                        ( u(k,j,i+2) - u(k,j,i-3) )           &
    10101016                                                     )
    10111017
    … …  
    10161022             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
    10171023             flux_l_u(k,j,tn) = u_comp_l * (                                   &
    1018                              37.0 * ( u(k,j,i) + u(k,j,i-1)   )                &
    1019                            -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )                &
    1020                            +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
     1024                             37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )             &
     1025                           -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )             &
     1026                           +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
    10211027             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                            &
    1022                              10.0 * ( u(k,j,i) - u(k,j,i-1)   )                &
    1023                            -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )                &
    1024                            +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
     1028                             10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )             &
     1029                           -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )             &
     1030                           +           ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
    10251031
    10261032          ENDDO
    … …  
    10281034       ENDIF
    10291035
    1030        flux_t(0) = 0.0
    1031        diss_t(0) = 0.0
    1032        flux_d    = 0.0
    1033        diss_d    = 0.0
     1036       flux_t(0) = 0.0_wp
     1037       diss_t(0) = 0.0_wp
     1038       flux_d    = 0.0_wp
     1039       diss_d    = 0.0_wp
    10341040!
    10351041!--    Now compute the fluxes tendency terms for the horizontal and
    … …  
    10421048
    10431049          u_comp(k) = u(k,j,i+1) + u(k,j,i)
    1044           flux_r(k) = ( u_comp(k) - gu ) * (                                 &
    1045                      ( 37.0 * ibit11 * adv_mom_5                             &
    1046                   +     7.0 * ibit10 * adv_mom_3                             &
    1047                   +           ibit9  * adv_mom_1                             &
    1048                      ) *                                                     &
    1049                                  ( u(k,j,i+1) + u(k,j,i)   )                 &
    1050               -      (  8.0 * ibit11 * adv_mom_5                             &
    1051                   +           ibit10 * adv_mom_3                             &
    1052                      ) *                                                     &
    1053                                  ( u(k,j,i+2) + u(k,j,i-1) )                 &
    1054               +      (        ibit11 * adv_mom_5                             &
    1055                      ) *                                                     &
    1056                                  ( u(k,j,i+3) + u(k,j,i-2) )                 &
     1050          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
     1051                     ( 37.0_wp * ibit11 * adv_mom_5                           &
     1052                  +     7.0_wp * ibit10 * adv_mom_3                           &
     1053                  +              ibit9  * adv_mom_1                           &
     1054                     ) *                                                      &
     1055                                    ( u(k,j,i+1) + u(k,j,i)   )               &
     1056              -      (  8.0_wp * ibit11 * adv_mom_5                           &
     1057                  +              ibit10 * adv_mom_3                           &
     1058                     ) *                                                      &
     1059                                    ( u(k,j,i+2) + u(k,j,i-1) )               &
     1060              +      (           ibit11 * adv_mom_5                           &
     1061                     ) *                                                      &
     1062                                    ( u(k,j,i+3) + u(k,j,i-2) )               &
    10571063                                           )
    10581064
    1059           diss_r(k) = - ABS( u_comp(k) - gu ) * (                            &
    1060                      ( 10.0 * ibit11 * adv_mom_5                             &
    1061                   +     3.0 * ibit10 * adv_mom_3                             &
    1062                   +           ibit9  * adv_mom_1                             &
    1063                      ) *                                                     &
    1064                                  ( u(k,j,i+1) - u(k,j,i)  )                  &
    1065               -      (  5.0 * ibit11 * adv_mom_5                             &
    1066                   +           ibit10 * adv_mom_3                             &
    1067                      ) *                                                     &
    1068                                  ( u(k,j,i+2) - u(k,j,i-1) )                 &
    1069               +      (        ibit11 * adv_mom_5                             &
    1070                      ) *                                                     &
    1071                                  ( u(k,j,i+3) - u(k,j,i-2) )                 &
     1065          diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
     1066                     ( 10.0_wp * ibit11 * adv_mom_5                           &
     1067                  +     3.0_wp * ibit10 * adv_mom_3                           &
     1068                  +              ibit9  * adv_mom_1                           &
     1069                     ) *                                                      &
     1070                                    ( u(k,j,i+1) - u(k,j,i)   )               &
     1071              -      (  5.0_wp * ibit11 * adv_mom_5                           &
     1072                  +              ibit10 * adv_mom_3                           &
     1073                     ) *                                                      &
     1074                                    ( u(k,j,i+2) - u(k,j,i-1) )               &
     1075              +      (           ibit11 * adv_mom_5                           &
     1076                     ) *                                                      &
     1077                                    ( u(k,j,i+3) - u(k,j,i-2) )               &
    10721078                                                )
    10731079
    … …  
    10771083
    10781084          v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    1079           flux_n(k) = v_comp * (                                             &
    1080                      ( 37.0 * ibit14 * adv_mom_5                             &
    1081                   +     7.0 * ibit13 * adv_mom_3                             &
    1082                   +           ibit12 * adv_mom_1                             &
    1083                      ) *                                                     &
    1084                                  ( u(k,j+1,i) + u(k,j,i)   )                 &
    1085               -      (  8.0 * ibit14 * adv_mom_5                             &
    1086                   +           ibit13 * adv_mom_3                             &
    1087                      ) *                                                     &
    1088                                  ( u(k,j+2,i) + u(k,j-1,i) )                 &
    1089               +      (        ibit14 * adv_mom_5                             &
    1090                      ) *                                                     &
    1091                                  ( u(k,j+3,i) + u(k,j-2,i) )                 &
     1085          flux_n(k) = v_comp * (                                              &
     1086                     ( 37.0_wp * ibit14 * adv_mom_5                           &
     1087                  +     7.0_wp * ibit13 * adv_mom_3                           &
     1088                  +              ibit12 * adv_mom_1                           &
     1089                     ) *                                                      &
     1090                                    ( u(k,j+1,i) + u(k,j,i)   )               &
     1091              -      (  8.0_wp * ibit14 * adv_mom_5                           &
     1092                  +              ibit13 * adv_mom_3                           &
     1093                     ) *                                                      &
     1094                                    ( u(k,j+2,i) + u(k,j-1,i) )               &
     1095              +      (           ibit14 * adv_mom_5                           &
     1096                     ) *                                                      &
     1097                                    ( u(k,j+3,i) + u(k,j-2,i) )               &
    10921098                               )
    10931099
    1094           diss_n(k) = - ABS ( v_comp ) * (                                   &
    1095                      ( 10.0 * ibit14 * adv_mom_5                             &
    1096                   +     3.0 * ibit13 * adv_mom_3                             &
    1097                   +           ibit12 * adv_mom_1                             &
    1098                      ) *                                                     &
    1099                                  ( u(k,j+1,i) - u(k,j,i)  )                  &
    1100               -      (  5.0 * ibit14 * adv_mom_5                             &
    1101                   +           ibit13 * adv_mom_3                             &
    1102                      ) *                                                     &
    1103                                  ( u(k,j+2,i) - u(k,j-1,i) )                 &
    1104               +      (        ibit14 * adv_mom_5                             &
    1105                      ) *                                                     &
    1106                                  ( u(k,j+3,i) - u(k,j-2,i) )                 &
     1100          diss_n(k) = - ABS ( v_comp ) * (                                    &
     1101                     ( 10.0_wp * ibit14 * adv_mom_5                           &
     1102                  +     3.0_wp * ibit13 * adv_mom_3                           &
     1103                  +              ibit12 * adv_mom_1                           &
     1104                     ) *                                                      &
     1105                                    ( u(k,j+1,i) - u(k,j,i)   )               &
     1106              -      (  5.0_wp * ibit14 * adv_mom_5                           &
     1107                  +              ibit13 * adv_mom_3                           &
     1108                     ) *                                                      &
     1109                                    ( u(k,j+2,i) - u(k,j-1,i) )               &
     1110              +      (           ibit14 * adv_mom_5                           &
     1111                     ) *                                                      &
     1112                                    ( u(k,j+3,i) - u(k,j-2,i) )               &
    11071113                                         )
    11081114!
    … …  
    11181124
    11191125          w_comp    = w(k,j,i) + w(k,j,i-1)
    1120           flux_t(k) = w_comp  * (                                            &
    1121                      ( 37.0 * ibit17 * adv_mom_5                             &
    1122                   +     7.0 * ibit16 * adv_mom_3                             &
    1123                   +           ibit15 * adv_mom_1                             &
    1124                      ) *                                                     &
    1125                              ( u(k+1,j,i)  + u(k,j,i)     )                  &
    1126               -      (  8.0 * ibit17 * adv_mom_5                             &
    1127                   +           ibit16 * adv_mom_3                             &
    1128                      ) *                                                     &
    1129                              ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
    1130               +      (        ibit17 * adv_mom_5                             &
    1131                      ) *                                                     &
    1132                              ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     1126          flux_t(k) = w_comp  * (                                             &
     1127                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     1128                  +     7.0_wp * ibit16 * adv_mom_3                           &
     1129                  +              ibit15 * adv_mom_1                           &
     1130                     ) *                                                      &
     1131                                ( u(k+1,j,i)  + u(k,j,i)     )                &
     1132              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     1133                  +              ibit16 * adv_mom_3                           &
     1134                     ) *                                                      &
     1135                                ( u(k_pp,j,i) + u(k-1,j,i)   )                &
     1136              +      (           ibit17 * adv_mom_5                           &
     1137                     ) *                                                      &
     1138                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
    11331139                                 )
    11341140
    1135           diss_t(k) = - ABS( w_comp ) * (                                    &
    1136                      ( 10.0 * ibit17 * adv_mom_5                             &
    1137                   +     3.0 * ibit16 * adv_mom_3                             &
    1138                   +           ibit15 * adv_mom_1                             &
    1139                      ) *                                                     &
    1140                              ( u(k+1,j,i)   - u(k,j,i)    )                  &
    1141               -      (  5.0 * ibit17 * adv_mom_5                             &
    1142                   +           ibit16 * adv_mom_3                             &
    1143                      ) *                                                     &
    1144                              ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
    1145               +      (        ibit17 * adv_mom_5                             &
    1146                      ) *                                                     &
    1147                              ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     1141          diss_t(k) = - ABS( w_comp ) * (                                     &
     1142                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     1143                  +     3.0_wp * ibit16 * adv_mom_3                           &
     1144                  +              ibit15 * adv_mom_1                           &
     1145                     ) *                                                      &
     1146                                ( u(k+1,j,i)   - u(k,j,i)    )                &
     1147              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     1148                  +              ibit16 * adv_mom_3                           &
     1149                     ) *                                                      &
     1150                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
     1151              +      (           ibit17 * adv_mom_5                           &
     1152                     ) *                                                      &
     1153                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
    11481154                                         )
    11491155!
    … …  
    11511157!--       correction is needed to overcome numerical instabilities introduced
    11521158!--       by a not sufficient reduction of divergences near topography.
    1153           div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx      &
    1154                +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy      &
    1155                +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)  &
    1156                 ) * 0.5
     1159          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
     1160               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
     1161               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)   &
     1162                ) * 0.5_wp
    11571163
    11581164          tend(k,j,i) = tend(k,j,i) - (                                       &
    … …  
    11761182           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
    11771183                              + ( flux_r(k) *                                 &
    1178                                 ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
     1184                                ( u_comp(k) - 2.0_wp * hom(k,1,1,0) )         &
    11791185                              / ( u_comp(k) - gu + 1.0E-20_wp   )             &
    11801186                              +   diss_r(k) *                                 &
    1181                                   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
     1187                                  ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) )    &
    11821188                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
    11831189                              *   weight_substep(intermediate_timestep_count)
    … …  
    11931199          u_comp(k) = u(k,j,i+1) + u(k,j,i)
    11941200          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
    1195                          37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
    1196                        -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
    1197                        +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
     1201                         37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                &
     1202                       -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
     1203                       +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    11981204             diss_r(k) = - ABS( u_comp(k) - gu ) * (                          &
    1199                          10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
    1200                        -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
    1201                        +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     1205                         10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
     1206                       -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
     1207                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    12021208
    12031209             v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    12041210             flux_n(k) = v_comp * (                                           &
    1205                          37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
    1206                        -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
    1207                        +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     1211                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
     1212                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
     1213                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    12081214             diss_n(k) = - ABS( v_comp ) * (                                  &
    1209                          10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
    1210                        -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
    1211                        +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     1215                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
     1216                       -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
     1217                       +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    12121218!
    12131219!--       k index has to be modified near bottom and top, else array
    … …  
    12221228
    12231229          w_comp    = w(k,j,i) + w(k,j,i-1)
    1224           flux_t(k) = w_comp  * (                                            &
    1225                      ( 37.0 * ibit17 * adv_mom_5                             &
    1226                   +     7.0 * ibit16 * adv_mom_3                             &
    1227                   +           ibit15 * adv_mom_1                             &
    1228                      ) *                                                     &
    1229                              ( u(k+1,j,i)  + u(k,j,i)     )                  &
    1230               -      (  8.0 * ibit17 * adv_mom_5                             &
    1231                   +           ibit16 * adv_mom_3                             &
    1232                      ) *                                                     &
    1233                              ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
    1234               +      (        ibit17 * adv_mom_5                             &
    1235                      ) *                                                     &
    1236                              ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
     1230          flux_t(k) = w_comp  * (                                             &
     1231                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     1232                  +     7.0_wp * ibit16 * adv_mom_3                           &
     1233                  +              ibit15 * adv_mom_1                           &
     1234                     ) *                                                      &
     1235                                ( u(k+1,j,i)  + u(k,j,i)     )                &
     1236              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     1237                  +              ibit16 * adv_mom_3                           &
     1238                     ) *                                                      &
     1239                                ( u(k_pp,j,i) + u(k-1,j,i)   )                &
     1240              +      (           ibit17 * adv_mom_5                           &
     1241                     ) *                                                      &
     1242                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
    12371243                                 )
    12381244
    1239           diss_t(k) = - ABS( w_comp ) * (                                    &
    1240                      ( 10.0 * ibit17 * adv_mom_5                             &
    1241                   +     3.0 * ibit16 * adv_mom_3                             &
    1242                   +           ibit15 * adv_mom_1                             &
    1243                      ) *                                                     &
    1244                              ( u(k+1,j,i)   - u(k,j,i)    )                  &
    1245               -      (  5.0 * ibit17 * adv_mom_5                             &
    1246                   +           ibit16 * adv_mom_3                             &
    1247                      ) *                                                     &
    1248                              ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
    1249               +      (        ibit17 * adv_mom_5                             &
    1250                      ) *                                                     &
    1251                              ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
     1245          diss_t(k) = - ABS( w_comp ) * (                                     &
     1246                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     1247                  +     3.0_wp * ibit16 * adv_mom_3                           &
     1248                  +              ibit15 * adv_mom_1                           &
     1249                     ) *                                                      &
     1250                                ( u(k+1,j,i)   - u(k,j,i)    )                &
     1251              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     1252                  +              ibit16 * adv_mom_3                           &
     1253                     ) *                                                      &
     1254                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
     1255              +      (           ibit17 * adv_mom_5                           &
     1256                     ) *                                                      &
     1257                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
    12521258                                         )
    12531259!
    … …  
    12551261!--       correction is needed to overcome numerical instabilities introduced
    12561262!--       by a not sufficient reduction of divergences near topography.
    1257           div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx      &
    1258                +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy      &
    1259                +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)  &
    1260                 ) * 0.5
     1263          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
     1264               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
     1265               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)   &
     1266                ) * 0.5_wp
    12611267
    12621268          tend(k,j,i) = tend(k,j,i) - (                                       &
    … …  
    12801286           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
    12811287                              + ( flux_r(k) *                                 &
    1282                                 ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
    1283                               / ( u_comp(k) - gu + 1.0E-20_wp   )             &
     1288                                ( u_comp(k) - 2.0_wp * hom(k,1,1,0)      )    &
     1289                              / ( u_comp(k) - gu + 1.0E-20_wp          )      &
    12841290                              +   diss_r(k) *                                 &
    1285                                   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
     1291                                  ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) )    &
    12861292                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
    12871293                              *   weight_substep(intermediate_timestep_count)
    … …  
    13631369       REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !:
    13641370
    1365        gu = 2.0 * u_gtrans
    1366        gv = 2.0 * v_gtrans
     1371       gu = 2.0_wp * u_gtrans
     1372       gv = 2.0_wp * v_gtrans
    13671373
    13681374!       
    … …  
    13771383
    13781384             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
    1379              flux_l_v(k,j,tn) = u_comp * (                                     &
    1380                               ( 37.0 * ibit20 * adv_mom_5                      &
    1381                            +     7.0 * ibit19 * adv_mom_3                      &
    1382                            +           ibit18 * adv_mom_1                      &
    1383                               ) *                                              &
    1384                                      ( v(k,j,i)   + v(k,j,i-1) )               &
    1385                        -      (  8.0 * ibit20 * adv_mom_5                      &
    1386                            +           ibit19 * adv_mom_3                      &
    1387                               ) *                                              &
    1388                                      ( v(k,j,i+1) + v(k,j,i-2) )               &
    1389                        +      (        ibit20 * adv_mom_5                      &
    1390                               ) *                                              &
    1391                                      ( v(k,j,i+2) + v(k,j,i-3) )               &
     1385             flux_l_v(k,j,tn) = u_comp * (                                    &
     1386                              ( 37.0_wp * ibit20 * adv_mom_5                  &
     1387                           +     7.0_wp * ibit19 * adv_mom_3                  &
     1388                           +              ibit18 * adv_mom_1                  &
     1389                              ) *                                             &
     1390                                        ( v(k,j,i)   + v(k,j,i-1) )           &
     1391                       -      (  8.0_wp * ibit20 * adv_mom_5                  &
     1392                           +              ibit19 * adv_mom_3                  &
     1393                              ) *                                             &
     1394                                        ( v(k,j,i+1) + v(k,j,i-2) )           &
     1395                       +      (           ibit20 * adv_mom_5                  &
     1396                              ) *                                             &
     1397                                        ( v(k,j,i+2) + v(k,j,i-3) )           &
    13921398                                         )
    13931399
    1394               diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
    1395                               ( 10.0 * ibit20 * adv_mom_5                      &
    1396                            +     3.0 * ibit19 * adv_mom_3                      &
    1397                            +           ibit18 * adv_mom_1                      &
    1398                               ) *                                              &
    1399                                      ( v(k,j,i)   - v(k,j,i-1) )               &
    1400                        -      (  5.0 * ibit20 * adv_mom_5                      &
    1401                            +           ibit19 * adv_mom_3                      &
    1402                               ) *                                              &
    1403                                      ( v(k,j,i+1) - v(k,j,i-2) )               &
    1404                        +      (        ibit20 * adv_mom_5                      &
    1405                               ) *                                              &
    1406                                      ( v(k,j,i+2) - v(k,j,i-3) )               &
     1400              diss_l_v(k,j,tn) = - ABS( u_comp ) * (                          &
     1401                              ( 10.0_wp * ibit20 * adv_mom_5                  &
     1402                           +     3.0_wp * ibit19 * adv_mom_3                  &
     1403                           +              ibit18 * adv_mom_1                  &
     1404                              ) *                                             &
     1405                                        ( v(k,j,i)   - v(k,j,i-1) )           &
     1406                       -      (  5.0_wp * ibit20 * adv_mom_5                  &
     1407                           +              ibit19 * adv_mom_3                  &
     1408                              ) *                                             &
     1409                                        ( v(k,j,i+1) - v(k,j,i-2) )           &
     1410                       +      (           ibit20 * adv_mom_5                  &
     1411                              ) *                                             &
     1412                                        ( v(k,j,i+2) - v(k,j,i-3) )           &
    14071413                                                   )
    14081414
    … …  
    14131419             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
    14141420             flux_l_v(k,j,tn) = u_comp * (                                    &
    1415                              37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
    1416                            -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
    1417                            +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
     1421                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )            &
     1422                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )            &
     1423                           +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    14181424             diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
    1419                              10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
    1420                            -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
    1421                            +        ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
     1425                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )            &
     1426                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )            &
     1427                           +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
    14221428
    14231429          ENDDO
    … …  
    14361442             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
    14371443             flux_s_v(k,tn) = v_comp_l * (                                    &
    1438                             ( 37.0 * ibit23 * adv_mom_5                       &
    1439                          +     7.0 * ibit22 * adv_mom_3                       &
    1440                          +           ibit21 * adv_mom_1                       &
     1444                            ( 37.0_wp * ibit23 * adv_mom_5                    &
     1445                         +     7.0_wp * ibit22 * adv_mom_3                    &
     1446                         +              ibit21 * adv_mom_1                    &
    14411447                            ) *                                               &
    1442                                      ( v(k,j,i)   + v(k,j-1,i) )              &
    1443                      -      (  8.0 * ibit23 * adv_mom_5                       &
    1444                          +           ibit22 * adv_mom_3                       &
     1448                                        ( v(k,j,i)   + v(k,j-1,i) )           &
     1449                     -      (  8.0_wp * ibit23 * adv_mom_5                    &
     1450                         +              ibit22 * adv_mom_3                    &
    14451451                            ) *                                               &
    1446                                      ( v(k,j+1,i) + v(k,j-2,i) )              &
    1447                      +      (        ibit23 * adv_mom_5                       &
     1452                                        ( v(k,j+1,i) + v(k,j-2,i) )           &
     1453                     +      (           ibit23 * adv_mom_5                    &
    14481454                            ) *                                               &
    1449                                      ( v(k,j+2,i) + v(k,j-3,i) )              &
     1455                                        ( v(k,j+2,i) + v(k,j-3,i) )           &
    14501456                                         )
    14511457
    14521458             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    1453                             ( 10.0 * ibit23 * adv_mom_5                       &
    1454                          +     3.0 * ibit22 * adv_mom_3                       &
    1455                          +           ibit21 * adv_mom_1                       &
     1459                            ( 10.0_wp * ibit23 * adv_mom_5                    &
     1460                         +     3.0_wp * ibit22 * adv_mom_3                    &
     1461                         +              ibit21 * adv_mom_1                    &
    14561462                            ) *                                               &
    1457                                      ( v(k,j,i)   - v(k,j-1,i) )              &
    1458                      -      (  5.0 * ibit23 * adv_mom_5                       &
    1459                          +           ibit22 * adv_mom_3                       &
     1463                                        ( v(k,j,i)   - v(k,j-1,i) )           &
     1464                     -      (  5.0_wp * ibit23 * adv_mom_5                    &
     1465                         +              ibit22 * adv_mom_3                    &
    14601466                            ) *                                               &
    1461                                      ( v(k,j+1,i) - v(k,j-2,i) )              &
    1462                      +      (        ibit23 * adv_mom_5                       &
     1467                                        ( v(k,j+1,i) - v(k,j-2,i) )           &
     1468                     +      (           ibit23 * adv_mom_5                    &
    14631469                            ) *                                               &
    1464                                      ( v(k,j+2,i) - v(k,j-3,i) )              &
    1465                                                  )
     1470                                        ( v(k,j+2,i) - v(k,j-3,i) )           &
     1471                                                  )
    14661472
    14671473          ENDDO
    … …  
    14711477             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
    14721478             flux_s_v(k,tn) = v_comp_l * (                                    &
    1473                            37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
    1474                          -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
    1475                          +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
     1479                           37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )              &
     1480                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )              &
     1481                         +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
    14761482             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    1477                            10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
    1478                          -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
    1479                          +        ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
     1483                           10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )              &
     1484                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )              &
     1485                         +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
    14801486
    14811487          ENDDO
    … …  
    14831489       ENDIF
    14841490
    1485        flux_t(0) = 0.0
    1486        diss_t(0) = 0.0
    1487        flux_d    = 0.0
    1488        diss_d    = 0.0
     1491       flux_t(0) = 0.0_wp
     1492       diss_t(0) = 0.0_wp
     1493       flux_d    = 0.0_wp
     1494       diss_d    = 0.0_wp
    14891495!
    14901496!--    Now compute the fluxes and tendency terms for the horizontal and
    … …  
    14971503 
    14981504          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1499           flux_r(k) = u_comp * (                                             &
    1500                      ( 37.0 * ibit20 * adv_mom_5                             &
    1501                   +     7.0 * ibit19 * adv_mom_3                             &
    1502                   +           ibit18 * adv_mom_1                             &
    1503                      ) *                                                     &
    1504                                  ( v(k,j,i+1) + v(k,j,i)   )                 &
    1505               -      (  8.0 * ibit20 * adv_mom_5                             &
    1506                   +           ibit19 * adv_mom_3                             &
    1507                      ) *                                                     &
    1508                                  ( v(k,j,i+2) + v(k,j,i-1) )                 &
    1509               +      (        ibit20 * adv_mom_5                             &
    1510                      ) *                                                     &
    1511                                  ( v(k,j,i+3) + v(k,j,i-2) )                 &
     1505          flux_r(k) = u_comp * (                                              &
     1506                     ( 37.0_wp * ibit20 * adv_mom_5                           &
     1507                  +     7.0_wp * ibit19 * adv_mom_3                           &
     1508                  +              ibit18 * adv_mom_1                           &
     1509                     ) *                                                      &
     1510                                    ( v(k,j,i+1) + v(k,j,i)   )               &
     1511              -      (  8.0_wp * ibit20 * adv_mom_5                           &
     1512                  +              ibit19 * adv_mom_3                           &
     1513                     ) *                                                      &
     1514                                    ( v(k,j,i+2) + v(k,j,i-1) )               &
     1515              +      (           ibit20 * adv_mom_5                           &
     1516                     ) *                                                      &
     1517                                    ( v(k,j,i+3) + v(k,j,i-2) )               &
    15121518                               )
    15131519
    1514           diss_r(k) = - ABS( u_comp ) * (                                    &
    1515                      ( 10.0 * ibit20 * adv_mom_5                             &
    1516                   +     3.0 * ibit19 * adv_mom_3                             &
    1517                   +           ibit18 * adv_mom_1                             &
    1518                      ) *                                                     &
    1519                                  ( v(k,j,i+1) - v(k,j,i)  )                  &
    1520               -      (  5.0 * ibit20 * adv_mom_5                             &
    1521                   +           ibit19 * adv_mom_3                             &
    1522                      ) *                                                     &
    1523                                  ( v(k,j,i+2) - v(k,j,i-1) )                 &
    1524               +      (        ibit20 * adv_mom_5                             &
    1525                      ) *                                                     &
    1526                                  ( v(k,j,i+3) - v(k,j,i-2) )                 &
     1520          diss_r(k) = - ABS( u_comp ) * (                                     &
     1521                     ( 10.0_wp * ibit20 * adv_mom_5                           &
     1522                  +     3.0_wp * ibit19 * adv_mom_3                           &
     1523                  +              ibit18 * adv_mom_1                           &
     1524                     ) *                                                      &
     1525                                    ( v(k,j,i+1) - v(k,j,i)  )                &
     1526              -      (  5.0_wp * ibit20 * adv_mom_5                           &
     1527                  +              ibit19 * adv_mom_3                           &
     1528                     ) *                                                      &
     1529                                    ( v(k,j,i+2) - v(k,j,i-1) )               &
     1530              +      (           ibit20 * adv_mom_5                           &
     1531                     ) *                                                      &
     1532                                    ( v(k,j,i+3) - v(k,j,i-2) )               &
    15271533                                        )
    15281534
    … …  
    15331539
    15341540          v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1535           flux_n(k) = ( v_comp(k) - gv ) * (                                 &
    1536                      ( 37.0 * ibit23 * adv_mom_5                             &
    1537                   +     7.0 * ibit22 * adv_mom_3                             &
    1538                   +           ibit21 * adv_mom_1                             &
    1539                      ) *                                                     &
    1540                                  ( v(k,j+1,i) + v(k,j,i)   )                 &
    1541               -      (  8.0 * ibit23 * adv_mom_5                             &
    1542                   +           ibit22 * adv_mom_3                             &
    1543                      ) *                                                     &
    1544                                  ( v(k,j+2,i) + v(k,j-1,i) )                 &
    1545               +      (        ibit23 * adv_mom_5                             &
    1546                      ) *                                                     &
    1547                                  ( v(k,j+3,i) + v(k,j-2,i) )                 &
    1548                                )
    1549 
    1550           diss_n(k) = - ABS( v_comp(k) - gv ) * (                            &
    1551                      ( 10.0 * ibit23 * adv_mom_5                             &
    1552                   +     3.0 * ibit22 * adv_mom_3                             &
    1553                   +           ibit21 * adv_mom_1                             &
    1554                      ) *                                                     &
    1555                                  ( v(k,j+1,i) - v(k,j,i)  )                  &
    1556               -      (  5.0 * ibit23 * adv_mom_5                             &
    1557                   +           ibit22 * adv_mom_3                             &
    1558                      ) *                                                     &
    1559                                  ( v(k,j+2,i) - v(k,j-1,i) )                 &
    1560               +      (        ibit23 * adv_mom_5                             &
    1561                      ) *                                                     &
    1562                                  ( v(k,j+3,i) - v(k,j-2,i) )                 &
    1563                                         )
     1541          flux_n(k) = ( v_comp(k) - gv ) * (                                  &
     1542                     ( 37.0_wp * ibit23 * adv_mom_5                           &
     1543                  +     7.0_wp * ibit22 * adv_mom_3                           &
     1544                  +              ibit21 * adv_mom_1                           &
     1545                     ) *                                                      &
     1546                                    ( v(k,j+1,i) + v(k,j,i)   )               &
     1547              -      (  8.0_wp * ibit23 * adv_mom_5                           &
     1548                  +              ibit22 * adv_mom_3                           &
     1549                     ) *                                                      &
     1550                                    ( v(k,j+2,i) + v(k,j-1,i) )               &
     1551              +      (           ibit23 * adv_mom_5                           &
     1552                     ) *                                                      &
     1553                                    ( v(k,j+3,i) + v(k,j-2,i) )               &
     1554                                           )
     1555
     1556          diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
     1557                     ( 10.0_wp * ibit23 * adv_mom_5                           &
     1558                  +     3.0_wp * ibit22 * adv_mom_3                           &
     1559                  +              ibit21 * adv_mom_1                           &
     1560                     ) *                                                      &
     1561                                    ( v(k,j+1,i) - v(k,j,i)   )               &
     1562              -      (  5.0_wp * ibit23 * adv_mom_5                           &
     1563                  +              ibit22 * adv_mom_3                           &
     1564                     ) *                                                      &
     1565                                    ( v(k,j+2,i) - v(k,j-1,i) )               &
     1566              +      (           ibit23 * adv_mom_5                           &
     1567                     ) *                                                      &
     1568                                    ( v(k,j+3,i) - v(k,j-2,i) )               &
     1569                                                )
    15641570!
    15651571!--       k index has to be modified near bottom and top, else array
    … …  
    15741580
    15751581          w_comp    = w(k,j-1,i) + w(k,j,i)
    1576           flux_t(k) = w_comp  * (                                            &
    1577                      ( 37.0 * ibit26 * adv_mom_5                             &
    1578                   +     7.0 * ibit25 * adv_mom_3                             &
    1579                   +           ibit24 * adv_mom_1                             &
    1580                      ) *                                                     &
    1581                              ( v(k+1,j,i)   + v(k,j,i)    )                  &
    1582               -      (  8.0 * ibit26 * adv_mom_5                             &
    1583                   +           ibit25 * adv_mom_3                             &
    1584                      ) *                                                     &
    1585                              ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
    1586               +      (        ibit26 * adv_mom_5                             &
    1587                      ) *                                                     &
    1588                              ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
     1582          flux_t(k) = w_comp  * (                                             &
     1583                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     1584                  +     7.0_wp * ibit25 * adv_mom_3                           &
     1585                  +              ibit24 * adv_mom_1                           &
     1586                     ) *                                                      &
     1587                                ( v(k+1,j,i)   + v(k,j,i)    )                &
     1588              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     1589                  +              ibit25 * adv_mom_3                           &
     1590                     ) *                                                      &
     1591                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
     1592              +      (           ibit26 * adv_mom_5                           &
     1593                     ) *                                                      &
     1594                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
    15891595                                 )
    15901596
    1591           diss_t(k) = - ABS( w_comp ) * (                                    &
    1592                      ( 10.0 * ibit26 * adv_mom_5                             &
    1593                   +     3.0 * ibit25 * adv_mom_3                             &
    1594                   +           ibit24 * adv_mom_1                             &
    1595                      ) *                                                     &
    1596                              ( v(k+1,j,i)   - v(k,j,i)    )                  &
    1597               -      (  5.0 * ibit26 * adv_mom_5                             &
    1598                   +           ibit25 * adv_mom_3                             &
    1599                      ) *                                                     &
    1600                              ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
    1601               +      (        ibit26 * adv_mom_5                             &
    1602                      ) *                                                     &
    1603                              ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
     1597          diss_t(k) = - ABS( w_comp ) * (                                     &
     1598                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     1599                  +     3.0_wp * ibit25 * adv_mom_3                           &
     1600                  +              ibit24 * adv_mom_1                           &
     1601                     ) *                                                      &
     1602                                ( v(k+1,j,i)   - v(k,j,i)    )                &
     1603              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     1604                  +              ibit25 * adv_mom_3                           &
     1605                     ) *                                                      &
     1606                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
     1607              +      (           ibit26 * adv_mom_5                           &
     1608                     ) *                                                      &
     1609                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
    16041610                                         )
    16051611!
    … …  
    16071613!--       correction is needed to overcome numerical instabilities introduced
    16081614!--       by a not sufficient reduction of divergences near topography.
    1609           div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx      &
    1610                +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy      &
    1611                +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)  &
    1612                 ) * 0.5
     1615          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
     1616               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
     1617               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)   &
     1618                ) * 0.5_wp
    16131619
    16141620          tend(k,j,i) = tend(k,j,i) - (                                       &
    … …  
    16331639           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
    16341640             + ( flux_n(k)                                                    &
    1635              * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
    1636              / ( v_comp(k) - gv + 1.0E-20_wp )                                &
     1641             * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)      )                     &
     1642             / ( v_comp(k) - gv + 1.0E-20_wp       )                          &
    16371643             +   diss_n(k)                                                    &
    1638              *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
     1644             *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) )                     &
    16391645             / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )                        &
    16401646             *   weight_substep(intermediate_timestep_count)
    … …  
    16501656
    16511657          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    1652           flux_r(k) = u_comp * (                                           &
    1653                       37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
    1654                     -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                   &
    1655                     +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    1656 
    1657           diss_r(k) = - ABS( u_comp ) * (                                  &
    1658                       10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
    1659                     -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
    1660                     +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
     1658          flux_r(k) = u_comp * (                                              &
     1659                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
     1660                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
     1661                    +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
     1662
     1663          diss_r(k) = - ABS( u_comp ) * (                                     &
     1664                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
     1665                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
     1666                    +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    16611667
    16621668
    16631669          v_comp(k) = v(k,j+1,i) + v(k,j,i)
    1664           flux_n(k) = ( v_comp(k) - gv ) * (                               &
    1665                       37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
    1666                     -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                   &
    1667                       +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    1668 
    1669           diss_n(k) = - ABS( v_comp(k) - gv ) * (                          &
    1670                       10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
    1671                     -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
    1672                     +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
     1670          flux_n(k) = ( v_comp(k) - gv ) * (                                  &
     1671                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                   &
     1672                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                   &
     1673                      +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
     1674
     1675          diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
     1676                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                   &
     1677                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                   &
     1678                    +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    16731679!
    16741680!--       k index has to be modified near bottom and top, else array
    … …  
    16831689
    16841690          w_comp    = w(k,j-1,i) + w(k,j,i)
    1685           flux_t(k) = w_comp  * (                                            &
    1686                      ( 37.0 * ibit26 * adv_mom_5                             &
    1687                   +     7.0 * ibit25 * adv_mom_3                             &
    1688                   +           ibit24 * adv_mom_1                             &
    1689                      ) *                                                     &
    1690                              ( v(k+1,j,i)   + v(k,j,i)    )                  &
    1691               -      (  8.0 * ibit26 * adv_mom_5                             &
    1692                   +           ibit25 * adv_mom_3                             &
    1693                      ) *                                                     &
    1694                              ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
    1695               +      (        ibit26 * adv_mom_5                             &
    1696                      ) *                                                     &
    1697                              ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
     1691          flux_t(k) = w_comp  * (                                             &
     1692                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     1693                  +     7.0_wp * ibit25 * adv_mom_3                           &
     1694                  +              ibit24 * adv_mom_1                           &
     1695                     ) *                                                      &
     1696                                ( v(k+1,j,i)   + v(k,j,i)    )                &
     1697              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     1698                  +              ibit25 * adv_mom_3                           &
     1699                     ) *                                                      &
     1700                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
     1701              +      (           ibit26 * adv_mom_5                           &
     1702                     ) *                                                      &
     1703                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
    16981704                                 )
    16991705
    1700           diss_t(k) = - ABS( w_comp ) * (                                    &
    1701                      ( 10.0 * ibit26 * adv_mom_5                             &
    1702                   +     3.0 * ibit25 * adv_mom_3                             &
    1703                   +           ibit24 * adv_mom_1                             &
    1704                      ) *                                                     &
    1705                              ( v(k+1,j,i)   - v(k,j,i)    )                  &
    1706               -      (  5.0 * ibit26 * adv_mom_5                             &
    1707                   +           ibit25 * adv_mom_3                             &
    1708                      ) *                                                     &
    1709                              ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
    1710               +      (        ibit26 * adv_mom_5                             &
    1711                      ) *                                                     &
    1712                              ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
     1706          diss_t(k) = - ABS( w_comp ) * (                                     &
     1707                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     1708                  +     3.0_wp * ibit25 * adv_mom_3                           &
     1709                  +              ibit24 * adv_mom_1                           &
     1710                     ) *                                                      &
     1711                                ( v(k+1,j,i)   - v(k,j,i)    )                &
     1712              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     1713                  +              ibit25 * adv_mom_3                           &
     1714                     ) *                                                      &
     1715                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
     1716              +      (           ibit26 * adv_mom_5                           &
     1717                     ) *                                                      &
     1718                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
    17131719                                         )
    17141720!
    … …  
    17161722!--       correction is needed to overcome numerical instabilities introduced
    17171723!--       by a not sufficient reduction of divergences near topography.
    1718           div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx      &
    1719                +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy      &
    1720                +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)  &
    1721                 ) * 0.5
     1724          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
     1725               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
     1726               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)   &
     1727                ) * 0.5_wp
    17221728
    17231729          tend(k,j,i) = tend(k,j,i) - (                                       &
    … …  
    17421748           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
    17431749             + ( flux_n(k)                                                    &
    1744              * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
    1745              / ( v_comp(k) - gv + 1.0E-20_wp )                                &
     1750             * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)      )                     &
     1751             / ( v_comp(k) - gv + 1.0E-20_wp       )                          &
    17461752             +   diss_n(k)                                                    &
    1747              *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
     1753             *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) )                     &
    17481754             / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )                        &
    17491755             *   weight_substep(intermediate_timestep_count)
    17501756!
    17511757!--        Statistical Evaluation of w'v'.
    1752            sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
    1753                               + ( flux_t(k) + diss_t(k) )                      &
     1758           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                        &
     1759                              + ( flux_t(k) + diss_t(k) )                     &
    17541760                              *   weight_substep(intermediate_timestep_count)
    17551761
    … …  
    17671773    SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
    17681774
    1769        USE arrays_3d,                                                          &
     1775       USE arrays_3d,                                                         &
    17701776           ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w
    17711777
    1772        USE constants,                                                          &
     1778       USE constants,                                                         &
    17731779           ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
    17741780
    1775        USE control_parameters,                                                 &
     1781       USE control_parameters,                                                &
    17761782           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    17771783
    1778        USE grid_variables,                                                     &
     1784       USE grid_variables,                                                    &
    17791785           ONLY:  ddx, ddy
    17801786
    1781        USE indices,                                                            &
    1782            ONLY:  nxl, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0,         &
     1787       USE indices,                                                           &
     1788           ONLY:  nxl, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0,        &
    17831789                  wall_flags_00
    17841790
    17851791       USE kinds
    17861792       
    1787        USE statistics,                                                         &
     1793       USE statistics,                                                        &
    17881794           ONLY:  hom, sums_ws2_ws_l, weight_substep
    17891795
    … …  
    18241830       REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_t !:
    18251831
    1826        gu = 2.0 * u_gtrans
    1827        gv = 2.0 * v_gtrans
     1832       gu = 2.0_wp * u_gtrans
     1833       gv = 2.0_wp * v_gtrans
    18281834
    18291835!
    … …  
    18381844             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
    18391845             flux_s_w(k,tn) = v_comp * (                                      &
    1840                             ( 37.0 * ibit32 * adv_mom_5                       &
    1841                          +     7.0 * ibit31 * adv_mom_3                       &
    1842                          +           ibit30 * adv_mom_1                       &
     1846                            ( 37.0_wp * ibit32 * adv_mom_5                    &
     1847                         +     7.0_wp * ibit31 * adv_mom_3                    &
     1848                         +              ibit30 * adv_mom_1                    &
    18431849                            ) *                                               &
    1844                                      ( w(k,j,i)   + w(k,j-1,i) )              &
    1845                      -      (  8.0 * ibit32 * adv_mom_5                       &
    1846                          +           ibit31 * adv_mom_3                       &
     1850                                        ( w(k,j,i)   + w(k,j-1,i) )           &
     1851                     -      (  8.0_wp * ibit32 * adv_mom_5                    &
     1852                         +              ibit31 * adv_mom_3                    &
    18471853                            ) *                                               &
    1848                                      ( w(k,j+1,i) + w(k,j-2,i) )              &
    1849                      +      (        ibit32 * adv_mom_5                       &
     1854                                        ( w(k,j+1,i) + w(k,j-2,i) )           &
     1855                     +      (           ibit32 * adv_mom_5                    &
    18501856                            ) *                                               &
    1851                                      ( w(k,j+2,i) + w(k,j-3,i) )              &
    1852                                          )
     1857                                        ( w(k,j+2,i) + w(k,j-3,i) )           &
     1858                                       )
    18531859
    18541860             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
    1855                             ( 10.0 * ibit32 * adv_mom_5                       &
    1856                          +     3.0 * ibit31 * adv_mom_3                       &
    1857                          +           ibit30 * adv_mom_1                       &
     1861                            ( 10.0_wp * ibit32 * adv_mom_5                    &
     1862                         +     3.0_wp * ibit31 * adv_mom_3                    &
     1863                         +              ibit30 * adv_mom_1                    &
    18581864                            ) *                                               &
    1859                                      ( w(k,j,i)   - w(k,j-1,i) )              &
    1860                      -      (  5.0 * ibit32 * adv_mom_5                       &
    1861                          +           ibit31 * adv_mom_3                       &
     1865                                        ( w(k,j,i)   - w(k,j-1,i) )           &
     1866                     -      (  5.0_wp * ibit32 * adv_mom_5                    &
     1867                         +              ibit31 * adv_mom_3                    &
    18621868                            ) *                                               &
    1863                                      ( w(k,j+1,i) - w(k,j-2,i) )              &
    1864                      +      (        ibit32 * adv_mom_5                       &
     1869                                        ( w(k,j+1,i) - w(k,j-2,i) )           &
     1870                     +      (           ibit32 * adv_mom_5                    &
    18651871                            ) *                                               &
    1866                                      ( w(k,j+2,i) - w(k,j-3,i) )              &
    1867                                                  )
     1872                                        ( w(k,j+2,i) - w(k,j-3,i) )           &
     1873                                                )
    18681874
    18691875          ENDDO
    … …  
    18731879             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
    18741880             flux_s_w(k,tn) = v_comp * (                                      &
    1875                            37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
    1876                          -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
    1877                          +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
     1881                           37.0_wp * ( w(k,j,i) + w(k,j-1,i)   )              &
     1882                         -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )              &
     1883                         +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    18781884             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
    1879                            10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
    1880                          -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
    1881                          +        ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
     1885                           10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )              &
     1886                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )              &
     1887                         +           ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
    18821888
    18831889          ENDDO
    … …  
    18951901
    18961902             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
    1897              flux_l_w(k,j,tn) = u_comp * (                                     &
    1898                              ( 37.0 * ibit29 * adv_mom_5                       &
    1899                           +     7.0 * ibit28 * adv_mom_3                       &
    1900                           +           ibit27 * adv_mom_1                       &
    1901                              ) *                                               &
    1902                                      ( w(k,j,i)   + w(k,j,i-1) )               &
    1903                       -      (  8.0 * ibit29 * adv_mom_5                       &
    1904                           +           ibit28 * adv_mom_3                       &
    1905                              ) *                                               &
    1906                                      ( w(k,j,i+1) + w(k,j,i-2) )               &
    1907                       +      (        ibit29 * adv_mom_5                       &
    1908                              ) *                                               &
    1909                                      ( w(k,j,i+2) + w(k,j,i-3) )               &
     1903             flux_l_w(k,j,tn) = u_comp * (                                    &
     1904                             ( 37.0_wp * ibit29 * adv_mom_5                   &
     1905                          +     7.0_wp * ibit28 * adv_mom_3                   &
     1906                          +              ibit27 * adv_mom_1                   &
     1907                             ) *                                              &
     1908                                        ( w(k,j,i)   + w(k,j,i-1) )           &
     1909                      -      (  8.0_wp * ibit29 * adv_mom_5                   &
     1910                          +              ibit28 * adv_mom_3                   &
     1911                             ) *                                              &
     1912                                        ( w(k,j,i+1) + w(k,j,i-2) )           &
     1913                      +      (           ibit29 * adv_mom_5                   &
     1914                             ) *                                              &
     1915                                        ( w(k,j,i+2) + w(k,j,i-3) )           &
    19101916                                         )
    19111917
    1912                diss_l_w(k,j,tn) = - ABS( u_comp ) * (                          &
    1913                              ( 10.0 * ibit29 * adv_mom_5                       &
    1914                           +     3.0 * ibit28 * adv_mom_3                       &
    1915                           +           ibit27 * adv_mom_1                       &
    1916                              ) *                                               &
    1917                                      ( w(k,j,i)   - w(k,j,i-1) )               &
    1918                       -      (  5.0 * ibit29 * adv_mom_5                       &
    1919                           +           ibit28 * adv_mom_3                       &
    1920                              ) *                                               &
    1921                                      ( w(k,j,i+1) - w(k,j,i-2) )               &
    1922                       +      (        ibit29 * adv_mom_5                       &
    1923                              ) *                                               &
    1924                                      ( w(k,j,i+2) - w(k,j,i-3) )               &
     1918               diss_l_w(k,j,tn) = - ABS( u_comp ) * (                         &
     1919                             ( 10.0_wp * ibit29 * adv_mom_5                   &
     1920                          +     3.0_wp * ibit28 * adv_mom_3                   &
     1921                          +              ibit27 * adv_mom_1                   &
     1922                             ) *                                              &
     1923                                        ( w(k,j,i)   - w(k,j,i-1) )           &
     1924                      -      (  5.0_wp * ibit29 * adv_mom_5                   &
     1925                          +              ibit28 * adv_mom_3                   &
     1926                             ) *                                              &
     1927                                        ( w(k,j,i+1) - w(k,j,i-2) )           &
     1928                      +      (           ibit29 * adv_mom_5                   &
     1929                             ) *                                              &
     1930                                        ( w(k,j,i+2) - w(k,j,i-3) )           &
    19251931                                                    )
    19261932
    … …  
    19311937             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
    19321938             flux_l_w(k,j,tn) = u_comp * (                                    &
    1933                             37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
    1934                           -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
    1935                           +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
     1939                            37.0_wp * ( w(k,j,i) + w(k,j,i-1)   )             &
     1940                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )             &
     1941                          +           ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    19361942             diss_l_w(k,j,tn) = - ABS( u_comp ) * (                           &
    1937                             10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
    1938                           -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
    1939                           +        ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
     1943                            10.0_wp * ( w(k,j,i) - w(k,j,i-1)   )             &
     1944                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )             &
     1945                          +           ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
    19401946
    19411947          ENDDO
    … …  
    19621968
    19631969          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    1964           flux_r(k) = u_comp * (                                             &
    1965                      ( 37.0 * ibit29 * adv_mom_5                             &
    1966                   +     7.0 * ibit28 * adv_mom_3                             &
    1967                   +           ibit27 * adv_mom_1                             &
    1968                      ) *                                                     &
    1969                                  ( w(k,j,i+1) + w(k,j,i)   )                 &
    1970               -      (  8.0 * ibit29 * adv_mom_5                             &
    1971                   +           ibit28 * adv_mom_3                             &
    1972                      ) *                                                     &
    1973                                  ( w(k,j,i+2) + w(k,j,i-1) )                 &
    1974               +      (        ibit29 * adv_mom_5                             &
    1975                      ) *                                                     &
    1976                                  ( w(k,j,i+3) + w(k,j,i-2) )                 &
    1977                                 )
    1978 
    1979           diss_r(k) = - ABS( u_comp ) * (                                    &
    1980                      ( 10.0 * ibit29 * adv_mom_5                             &
    1981                   +     3.0 * ibit28 * adv_mom_3                             &
    1982                   +           ibit27 * adv_mom_1                             &
    1983                      ) *                                                     &
    1984                                  ( w(k,j,i+1) - w(k,j,i)  )                  &
    1985               -      (  5.0 * ibit29 * adv_mom_5                             &
    1986                   +           ibit28 * adv_mom_3                             &
    1987                      ) *                                                     &
    1988                                  ( w(k,j,i+2) - w(k,j,i-1) )                 &
    1989               +      (        ibit29 * adv_mom_5                             &
    1990                      ) *                                                     &
    1991                                  ( w(k,j,i+3) - w(k,j,i-2) )                 &
     1970          flux_r(k) = u_comp * (                                              &
     1971                     ( 37.0_wp * ibit29 * adv_mom_5                           &
     1972                  +     7.0_wp * ibit28 * adv_mom_3                           &
     1973                  +              ibit27 * adv_mom_1                           &
     1974                     ) *                                                      &
     1975                                    ( w(k,j,i+1) + w(k,j,i)   )               &
     1976              -      (  8.0_wp * ibit29 * adv_mom_5                           &
     1977                  +              ibit28 * adv_mom_3                           &
     1978                     ) *                                                      &
     1979                                    ( w(k,j,i+2) + w(k,j,i-1) )               &
     1980              +      (           ibit29 * adv_mom_5                           &
     1981                     ) *                                                      &
     1982                                    ( w(k,j,i+3) + w(k,j,i-2) )               &
     1983                               )
     1984
     1985          diss_r(k) = - ABS( u_comp ) * (                                     &
     1986                     ( 10.0_wp * ibit29 * adv_mom_5                           &
     1987                  +     3.0_wp * ibit28 * adv_mom_3                           &
     1988                  +              ibit27 * adv_mom_1                           &
     1989                     ) *                                                      &
     1990                                    ( w(k,j,i+1) - w(k,j,i)   )               &
     1991              -      (  5.0_wp * ibit29 * adv_mom_5                           &
     1992                  +              ibit28 * adv_mom_3                           &
     1993                     ) *                                                      &
     1994                                    ( w(k,j,i+2) - w(k,j,i-1) )               &
     1995              +      (           ibit29 * adv_mom_5                           &
     1996                     ) *                                                      &
     1997                                    ( w(k,j,i+3) - w(k,j,i-2) )               &
    19921998                                        )
    19931999
    … …  
    19972003
    19982004          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    1999           flux_n(k) = v_comp * (                                             &
    2000                      ( 37.0 * ibit32 * adv_mom_5                             &
    2001                   +     7.0 * ibit31 * adv_mom_3                             &
    2002                   +           ibit30 * adv_mom_1                             &
    2003                      ) *                                                     &
    2004                                  ( w(k,j+1,i) + w(k,j,i)   )                 &
    2005               -      (  8.0 * ibit32 * adv_mom_5                             &
    2006                   +           ibit31 * adv_mom_3                             &
    2007                      ) *                                                     &
    2008                                  ( w(k,j+2,i) + w(k,j-1,i) )                 &
    2009               +      (        ibit32 * adv_mom_5                             &
    2010                      ) *                                                     &
    2011                                  ( w(k,j+3,i) + w(k,j-2,i) )                 &
     2005          flux_n(k) = v_comp * (                                              &
     2006                     ( 37.0_wp * ibit32 * adv_mom_5                           &
     2007                  +     7.0_wp * ibit31 * adv_mom_3                           &
     2008                  +              ibit30 * adv_mom_1                           &
     2009                     ) *                                                      &
     2010                                    ( w(k,j+1,i) + w(k,j,i)   )               &
     2011              -      (  8.0_wp * ibit32 * adv_mom_5                           &
     2012                  +              ibit31 * adv_mom_3                           &
     2013                     ) *                                                      &
     2014                                    ( w(k,j+2,i) + w(k,j-1,i) )               &
     2015              +      (           ibit32 * adv_mom_5                           &
     2016                     ) *                                                      &
     2017                                    ( w(k,j+3,i) + w(k,j-2,i) )               &
    20122018                               )
    20132019
    2014           diss_n(k) = - ABS( v_comp ) * (                                    &
    2015                      ( 10.0 * ibit32 * adv_mom_5                             &
    2016                   +     3.0 * ibit31 * adv_mom_3                             &
    2017                   +           ibit30 * adv_mom_1                             &
    2018                      ) *                                                     &
    2019                                  ( w(k,j+1,i) - w(k,j,i)  )                  &
    2020               -      (  5.0 * ibit32 * adv_mom_5                             &
    2021                   +           ibit31 * adv_mom_3                             &
    2022                      ) *                                                     &
    2023                                  ( w(k,j+2,i) - w(k,j-1,i) )                 &
    2024               +      (        ibit32 * adv_mom_5                             &
    2025                      ) *                                                     &
    2026                                  ( w(k,j+3,i) - w(k,j-2,i) )                 &
     2020          diss_n(k) = - ABS( v_comp ) * (                                     &
     2021                     ( 10.0_wp * ibit32 * adv_mom_5                           &
     2022                  +     3.0_wp * ibit31 * adv_mom_3                           &
     2023                  +              ibit30 * adv_mom_1                           &
     2024                     ) *                                                      &
     2025                                    ( w(k,j+1,i) - w(k,j,i)  )                &
     2026              -      (  5.0_wp * ibit32 * adv_mom_5                           &
     2027                  +              ibit31 * adv_mom_3                           &
     2028                     ) *                                                      &
     2029                                   ( w(k,j+2,i) - w(k,j-1,i) )                &
     2030              +      (           ibit32 * adv_mom_5                           &
     2031                     ) *                                                      &
     2032                                   ( w(k,j+3,i) - w(k,j-2,i) )                &
    20272033                                        )
    20282034!
    … …  
    20382044
    20392045          w_comp    = w(k+1,j,i) + w(k,j,i)
    2040           flux_t(k) = w_comp  * (                                            &
    2041                      ( 37.0 * ibit35 * adv_mom_5                             &
    2042                   +     7.0 * ibit34 * adv_mom_3                             &
    2043                   +           ibit33 * adv_mom_1                             &
    2044                      ) *                                                     &
    2045                              ( w(k+1,j,i)  + w(k,j,i)     )                  &
    2046               -      (  8.0 * ibit35 * adv_mom_5                             &
    2047                   +           ibit34 * adv_mom_3                             &
    2048                      ) *                                                     &
    2049                              ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
    2050               +      (        ibit35 * adv_mom_5                             &
    2051                      ) *                                                     &
    2052                              ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
    2053                                  )
    2054 
    2055           diss_t(k) = - ABS( w_comp ) * (                                    &
    2056                      ( 10.0 * ibit35 * adv_mom_5                             &
    2057                   +     3.0 * ibit34 * adv_mom_3                             &
    2058                   +           ibit33 * adv_mom_1                             &
    2059                      ) *                                                     &
    2060                              ( w(k+1,j,i)   - w(k,j,i)    )                  &
    2061               -      (  5.0 * ibit35 * adv_mom_5                             &
    2062                   +           ibit34 * adv_mom_3                             &
    2063                      ) *                                                     &
    2064                              ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
    2065               +      (        ibit35 * adv_mom_5                             &
    2066                      ) *                                                     &
    2067                              ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
    2068                                          )
     2046          flux_t(k) = w_comp  * (                                             &
     2047                     ( 37.0_wp * ibit35 * adv_mom_5                           &
     2048                  +     7.0_wp * ibit34 * adv_mom_3                           &
     2049                  +              ibit33 * adv_mom_1                           &
     2050                     ) *                                                      &
     2051                                ( w(k+1,j,i)  + w(k,j,i)     )                &
     2052              -      (  8.0_wp * ibit35 * adv_mom_5                           &
     2053                  +              ibit34 * adv_mom_3                           &
     2054                     ) *                                                      &
     2055                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
     2056              +      (           ibit35 * adv_mom_5                           &
     2057                     ) *                                                      &
     2058                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
     2059                                )
     2060
     2061          diss_t(k) = - ABS( w_comp ) * (                                     &
     2062                     ( 10.0_wp * ibit35 * adv_mom_5                           &
     2063                  +     3.0_wp * ibit34 * adv_mom_3                           &
     2064                  +              ibit33 * adv_mom_1                           &
     2065                     ) *                                                      &
     2066                                ( w(k+1,j,i)   - w(k,j,i)    )                &
     2067              -      (  5.0_wp * ibit35 * adv_mom_5                           &
     2068                  +              ibit34 * adv_mom_3                           &
     2069                     ) *                                                      &
     2070                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
     2071              +      (           ibit35 * adv_mom_5                           &
     2072                     ) *                                                      &
     2073                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
     2074                                        )
    20692075
    20702076!
    … …  
    20752081              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
    20762082              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
    2077                 ) * 0.5
     2083                ) * 0.5_wp
    20782084
    20792085          tend(k,j,i) = tend(k,j,i) - (                                       &
    … …  
    20942100!
    20952101!--        Statistical Evaluation of w'w'.
    2096           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
    2097                              + ( flux_t(k) + diss_t(k) )                     &
     2102          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
     2103                             + ( flux_t(k) + diss_t(k) )                      &
    20982104                             *   weight_substep(intermediate_timestep_count)
    20992105
    … …  
    21032109
    21042110          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    2105           flux_r(k) = u_comp * (                                            &
    2106                       37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
    2107                     -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
    2108                     +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    2109 
    2110           diss_r(k) = - ABS( u_comp ) * (                                   &
    2111                       10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
    2112                     -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
    2113                     +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
     2111          flux_r(k) = u_comp * (                                              &
     2112                      37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
     2113                    -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
     2114                    +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
     2115
     2116          diss_r(k) = - ABS( u_comp ) * (                                     &
     2117                      10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
     2118                    -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
     2119                    +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    21142120
    21152121          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    2116           flux_n(k) = v_comp * (                                            &
    2117                       37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
    2118                     -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
    2119                     +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    2120 
    2121           diss_n(k) = - ABS( v_comp ) * (                                   &
    2122                       10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
    2123                     -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
    2124                     +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
     2122          flux_n(k) = v_comp * (                                              &
     2123                      37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
     2124                    -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
     2125                    +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
     2126
     2127          diss_n(k) = - ABS( v_comp ) * (                                     &
     2128                      10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
     2129                    -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
     2130                    +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    21252131!
    21262132!--       k index has to be modified near bottom and top, else array
    … …  
    21352141
    21362142          w_comp    = w(k+1,j,i) + w(k,j,i)
    2137           flux_t(k) = w_comp  * (                                            &
    2138                      ( 37.0 * ibit35 * adv_mom_5                             &
    2139                   +     7.0 * ibit34 * adv_mom_3                             &
    2140                   +           ibit33 * adv_mom_1                             &
    2141                      ) *                                                     &
    2142                              ( w(k+1,j,i)  + w(k,j,i)     )                  &
    2143               -      (  8.0 * ibit35 * adv_mom_5                             &
    2144                   +           ibit34 * adv_mom_3                             &
    2145                      ) *                                                     &
    2146                              ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
    2147               +      (        ibit35 * adv_mom_5                             &
    2148                      ) *                                                     &
    2149                              ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
    2150                                  )
    2151 
    2152           diss_t(k) = - ABS( w_comp ) * (                                    &
    2153                      ( 10.0 * ibit35 * adv_mom_5                             &
    2154                   +     3.0 * ibit34 * adv_mom_3                             &
    2155                   +           ibit33 * adv_mom_1                             &
    2156                      ) *                                                     &
    2157                              ( w(k+1,j,i)   - w(k,j,i)    )                  &
    2158               -      (  5.0 * ibit35 * adv_mom_5                             &
    2159                   +           ibit34 * adv_mom_3                             &
    2160                      ) *                                                     &
    2161                              ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
    2162               +      (        ibit35 * adv_mom_5                             &
    2163                      ) *                                                     &
    2164                              ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
    2165                                          )
     2143          flux_t(k) = w_comp  * (                                             &
     2144                     ( 37.0_wp * ibit35 * adv_mom_5                           &
     2145                  +     7.0_wp * ibit34 * adv_mom_3                           &
     2146                  +              ibit33 * adv_mom_1                           &
     2147                     ) *                                                      &
     2148                                ( w(k+1,j,i)  + w(k,j,i)     )                &
     2149              -      (  8.0_wp * ibit35 * adv_mom_5                           &
     2150                  +              ibit34 * adv_mom_3                           &
     2151                     ) *                                                      &
     2152                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
     2153              +      (           ibit35 * adv_mom_5                           &
     2154                     ) *                                                      &
     2155                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
     2156                                )
     2157
     2158          diss_t(k) = - ABS( w_comp ) * (                                     &
     2159                     ( 10.0_wp * ibit35 * adv_mom_5                           &
     2160                  +     3.0_wp * ibit34 * adv_mom_3                           &
     2161                  +              ibit33 * adv_mom_1                           &
     2162                     ) *                                                      &
     2163                                ( w(k+1,j,i)   - w(k,j,i)    )                &
     2164              -      (  5.0_wp * ibit35 * adv_mom_5                           &
     2165                  +              ibit34 * adv_mom_3                           &
     2166                     ) *                                                      &
     2167                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
     2168              +      (           ibit35 * adv_mom_5                           &
     2169                     ) *                                                      &
     2170                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
     2171                                        )
    21662172!
    21672173!--       Calculate the divergence of the velocity field. A respective
    … …  
    21712177              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
    21722178              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &