Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 1961)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 1996)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix concerning calculation of turbulent of turbulent fluxes
 ! 
 ! Former revisions:
@@ -916,5 +916,5 @@
 
        USE statistics,                                                         &
-           ONLY:  sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l,           &
+           ONLY:  hom, sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l,      &
                   sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsss_ws_l,             &
                   weight_substep
@@ -1649,7 +1649,10 @@
 
              DO  k = nzb, nzt
-                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +               &
-                                       ( flux_t(k) + diss_t(k) )              &
-                               * weight_substep(intermediate_timestep_count)
+                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
              ENDDO
             
@@ -1657,7 +1660,10 @@
 
              DO  k = nzb, nzt
-                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +               &
-                                       ( flux_t(k) + diss_t(k) )              &
-                               * weight_substep(intermediate_timestep_count)
+                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
              ENDDO
             
@@ -1665,7 +1671,10 @@
 
              DO  k = nzb, nzt
-                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                &
-                                      ( flux_t(k) + diss_t(k) )               &
-                               * weight_substep(intermediate_timestep_count)
+                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                 &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
              ENDDO
 
@@ -1673,7 +1682,10 @@
 
              DO  k = nzb, nzt
-                sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +              &
-                                      ( flux_t(k) + diss_t(k) )               &
-                               * weight_substep(intermediate_timestep_count)
+                sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +               &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
              ENDDO
 
@@ -1681,7 +1693,10 @@
 
              DO  k = nzb, nzt
-                sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +              &
-                                      ( flux_t(k) + diss_t(k) )               &
-                               * weight_substep(intermediate_timestep_count)
+                sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +               &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
              ENDDO
              
@@ -1689,7 +1704,10 @@
           
              DO  k = nzb, nzt
-                sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn) +                &
-                                      ( flux_t(k) + diss_t(k) )               &
-                               * weight_substep(intermediate_timestep_count)
+                sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn) +                 &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
              ENDDO
 
@@ -2049,17 +2067,22 @@
 !--        Statistical Evaluation of u'u'. The factor has to be applied for
 !--        right evaluation when gallilei_trans = .T. .
-           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
-                              + ( flux_r(k) *                                 &
-                                ( u_comp(k) - 2.0_wp * hom(k,1,1,0) )         &
-                              / ( u_comp(k) - gu + 1.0E-20_wp   )             &
-                              +   diss_r(k) *                                 &
-                                  ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) )    &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
-                              *   weight_substep(intermediate_timestep_count)
+           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                           &
+                + ( flux_r(k)                                                  &
+                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
+                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
+                  + diss_r(k)                                                  &
+                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
+                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
 !
 !--        Statistical Evaluation of w'u'.
-           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                        &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
+           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                  + diss_t(k)                                                  &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
        ENDDO
 
@@ -2071,15 +2094,15 @@
                        -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
                        +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-             diss_r(k) = - ABS( u_comp(k) - gu ) * (                          &
+          diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
                          10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
                        -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
                        +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
 
-             v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-             flux_n(k) = v_comp * (                                           &
+          v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
+          flux_n(k) = v_comp * (                                              &
                          37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
                        -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
                        +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
-             diss_n(k) = - ABS( v_comp ) * (                                  &
+          diss_n(k) = - ABS( v_comp ) * (                                     &
                          10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
                        -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
@@ -2153,17 +2176,22 @@
 !--        Statistical Evaluation of u'u'. The factor has to be applied for
 !--        right evaluation when gallilei_trans = .T. .
-           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
-                              + ( flux_r(k) *                                 &
-                                ( u_comp(k) - 2.0_wp * hom(k,1,1,0)      )    &
-                              / ( u_comp(k) - gu + 1.0E-20_wp          )      &
-                              +   diss_r(k) *                                 &
-                                  ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) )    &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
-                              *   weight_substep(intermediate_timestep_count)
+           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                           &
+                + ( flux_r(k)                                                  &
+                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
+                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
+                  + diss_r(k)                                                  &
+                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
+                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
 !
 !--        Statistical Evaluation of w'u'.
-           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                        &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
+           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                  + diss_t(k)                                                  &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
        ENDDO
 
@@ -2530,17 +2558,22 @@
 !--        Statistical Evaluation of v'v'. The factor has to be applied for
 !--        right evaluation when gallilei_trans = .T. .
-           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
-             + ( flux_n(k)                                                    &
-             * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)      )                     &
-             / ( v_comp(k) - gv + 1.0E-20_wp       )                          &
-             +   diss_n(k)                                                    &
-             *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) )                     &
-             / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )                        &
-             *   weight_substep(intermediate_timestep_count)
-!
-!--        Statistical Evaluation of w'v'.
-           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                        &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
+           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                           &
+                + ( flux_n(k)                                                  &
+                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
+                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
+                  + diss_n(k)                                                  &
+                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
+                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
+!
+!--        Statistical Evaluation of w'u'.
+           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                  + diss_t(k)                                                  &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
 
        ENDDO
@@ -2639,17 +2672,22 @@
 !--        Statistical Evaluation of v'v'. The factor has to be applied for
 !--        right evaluation when gallilei_trans = .T. .
-           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
-             + ( flux_n(k)                                                    &
-             * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)      )                     &
-             / ( v_comp(k) - gv + 1.0E-20_wp       )                          &
-             +   diss_n(k)                                                    &
-             *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) )                     &
-             / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )                        &
-             *   weight_substep(intermediate_timestep_count)
-!
-!--        Statistical Evaluation of w'v'.
-           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                        &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
+           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                           &
+                + ( flux_n(k)                                                  &
+                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
+                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
+                  + diss_n(k)                                                  &
+                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
+                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
+!
+!--        Statistical Evaluation of w'u'.
+           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                  + diss_t(k)                                                  &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
 
        ENDDO
@@ -3013,8 +3051,13 @@
           diss_d           = diss_t(k)
 !
-!--        Statistical Evaluation of w'w'.
+!--       Statistical Evaluation of w'w'.
           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
-                             + ( flux_t(k) + diss_t(k) )                      &
-                             *   weight_substep(intermediate_timestep_count)
+                      + ( flux_t(k)                                           &
+                       * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
+                       / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
+                        + diss_t(k)                                           &
+                       *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
+                       / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
+                        ) *   weight_substep(intermediate_timestep_count)
 
        ENDDO
@@ -3109,8 +3152,13 @@
           diss_d           = diss_t(k)
 !
-!--        Statistical Evaluation of w'w'.
+!--       Statistical Evaluation of w'w'.
           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
-                             + ( flux_t(k) + diss_t(k) )                      &
-                             *   weight_substep(intermediate_timestep_count)
+                      + ( flux_t(k)                                           &
+                       * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
+                       / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
+                        + diss_t(k)                                           &
+                       *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
+                       / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
+                        ) *   weight_substep(intermediate_timestep_count)
 
        ENDDO
@@ -3147,5 +3195,5 @@
        
        USE statistics,                                                        &
-           ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,           &
+           ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,      &
                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, sums_wsss_ws_l,           &
                   weight_substep
@@ -3867,37 +3915,67 @@
                  CASE ( 'pt' )
                     DO  k = nzb, nzt
-                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)          &
-                        + ( flux_t(k) + diss_t(k) )                           &
-                        *   weight_substep(intermediate_timestep_count)
+                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)           &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
                     ENDDO
                  CASE ( 'sa' )
                     DO  k = nzb, nzt
