Changeset 2731


Ignore:
Timestamp:
Jan 9, 2018 5:44:02 PM (7 years ago)
Author:
suehring
Message:

Enable vectorization by splitting the k-loop

File:
1 edited

Legend:

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

    r2718 r2731  
    2525! -----------------
    2626! $Id$
     27! Enable loop vectorization by splitting the k-loop
     28!
     29! 2718 2018-01-02 08:49:38Z maronga
    2730! Corrected "Former revisions" section
    2831!
     
    13531356       flux_t(0) = 0.0_wp
    13541357       diss_t(0) = 0.0_wp
    1355        flux_d    = 0.0_wp
    1356        diss_d    = 0.0_wp
    13571358!       
    13581359!--    Now compute the fluxes and tendency terms for the horizontal and
     
    14731474                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    14741475                                         )
     1476       ENDDO
     1477
     1478       DO  k = nzb+1, nzb_max
     1479
     1480          flux_d    = flux_t(k-1)
     1481          diss_d    = diss_t(k-1)
    14751482!
    14761483!--       Calculate the divergence of the velocity field. A respective
     
    15131520          swap_flux_x_local(k,j,tn) = flux_r(k)
    15141521          swap_diss_x_local(k,j,tn) = diss_r(k)
    1515           flux_d                    = flux_t(k)
    1516           diss_d                    = diss_t(k)
    15171522
    15181523       ENDDO
     
    15811586                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    15821587                                         )
     1588
     1589       ENDDO
     1590
     1591       DO  k = nzb_max+1, nzt
     1592
     1593          flux_d    = flux_t(k-1)
     1594          diss_d    = diss_t(k-1)
     1595
    15831596!
    15841597!--       Calculate the divergence of the velocity field. A respective
     
    16061619          swap_flux_x_local(k,j,tn) = flux_r(k)
    16071620          swap_diss_x_local(k,j,tn) = diss_r(k)
    1608           flux_d                    = flux_t(k)
    1609           diss_d                    = diss_t(k)
    16101621
    16111622       ENDDO
     
    17651776       REAL(wp)    ::  gv       !<
    17661777       REAL(wp)    ::  u_comp_l !<
    1767        REAL(wp)    ::  v_comp   !<
    1768        REAL(wp)    ::  w_comp   !<
    17691778       
    17701779       REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_n !<
     
    17751784       REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_t !<
    17761785       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_comp !<
     1786       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_comp !<
     1787       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_comp !<
    17771788
    17781789       gu = 2.0_wp * u_gtrans
     
    17881799             ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
    17891800
    1790              v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    1791              flux_s_u(k,tn) = v_comp * (                                      &
     1801             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
     1802             flux_s_u(k,tn) = v_comp(k) * (                                   &
    17921803                            ( 37.0_wp * ibit14 * adv_mom_5                    &
    17931804                         +     7.0_wp * ibit13 * adv_mom_3                    &
     
    18021813                            ) *                                               &
    18031814                                        ( u(k,j+2,i) + u(k,j-3,i) )           &
    1804                                        )
    1805 
    1806              diss_s_u(k,tn) = - ABS ( v_comp ) * (                            &
     1815                                          )
     1816
     1817             diss_s_u(k,tn) = - ABS ( v_comp(k) ) * (                         &
    18071818                            ( 10.0_wp * ibit14 * adv_mom_5                    &
    18081819                         +     3.0_wp * ibit13 * adv_mom_3                    &
     
    18171828                            ) *                                               &
    18181829                                        ( u(k,j+2,i) - u(k,j-3,i) )           &
    1819                                                  )
     1830                                                    )
    18201831
    18211832          ENDDO
     
    18231834          DO  k = nzb_max+1, nzt
    18241835
    1825              v_comp         = v(k,j,i) + v(k,j,i-1) - gv
    1826              flux_s_u(k,tn) = v_comp * (                                      &
     1836             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
     1837             flux_s_u(k,tn) = v_comp(k) * (                                   &
    18271838                           37.0_wp * ( u(k,j,i) + u(k,j-1,i)   )              &
    18281839                         -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )              &
    18291840                         +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    1830              diss_s_u(k,tn) = - ABS(v_comp) * (                               &
     1841             diss_s_u(k,tn) = - ABS(v_comp(k)) * (                            &
    18311842                           10.0_wp * ( u(k,j,i) - u(k,j-1,i)   )              &
    18321843                         -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )              &
     
    18971908       flux_t(0) = 0.0_wp
    18981909       diss_t(0) = 0.0_wp
    1899        flux_d    = 0.0_wp
    1900        diss_d    = 0.0_wp
     1910       w_comp(0) = 0.0_wp
    19011911!
    19021912!--    Now compute the fluxes tendency terms for the horizontal and
     
    19431953          ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
    19441954
    1945           v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    1946           flux_n(k) = v_comp * (                                              &
     1955          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
     1956          flux_n(k) = v_comp(k) * (                                           &
    19471957                     ( 37.0_wp * ibit14 * adv_mom_5                           &
    19481958                  +     7.0_wp * ibit13 * adv_mom_3                           &
     
    19571967                     ) *                                                      &
    19581968                                    ( u(k,j+3,i) + u(k,j-2,i) )               &
    1959                                )
    1960 
    1961           diss_n(k) = - ABS ( v_comp ) * (                                    &
     1969                                  )
     1970
     1971          diss_n(k) = - ABS ( v_comp(k) ) * (                                 &
    19621972                     ( 10.0_wp * ibit14 * adv_mom_5                           &
    19631973                  +     3.0_wp * ibit13 * adv_mom_3                           &
     
    19721982                     ) *                                                      &
    19731983                                    ( u(k,j+3,i) - u(k,j-2,i) )               &
    1974                                          )
     1984                                            )
    19751985!
    19761986!--       k index has to be modified near bottom and top, else array
     
    19841994          k_mm  = k - 2 * ibit17
    19851995
    1986           w_comp    = w(k,j,i) + w(k,j,i-1)
    1987           flux_t(k) = w_comp * rho_air_zw(k) * (                              &
     1996          w_comp(k) = w(k,j,i) + w(k,j,i-1)
     1997          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    19881998                     ( 37.0_wp * ibit17 * adv_mom_5                           &
    19891999                  +     7.0_wp * ibit16 * adv_mom_3                           &
     
    19982008                     ) *                                                      &
    19992009                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
    2000                                  )
    2001 
    2002           diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
     2010                                                  )
     2011
     2012          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    20032013                     ( 10.0_wp * ibit17 * adv_mom_5                           &
    20042014                  +     3.0_wp * ibit16 * adv_mom_3                           &
     
    20132023                     ) *                                                      &
    20142024                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
    2015                                          )
     2025                                                           )
     2026       ENDDO
     2027
     2028       DO  k = nzb+1, nzb_max
     2029
     2030          flux_d    = flux_t(k-1)
     2031          diss_d    = diss_t(k-1)
    20162032!
    20172033!--       Calculate the divergence of the velocity field. A respective
     
    20252041                                      )                                       &
    20262042                  ) * ddx                                                     &
    2027                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
     2043               +  ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 )         &
    20282044                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
    20292045                                    * ( IBITS(advc_flags_1(k,j-1,i),12,1)     &
     
    20322048                                      )                                       &
    20332049                  ) * ddy                                                     &
    2034                +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
    2035                 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
     2050               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
     2051                -   w_comp(k-1) * rho_air_zw(k-1)                             &
    20362052                                    * ( IBITS(advc_flags_1(k-1,j,i),15,1)     &
    20372053                                      + IBITS(advc_flags_1(k-1,j,i),16,1)     &
     
    20562072           flux_s_u(k,tn)   = flux_n(k)
    20572073           diss_s_u(k,tn)   = diss_n(k)
    2058            flux_d           = flux_t(k)
    2059            diss_d           = diss_t(k)
    20602074!
    20612075!--        Statistical Evaluation of u'u'. The factor has to be applied for
     
    20732087           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
    20742088                + ( flux_t(k)                                                  &
    2075                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
    2076                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
     2089                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                )     &
     2090                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        )     &
    20772091                  + diss_t(k)                                                  &
    2078                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
    2079                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
     2092                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           )     &
     2093                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                    )     &
    20802094                  ) *   weight_substep(intermediate_timestep_count)
    20812095       ENDDO
     
    20932107                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    20942108
    2095           v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2096           flux_n(k) = v_comp * (                                              &
     2109          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
     2110          flux_n(k) = v_comp(k) * (                                           &
    20972111                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
    20982112                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
    20992113                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    2100           diss_n(k) = - ABS( v_comp ) * (                                     &
     2114          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    21012115                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
    21022116                       -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
     
    21132127          k_mm  = k - 2 * ibit17
    21142128
    2115           w_comp    = w(k,j,i) + w(k,j,i-1)
    2116           flux_t(k) = w_comp * rho_air_zw(k) * (                              &
     2129          w_comp(k) = w(k,j,i) + w(k,j,i-1)
     2130          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    21172131                     ( 37.0_wp * ibit17 * adv_mom_5                           &
    21182132                  +     7.0_wp * ibit16 * adv_mom_3                           &
     
    21272141                     ) *                                                      &
    21282142                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
    2129                                  )
    2130 
    2131           diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
     2143                                                  )
     2144
     2145          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    21322146                     ( 10.0_wp * ibit17 * adv_mom_5                           &
    21332147                  +     3.0_wp * ibit16 * adv_mom_3                           &
     
    21422156                     ) *                                                      &
    21432157                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
    2144                                          )
     2158                                                           )
     2159
     2160       ENDDO
     2161
     2162       DO  k = nzb_max+1, nzt
     2163
     2164          flux_d    = flux_t(k-1)
     2165          diss_d    = diss_t(k-1)
    21452166!
    21462167!--       Calculate the divergence of the velocity field. A respective
    21472168!--       correction is needed to overcome numerical instabilities introduced
    21482169!--       by a not sufficient reduction of divergences near topography.
    2149           div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
    2150                +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
    2151                +  (   w_comp                      * rho_air_zw(k)   -         &
    2152                     ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)           &
     2170          div = ( ( u_comp(k)      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx    &
     2171               +  ( v_comp(k) + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy    &
     2172               +  ( w_comp(k)   * rho_air_zw(k)   -                           &
     2173                    w_comp(k-1) * rho_air_zw(k-1)                             &
    21532174                  ) * drho_air(k) * ddzw(k)                                   &
    21542175                ) * 0.5_wp
     
    21682189           flux_s_u(k,tn)   = flux_n(k)
    21692190           diss_s_u(k,tn)   = diss_n(k)
    2170            flux_d           = flux_t(k)
    2171            diss_d           = diss_t(k)
    21722191!
    21732192!--        Statistical Evaluation of u'u'. The factor has to be applied for
     
    21852204           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
    21862205                + ( flux_t(k)                                                  &
    2187                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
    2188                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
     2206                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     2207                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
    21892208                  + diss_t(k)                                                  &
    2190                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
    2191                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
     2209                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
     2210                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
    21922211                  ) *   weight_substep(intermediate_timestep_count)
    21932212       ENDDO
     
    22542273       REAL(wp)     ::  gu       !<
    22552274       REAL(wp)     ::  gv       !<
    2256        REAL(wp)     ::  u_comp   !<
    22572275       REAL(wp)     ::  v_comp_l !<
    2258        REAL(wp)     ::  w_comp   !<
    22592276       
    22602277       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !<
     
    22642281       REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_r !<
    22652282       REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_t !<
     2283       REAL(wp), DIMENSION(nzb:nzt+1)  ::  u_comp !<
    22662284       REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !<
     2285       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !<
    22672286
    22682287       gu = 2.0_wp * u_gtrans
     
    22792298             ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
    22802299
    2281              u_comp           = u(k,j-1,i) + u(k,j,i) - gu
    2282              flux_l_v(k,j,tn) = u_comp * (                                    &
     2300             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
     2301             flux_l_v(k,j,tn) = u_comp(k) * (                                 &
    22832302                              ( 37.0_wp * ibit20 * adv_mom_5                  &
    22842303                           +     7.0_wp * ibit19 * adv_mom_3                  &
     
    22932312                              ) *                                             &
    22942313                                        ( v(k,j,i+2) + v(k,j,i-3) )           &
    2295                                          )
    2296 
    2297               diss_l_v(k,j,tn) = - ABS( u_comp ) * (                          &
     2314                                            )
     2315
     2316              diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                       &
    22982317                              ( 10.0_wp * ibit20 * adv_mom_5                  &
    22992318                           +     3.0_wp * ibit19 * adv_mom_3                  &
     
    23082327                              ) *                                             &
    23092328                                        ( v(k,j,i+2) - v(k,j,i-3) )           &
    2310                                                    )
     2329                                                      )
    23112330
    23122331          ENDDO
     
    23142333          DO  k = nzb_max+1, nzt
    23152334
    2316              u_comp           = u(k,j-1,i) + u(k,j,i) - gu
    2317              flux_l_v(k,j,tn) = u_comp * (                                    &
     2335             u_comp(k)           = u(k,j-1,i) + u(k,j,i) - gu
     2336             flux_l_v(k,j,tn) = u_comp(k) * (                                 &
    23182337                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )            &
    23192338                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )            &
    23202339                           +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    2321              diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
     2340             diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    23222341                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )            &
    23232342                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )            &
     
    23882407       flux_t(0) = 0.0_wp
    23892408       diss_t(0) = 0.0_wp
    2390        flux_d    = 0.0_wp
    2391        diss_d    = 0.0_wp
     2409       w_comp(0) = 0.0_wp
    23922410!
    23932411!--    Now compute the fluxes and tendency terms for the horizontal and
     
    23992417          ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
    24002418 
    2401           u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    2402           flux_r(k) = u_comp * (                                              &
     2419          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
     2420          flux_r(k) = u_comp(k) * (                                           &
    24032421                     ( 37.0_wp * ibit20 * adv_mom_5                           &
    24042422                  +     7.0_wp * ibit19 * adv_mom_3                           &
     
    24132431                     ) *                                                      &
    24142432                                    ( v(k,j,i+3) + v(k,j,i-2) )               &
    2415                                )
    2416 
    2417           diss_r(k) = - ABS( u_comp ) * (                                     &
     2433                                  )
     2434
     2435          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    24182436                     ( 10.0_wp * ibit20 * adv_mom_5                           &
    24192437                  +     3.0_wp * ibit19 * adv_mom_3                           &
     
    24282446                     ) *                                                      &
    24292447                                    ( v(k,j,i+3) - v(k,j,i-2) )               &
    2430                                         )
     2448                                           )
    24312449
    24322450          ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
     
    24762494          k_mm  = k - 2 * ibit26
    24772495
    2478           w_comp    = w(k,j-1,i) + w(k,j,i)
    2479           flux_t(k) = w_comp * rho_air_zw(k) * (                              &
     2496          w_comp(k) = w(k,j-1,i) + w(k,j,i)
     2497          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    24802498                     ( 37.0_wp * ibit26 * adv_mom_5                           &
    24812499                  +     7.0_wp * ibit25 * adv_mom_3                           &
     
    24902508                     ) *                                                      &
    24912509                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
    2492                                  )
    2493 
    2494           diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
     2510                                                  )
     2511
     2512          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    24952513                     ( 10.0_wp * ibit26 * adv_mom_5                           &
    24962514                  +     3.0_wp * ibit25 * adv_mom_3                           &
     
    25052523                     ) *                                                      &
    25062524                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
    2507                                          )
     2525                                                           )
     2526       ENDDO
     2527
     2528       DO  k = nzb+1, nzb_max
     2529
     2530          flux_d    = flux_t(k-1)
     2531          diss_d    = diss_t(k-1)
    25082532!
    25092533!--       Calculate the divergence of the velocity field. A respective
    25102534!--       correction is needed to overcome numerical instabilities introduced
    25112535!--       by a not sufficient reduction of divergences near topography.
    2512           div = ( ( ( u_comp     + gu )                                       &
     2536          div = ( ( ( u_comp(k)     + gu )                                    &
    25132537                                       * ( ibit18 + ibit19 + ibit20 )         &
    25142538                  - ( u(k,j-1,i) + u(k,j,i) )                                 &
     
    25262550                                         )                                    &
    25272551                  ) * ddy                                                     &
    2528                +  ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )     &
    2529                 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
     2552               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&
     2553                -   w_comp(k-1) * rho_air_zw(k-1)                             &
    25302554                                       * ( IBITS(advc_flags_1(k-1,j,i),24,1)  &
    25312555                                         + IBITS(advc_flags_1(k-1,j,i),25,1)  &
     
    25502574           flux_s_v(k,tn)   = flux_n(k)
    25512575           diss_s_v(k,tn)   = diss_n(k)
    2552            flux_d           = flux_t(k)
    2553            diss_d           = diss_t(k)
    25542576
    25552577!
     
    25682590           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
    25692591                + ( flux_t(k)                                                  &
    2570                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
    2571                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
     2592                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     2593                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
    25722594                  + diss_t(k)                                                  &
    2573                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
    2574                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
     2595                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
     2596                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
    25752597                  ) *   weight_substep(intermediate_timestep_count)
    25762598
     
    25792601       DO  k = nzb_max+1, nzt
    25802602
    2581           u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    2582           flux_r(k) = u_comp * (                                              &
     2603          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
     2604          flux_r(k) = u_comp(k) * (                                           &
    25832605                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
    25842606                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
    25852607                    +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    25862608
    2587           diss_r(k) = - ABS( u_comp ) * (                                     &
     2609          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    25882610                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
    25892611                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
     
    26122634          k_mm  = k - 2 * ibit26
    26132635
    2614           w_comp    = w(k,j-1,i) + w(k,j,i)
    2615           flux_t(k) = w_comp * rho_air_zw(k) * (                              &
     2636          w_comp(k) = w(k,j-1,i) + w(k,j,i)
     2637          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    26162638                     ( 37.0_wp * ibit26 * adv_mom_5                           &
    26172639                  +     7.0_wp * ibit25 * adv_mom_3                           &
     
    26262648                     ) *                                                      &
    26272649                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
    2628                                  )
    2629 
    2630           diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (                     &
     2650                                                  )
     2651
     2652          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    26312653                     ( 10.0_wp * ibit26 * adv_mom_5                           &
    26322654                  +     3.0_wp * ibit25 * adv_mom_3                           &
     
    26412663                     ) *                                                      &
    26422664                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
    2643                                          )
     2665                                                           )
     2666       ENDDO
     2667
     2668       DO  k = nzb_max+1, nzt
     2669
     2670          flux_d    = flux_t(k-1)
     2671          diss_d    = diss_t(k-1)
    26442672!
    26452673!--       Calculate the divergence of the velocity field. A respective
    26462674!--       correction is needed to overcome numerical instabilities introduced
    26472675!--       by a not sufficient reduction of divergences near topography.
    2648           div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
    2649                +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
    2650                +  (   w_comp                      * rho_air_zw(k)   -         &
    2651                     ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)           &
     2676          div = ( ( u_comp(k) + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx    &
     2677               +  ( v_comp(k)      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy    &
     2678               +  ( w_comp(k)   * rho_air_zw(k)   -                           &
     2679                    w_comp(k-1) * rho_air_zw(k-1)                             &
    26522680                  ) * drho_air(k) * ddzw(k)                                   &
    26532681                ) * 0.5_wp
     
    26672695           flux_s_v(k,tn)   = flux_n(k)
    26682696           diss_s_v(k,tn)   = diss_n(k)
    2669            flux_d           = flux_t(k)
    2670            diss_d           = diss_t(k)
    2671 
    26722697!
    26732698!--        Statistical Evaluation of v'v'. The factor has to be applied for
     
    26852710           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
    26862711                + ( flux_t(k)                                                  &
    2687                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
    2688                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
     2712                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
     2713                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
    26892714                  + diss_t(k)                                                  &
    2690                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
    2691                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
     2715                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
     2716                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
    26922717                  ) *   weight_substep(intermediate_timestep_count)
    26932718
     
    27542779       REAL(wp)    ::  gu      !<
    27552780       REAL(wp)    ::  gv      !<
    2756        REAL(wp)    ::  u_comp  !<
    2757        REAL(wp)    ::  v_comp  !<
    2758        REAL(wp)    ::  w_comp  !<
    27592781       
    27602782       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !<
     
    27642786       REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_r !<
    27652787       REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_t !<
     2788       REAL(wp), DIMENSION(nzb:nzt+1)  ::  u_comp !<
     2789       REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !<
     2790       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !<
    27662791
    27672792       gu = 2.0_wp * u_gtrans
     
    27772802             ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
    27782803
    2779              v_comp         = v(k+1,j,i) + v(k,j,i) - gv
    2780              flux_s_w(k,tn) = v_comp * (                                      &
     2804             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
     2805             flux_s_w(k,tn) = v_comp(k) * (                                   &
    27812806                            ( 37.0_wp * ibit32 * adv_mom_5                    &
    27822807                         +     7.0_wp * ibit31 * adv_mom_3                    &
     
    27912816                            ) *                                               &
    27922817                                        ( w(k,j+2,i) + w(k,j-3,i) )           &
    2793                                        )
    2794 
    2795              diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
     2818                                          )
     2819
     2820             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
    27962821                            ( 10.0_wp * ibit32 * adv_mom_5                    &
    27972822                         +     3.0_wp * ibit31 * adv_mom_3                    &
     
    28062831                            ) *                                               &
    28072832                                        ( w(k,j+2,i) - w(k,j-3,i) )           &
    2808                                                 )
     2833                                                   )
    28092834
    28102835          ENDDO
     
    28122837          DO  k = nzb_max+1, nzt
    28132838
    2814              v_comp         = v(k+1,j,i) + v(k,j,i) - gv
    2815              flux_s_w(k,tn) = v_comp * (                                      &
     2839             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
     2840             flux_s_w(k,tn) = v_comp(k) * (                                   &
    28162841                           37.0_wp * ( w(k,j,i) + w(k,j-1,i)   )              &
    28172842                         -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )              &
    28182843                         +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    2819              diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
     2844             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                         &
    28202845                           10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )              &
    28212846                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )              &
     
    28352860             ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
    28362861
    2837              u_comp           = u(k+1,j,i) + u(k,j,i) - gu
    2838              flux_l_w(k,j,tn) = u_comp * (                                    &
     2862             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
     2863             flux_l_w(k,j,tn) = u_comp(k) * (                                 &
    28392864                             ( 37.0_wp * ibit29 * adv_mom_5                   &
    28402865                          +     7.0_wp * ibit28 * adv_mom_3                   &
     
    28492874                             ) *                                              &
    28502875                                        ( w(k,j,i+2) + w(k,j,i-3) )           &
    2851                                          )
    2852 
    2853                diss_l_w(k,j,tn) = - ABS( u_comp ) * (                         &
     2876                                            )
     2877
     2878               diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                      &
    28542879                             ( 10.0_wp * ibit29 * adv_mom_5                   &
    28552880                          +     3.0_wp * ibit28 * adv_mom_3                   &
     
    28642889                             ) *                                              &
    28652890                                        ( w(k,j,i+2) - w(k,j,i-3) )           &
    2866                                                     )
     2891                                                       )
    28672892
    28682893          ENDDO
     
    28702895          DO  k = nzb_max+1, nzt
    28712896
    2872              u_comp           = u(k+1,j,i) + u(k,j,i) - gu
    2873              flux_l_w(k,j,tn) = u_comp * (                                    &
     2897             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
     2898             flux_l_w(k,j,tn) = u_comp(k) * (                                 &
    28742899                            37.0_wp * ( w(k,j,i) + w(k,j,i-1)   )             &
    28752900                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )             &
    28762901                          +           ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    2877              diss_l_w(k,j,tn) = - ABS( u_comp ) * (                           &
     2902             diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    28782903                            10.0_wp * ( w(k,j,i) - w(k,j,i-1)   )             &
    28792904                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )             &
     
    28882913!--    advc_flags_1.
    28892914       k         = nzb + 1
    2890        w_comp    = w(k,j,i) + w(k-1,j,i)
    2891        flux_t(0) = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
    2892        diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
    2893        flux_d    = flux_t(0)
    2894        diss_d    = diss_t(0)
     2915       w_comp(k) = w(k,j,i) + w(k-1,j,i)
     2916       flux_t(0) = w_comp(k)       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
     2917       diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
    28952918!
    28962919!--    Now compute the fluxes and tendency terms for the horizontal
     
    29022925          ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
    29032926
    2904           u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    2905           flux_r(k) = u_comp * (                                              &
     2927          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
     2928          flux_r(k) = u_comp(k) * (                                           &
    29062929                     ( 37.0_wp * ibit29 * adv_mom_5                           &
    29072930                  +     7.0_wp * ibit28 * adv_mom_3                           &
     
    29162939                     ) *                                                      &
    29172940                                    ( w(k,j,i+3) + w(k,j,i-2) )               &
    2918                                )
    2919 
    2920           diss_r(k) = - ABS( u_comp ) * (                                     &
     2941                                  )
     2942
     2943          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    29212944                     ( 10.0_wp * ibit29 * adv_mom_5                           &
    29222945                  +     3.0_wp * ibit28 * adv_mom_3                           &
     
    29312954                     ) *                                                      &
    29322955                                    ( w(k,j,i+3) - w(k,j,i-2) )               &
    2933                                         )
     2956                                           )
    29342957
    29352958          ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
     
    29372960          ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
    29382961
    2939           v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    2940           flux_n(k) = v_comp * (                                              &
     2962          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
     2963          flux_n(k) = v_comp(k) * (                                           &
    29412964                     ( 37.0_wp * ibit32 * adv_mom_5                           &
    29422965                  +     7.0_wp * ibit31 * adv_mom_3                           &
     
    29512974                     ) *                                                      &
    29522975                                    ( w(k,j+3,i) + w(k,j-2,i) )               &
    2953                                )
    2954 
    2955           diss_n(k) = - ABS( v_comp ) * (                                     &
     2976                                  )
     2977
     2978          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    29562979                     ( 10.0_wp * ibit32 * adv_mom_5                           &
    29572980                  +     3.0_wp * ibit31 * adv_mom_3                           &
     
    29662989                     ) *                                                      &
    29672990                                   ( w(k,j+3,i) - w(k,j-2,i) )                &
    2968                                         )
     2991                                           )
    29692992!
    29702993!--       k index has to be modified near bottom and top, else array
     
    29783001          k_mm  = k - 2 * ibit35
    29793002
    2980           w_comp    = w(k+1,j,i) + w(k,j,i)
    2981           flux_t(k) = w_comp * rho_air(k+1) * (                               &
     3003          w_comp(k) = w(k+1,j,i) + w(k,j,i)
     3004          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    29823005                     ( 37.0_wp * ibit35 * adv_mom_5                           &
    29833006                  +     7.0_wp * ibit34 * adv_mom_3                           &
     
    29923015                     ) *                                                      &
    29933016                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
    2994                                 )
    2995 
    2996           diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
     3017                                                 )
     3018
     3019          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    29973020                     ( 10.0_wp * ibit35 * adv_mom_5                           &
    29983021                  +     3.0_wp * ibit34 * adv_mom_3                           &
     
    30073030                     ) *                                                      &
    30083031                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
    3009                                         )
    3010 
     3032                                                          )
     3033       ENDDO
     3034
     3035       DO  k = nzb+1, nzb_max
     3036
     3037          flux_d    = flux_t(k-1)
     3038          diss_d    = diss_t(k-1)
    30113039!
    30123040!--       Calculate the divergence of the velocity field. A respective
    30133041!--       correction is needed to overcome numerical instabilities introduced
    30143042!--       by a not sufficient reduction of divergences near topography.
    3015           div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )            &
     3043          div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 )         &
    30163044                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
    30173045                                    * ( IBITS(advc_flags_1(k,j,i-1),27,1)     &
     
    30203048                                      )                                       &
    30213049                  ) * ddx                                                     &
    3022               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
     3050              +   ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 )         &
    30233051                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
    30243052                                    * ( IBITS(advc_flags_1(k,j-1,i),30,1)     &
     
    30273055                                      )                                       &
    30283056                  ) * ddy                                                     &
    3029               +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
    3030                 - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
     3057              +   ( w_comp(k)               * rho_air(k+1)                    &
     3058                                            * ( ibit33 + ibit34 + ibit35 )    &
     3059                - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
    30313060                                    * ( IBITS(advc_flags_2(k-1,j,i),1,1)      &
    30323061                                      + IBITS(advc_flags_2(k-1,j,i),2,1)      &
     
    30513080          flux_s_w(k,tn)   = flux_n(k)
    30523081          diss_s_w(k,tn)   = diss_n(k)
    3053           flux_d           = flux_t(k)
    3054           diss_d           = diss_t(k)
    30553082!
    30563083!--       Statistical Evaluation of w'w'.
    30573084          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
    30583085                      + ( flux_t(k)                                           &
    3059                        * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
    3060                        / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
     3086                       * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                ) &
     3087                       / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        ) &
    30613088                        + diss_t(k)                                           &
    3062                        *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
    3063                        / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
     3089                       *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           ) &
     3090                       / ( ABS( w_comp(k) ) + 1.0E-20_wp                    ) &
    30643091                        ) *   weight_substep(intermediate_timestep_count)
    30653092
     
    30683095       DO  k = nzb_max+1, nzt
    30693096
    3070           u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    3071           flux_r(k) = u_comp * (                                              &
     3097          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
     3098          flux_r(k) = u_comp(k) * (                                           &
    30723099                      37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
    30733100                    -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
    30743101                    +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    30753102
    3076           diss_r(k) = - ABS( u_comp ) * (                                     &
     3103          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    30773104                      10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
    30783105                    -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
    30793106                    +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    30803107
    3081           v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    3082           flux_n(k) = v_comp * (                                              &
     3108          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
     3109          flux_n(k) = v_comp(k) * (                                           &
    30833110                      37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
    30843111                    -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
    30853112                    +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    30863113
    3087           diss_n(k) = - ABS( v_comp ) * (                                     &
     3114          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    30883115                      10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
    30893116                    -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
     
    31003127          k_mm  = k - 2 * ibit35
    31013128
    3102           w_comp    = w(k+1,j,i) + w(k,j,i)
    3103           flux_t(k) = w_comp * rho_air(k+1) * (                               &
     3129          w_comp(k) = w(k+1,j,i) + w(k,j,i)
     3130          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    31043131                     ( 37.0_wp * ibit35 * adv_mom_5                           &
    31053132                  +     7.0_wp * ibit34 * adv_mom_3                           &
     
    31143141                     ) *                                                      &
    31153142                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
    3116                                 )
    3117 
    3118           diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (                      &
     3143                                                 )
     3144
     3145          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    31193146                     ( 10.0_wp * ibit35 * adv_mom_5                           &
    31203147                  +     3.0_wp * ibit34 * adv_mom_3                           &
     
    31293156                     ) *                                                      &
    31303157                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
    3131                                         )
     3158                                                          )
     3159       ENDDO
     3160
     3161       DO  k = nzb_max+1, nzt
     3162
     3163          flux_d    = flux_t(k-1)
     3164          diss_d    = diss_t(k-1)
    31323165!
    31333166!--       Calculate the divergence of the velocity field. A respective
    31343167!--       correction is needed to overcome numerical instabilities introduced
    31353168!--       by a not sufficient reduction of divergences near topography.
    3136           div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
    3137               +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
    3138               +   (   w_comp                    * rho_air(k+1) -              &
    3139                     ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)                  &
     3169          div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx      &
     3170              +   ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy      &
     3171              +   ( w_comp(k)               * rho_air(k+1) -                  &
     3172                  ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
    31403173                  ) * drho_air_zw(k) * ddzu(k+1)                              &
    31413174                ) * 0.5_wp
     
    31553188          flux_s_w(k,tn)   = flux_n(k)
    31563189          diss_s_w(k,tn)   = diss_n(k)
    3157           flux_d           = flux_t(k)
    3158           diss_d           = diss_t(k)
    31593190!
    31603191!--       Statistical Evaluation of w'w'.
    31613192          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
    31623193                      + ( flux_t(k)                                           &
    3163                        * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
    3164                        / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
     3194                       * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                ) &
     3195                       / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        ) &
    31653196                        + diss_t(k)                                           &
    3166                        *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
    3167                        / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
     3197                       *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           ) &
     3198                       / ( ABS( w_comp(k) ) + 1.0E-20_wp                    ) &
    31683199                        ) *   weight_substep(intermediate_timestep_count)
    31693200
Note: See TracChangeset for help on using the changeset viewer.