-                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)          &
-                        + ( flux_t(k) + diss_t(k) )                           &
-                        *   weight_substep(intermediate_timestep_count)
+                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)           &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
                     ENDDO
                  CASE ( 'q' )
                     DO  k = nzb, nzt
-                       sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)            &
-                        + ( flux_t(k) + diss_t(k) )                           &
-                        *   weight_substep(intermediate_timestep_count)
+                       sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn)            &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
                     ENDDO
                  CASE ( 'qr' )
                     DO  k = nzb, nzt
-                       sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)          &
-                        + ( flux_t(k) + diss_t(k) )                           &
-                        *   weight_substep(intermediate_timestep_count)
+                       sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn)          &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
                     ENDDO
                  CASE ( 'nr' )
                     DO  k = nzb, nzt
-                       sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)          &
-                        + ( flux_t(k) + diss_t(k) )                           &
-                        *   weight_substep(intermediate_timestep_count)
+                       sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn)          &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
                     ENDDO
                  CASE ( 's' )
                     DO  k = nzb, nzt
-                       sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn)          &
-                        + ( flux_t(k) + diss_t(k) )                           &
-                        *   weight_substep(intermediate_timestep_count)
+                       sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn)            &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
                     ENDDO   
                                    
@@ -4772,17 +4850,23 @@
 !--             Statistical Evaluation of u'u'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
-                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                     &
-                              + ( flux_r(k) *                                 &
-                                ( u_comp(k) - 2.0_wp * hom(k,1,1,0) )            &
-                              / ( u_comp(k) - gu + 1.0E-20_wp    )            &
-                              +   diss_r(k) *                                 &
-                                  ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) )       &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
-                              *   weight_substep(intermediate_timestep_count)
+                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                      &
+                + ( flux_r(k)                                                  &
+                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
+                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
+                  + diss_r(k)                                                  &
+                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
+                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
 !
 !--             Statistical Evaluation of w'u'.
-                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                   &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
+                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                    &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                  + diss_t(k)                                                  &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
+
              ENDDO
 
@@ -4859,5 +4943,5 @@
                       ) * 0.5_wp
 
-                 tend(k,j,i) = tend(k,j,i) - (                                 &
+                tend(k,j,i) = tend(k,j,i) - (                                  &
                  ( flux_r(k) + diss_r(k)                                       &
                -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
@@ -4869,27 +4953,32 @@
                                            ) + div * u(k,j,i)
 
-                 swap_flux_x_local_u(k,j) = flux_r(k)
-                 swap_diss_x_local_u(k,j) = diss_r(k)
-                 swap_flux_y_local_u(k)   = flux_n(k)
-                 swap_diss_y_local_u(k)   = diss_n(k)
-                 flux_d                   = flux_t(k)
-                 diss_d                   = diss_t(k)
-!
-!--              Statistical Evaluation of u'u'. The factor has to be applied
-!--              for right evaluation when gallilei_trans = .T. .
-                 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                    &
-                              + ( flux_r(k) *                                 &
-                                ( u_comp(k) - 2.0_wp * hom(k,1,1,0) )            &
-                              / ( u_comp(k) - gu + 1.0E-20_wp   )             &
-                              +   diss_r(k) *                                 &
-                                  ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) )       &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
-                              *   weight_substep(intermediate_timestep_count)
-!
-!--              Statistical Evaluation of w'u'.
-                 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                  &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
-       ENDDO
+                swap_flux_x_local_u(k,j) = flux_r(k)
+                swap_diss_x_local_u(k,j) = diss_r(k)
+                swap_flux_y_local_u(k)   = flux_n(k)
+                swap_diss_y_local_u(k)   = diss_n(k)
+                flux_d                   = flux_t(k)
+                diss_d                   = diss_t(k)
+!
+!--             Statistical Evaluation of u'u'. The factor has to be applied
+!--             for right evaluation when gallilei_trans = .T. .
+                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                      &
+                + ( flux_r(k)                                                  &
+                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
+                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
+                  + diss_r(k)                                                  &
+                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
+                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
+!
+!--             Statistical Evaluation of w'u'.
+                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                    &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+                  + diss_t(k)                                                  &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
+             ENDDO
           ENDDO
        ENDDO
@@ -5619,17 +5708,22 @@
 !--             Statistical Evaluation of v'v'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
-                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                     &
-                      + ( flux_n(k)                                           &
-                      * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) )                    &
-                      / ( v_comp(k) - gv + 1.0E-20_wp )                       &
-                      +   diss_n(k)                                           &
-                      *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) )               &
-                      / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )               &
-                      *   weight_substep(intermediate_timestep_count)
-!
-!--              Statistical Evaluation of w'v'.
-                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                  &
-                              + ( flux_t(k) + diss_t(k) )                     &
-                              *   weight_substep(intermediate_timestep_count)
+                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                      &
+                + ( flux_n(k)                                                  &
+                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
+                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
+               +   diss_n(k)                                                   &
+                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
+                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
+!
+!--             Statistical Evaluation of w'u'.
+                sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                    &
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+               +   diss_t(k)                                                   &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
 
              ENDDO
@@ -5732,17 +5826,22 @@
 !--             Statistical Evaluation of v'v'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
-                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                     &
-                         + ( flux_n(k)                                        &
-                         * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) )                 &
-                         / ( v_comp(k) - gv + 1.0E-20_wp )                    &
-                         +   diss_n(k)                                        &
-                         *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) )            &
-                         / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )            &
-                         *   weight_substep(intermediate_timestep_count)
-!
-!--             Statistical Evaluation of w'v'.
+                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                      &
+                + ( flux_n(k)                                                  &
+                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
+                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
+               +   diss_n(k)                                                   &
+                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
+                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
+                  ) *   weight_substep(intermediate_timestep_count)
+!
+!--             Statistical Evaluation of w'u'.
                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                    &
-                              + ( flux_t(k) + diss_t(k) )                      &
-                              *   weight_substep(intermediate_timestep_count)
+                + ( flux_t(k)                                                  &
+                    * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
+                    / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
+               +   diss_t(k)                                                   &
+                    *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
+                    / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
+                  ) *   weight_substep(intermediate_timestep_count)
 
              ENDDO
@@ -6481,6 +6580,11 @@
 
                 sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
-                               + ( flux_t(k) + diss_t(k) )                    &
-                               *   weight_substep(intermediate_timestep_count)
+                      + ( flux_t(k)                                           &
+                       * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
+                       / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
+                        + diss_t(k)                                           &
+                       *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
+                       / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
+                        ) *   weight_substep(intermediate_timestep_count)
 
              ENDDO
@@ -6580,6 +6684,11 @@
 
                 sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
-                               + ( flux_t(k) + diss_t(k) )                    &
-                               *   weight_substep(intermediate_timestep_count)
+                      + ( flux_t(k)                                           &
+                       * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
+                       / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
+                        + diss_t(k)                                           &
+                       *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
+                       / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
+                        ) *   weight_substep(intermediate_timestep_count)
 
              ENDDO
