Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 800)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 801)
@@ -4,4 +4,8 @@
 ! Current revisions:
 ! ------------------
+! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
+! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
+! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
+! dimension.
 !
 ! Former revisions:
@@ -135,7 +139,9 @@
        IF ( ws_scheme_mom )  THEN
 
-          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1), sums_wsvs_ws_l(nzb:nzt+1),  &
-                    sums_us2_ws_l(nzb:nzt+1),  sums_vs2_ws_l(nzb:nzt+1),   &
-                    sums_ws2_ws_l(nzb:nzt+1) )
+          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1),  &
+                    sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1),  &
+                    sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1),   &
+                    sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1),   &
+                    sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
 
           sums_wsus_ws_l = 0.0
@@ -149,14 +155,14 @@
        IF ( ws_scheme_sca )  THEN
 
-          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1) )
+          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
           sums_wspts_ws_l = 0.0
 
           IF ( humidity .OR. passive_scalar )  THEN
-             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1) )
+             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              sums_wsqs_ws_l = 0.0
           ENDIF
 
           IF ( ocean )  THEN
-             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1) )
+             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              sums_wssas_ws_l = 0.0
           ENDIF
@@ -643,5 +649,5 @@
 
              DO  k = nzb_s_inner(j,i), nzt
-                sums_wspts_ws_l(k) = sums_wspts_ws_l(k) +                      &
+                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                      &
                                        ( flux_t(k) + diss_t(k) )               &
                                  * weight_substep(intermediate_timestep_count)
@@ -651,5 +657,5 @@
 
              DO  k = nzb_s_inner(j,i), nzt
-                sums_wssas_ws_l(k) = sums_wssas_ws_l(k) +                      &
+                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                      &
                                        ( flux_t(k) + diss_t(k) )               &
                                  * weight_substep(intermediate_timestep_count)
@@ -659,5 +665,5 @@
 
              DO  k = nzb_s_inner(j,i), nzt
-                sums_wsqs_ws_l(k)  = sums_wsqs_ws_l(k) +                       &
+                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                       &
                                       ( flux_t(k) + diss_t(k) )                &
                                  * weight_substep(intermediate_timestep_count)
@@ -905,5 +911,5 @@
 !--        Statistical Evaluation of u'u'. The factor has to be applied for 
 !--        right evaluation when gallilei_trans = .T. .
-           sums_us2_ws_l(k) = sums_us2_ws_l(k)                                &
+           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                                &
              + ( flux_r(k) *                                                  &
                ( u_comp(k) - 2.0 * hom(k,1,1,0) )                             &
@@ -914,5 +920,5 @@
              *   weight_substep(intermediate_timestep_count)
        ENDDO 
-       sums_us2_ws_l(nzb_u_inner(j,i)) = sums_us2_ws_l(nzb_u_inner(j,i)+1)
+       sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
                                            
 
@@ -1000,5 +1006,5 @@
 !--    sum up the vertical momentum fluxes
        DO  k = nzb_u_inner(j,i), nzt
-          sums_wsus_ws_l(k) = sums_wsus_ws_l(k)                               &
+          sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                               &
               + ( flux_t(k) + diss_t(k) )                                     &
               *   weight_substep(intermediate_timestep_count)
@@ -1243,5 +1249,5 @@
 !--        right evaluation when gallilei_trans = .T. .
 
-           sums_vs2_ws_l(k) = sums_vs2_ws_l(k)                                &
+           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                                &
              + ( flux_n(k)                                                    &
              * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
@@ -1253,5 +1259,5 @@
 
        ENDDO
-       sums_vs2_ws_l(nzb_v_inner(j,i)) = sums_vs2_ws_l(nzb_v_inner(j,i)+1)  
+       sums_vs2_ws_l(nzb_v_inner(j,i),tn) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,tn)  
                                            
 !
@@ -1335,5 +1341,5 @@
               
        DO  k = nzb_v_inner(j,i), nzt
-          sums_wsvs_ws_l(k) = sums_wsvs_ws_l(k)                               &
+          sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                               &
                  + ( flux_t(k) + diss_t(k) )                                  &
                  *   weight_substep(intermediate_timestep_count) 
@@ -1667,5 +1673,5 @@
        
        DO  k = nzb_w_inner(j,i), nzt
-           sums_ws2_ws_l(k)  = sums_ws2_ws_l(k)                               &
+           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                               &
                  + ( flux_t(k) + diss_t(k) )                                  &
                  *   weight_substep(intermediate_timestep_count)
@@ -1689,6 +1695,5 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
-
+       INTEGER ::  i, j, k, tn = 0
        REAL, DIMENSION(:,:,:), POINTER ::  sk
        REAL    :: flux_d, diss_d, u_comp, v_comp
@@ -2039,5 +2044,5 @@
                CASE ( 'pt' )
                  DO  k = nzb_s_inner(j,i), nzt
-                   sums_wspts_ws_l(k) = sums_wspts_ws_l(k)                    &
+                   sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)                    &
                       + ( flux_t(k) + diss_t(k) )                             &
                       *   weight_substep(intermediate_timestep_count)
@@ -2045,5 +2050,5 @@
                CASE ( 'sa' )
                  DO  k = nzb_s_inner(j,i), nzt
-                   sums_wssas_ws_l(k) = sums_wssas_ws_l(k)                    &
+                   sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)                    &
                       + ( flux_t(k) + diss_t(k) )                             &
                       *   weight_substep(intermediate_timestep_count)
@@ -2051,5 +2056,5 @@
                CASE ( 'q' )
                  DO  k = nzb_s_inner(j,i), nzt
-                   sums_wsqs_ws_l(k) = sums_wsqs_ws_l(k)                      &
+                   sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)                      &
                       + ( flux_t(k) + diss_t(k) )                             &
                       *   weight_substep(intermediate_timestep_count)
@@ -2078,5 +2083,5 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
+       INTEGER ::  i, j, k, tn = 0
        REAL    ::  gu, gv, flux_d, diss_d, v_comp, w_comp
        REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u
@@ -2332,5 +2337,5 @@
                swap_diss_y_local_u(k)   = diss_n(k)
                       
-               sums_us2_ws_l(k)  = sums_us2_ws_l(k)                            &
+               sums_us2_ws_l(k,tn)  = sums_us2_ws_l(k,tn)                            &
                  + ( flux_r(k)                                                 &
                  * ( u_comp(k) - 2.0 * hom(k,1,1,0) )                          &
@@ -2341,5 +2346,5 @@
                  *   weight_substep(intermediate_timestep_count)
             ENDDO
-            sums_us2_ws_l(nzb_u_inner(j,i)) = sums_us2_ws_l(nzb_u_inner(j,i)+1)
+            sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
          ENDDO
        ENDDO
@@ -2434,5 +2439,5 @@
 !-- at last vertical momentum flux is accumulated
             DO  k = nzb_u_inner(j,i), nzt
-               sums_wsus_ws_l(k) = sums_wsus_ws_l(k)                           &
+               sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                           &
                               + ( flux_t(k) + diss_t(k) )                      &
                               *   weight_substep(intermediate_timestep_count)
@@ -2460,5 +2465,5 @@
 
 
-       INTEGER ::  i, j, k
+       INTEGER ::  i, j, k, tn = 0
        REAL    ::  gu, gv, flux_l, flux_s, flux_d, diss_l, diss_s, diss_d,    &
                    u_comp, w_comp
@@ -2711,5 +2716,5 @@
                swap_diss_y_local_v(k)   = diss_n(k)   
 
-               sums_vs2_ws_l(k) = sums_vs2_ws_l(k)                             &
+               sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                             &
                   + ( flux_n(k) * ( v_comp(k) - 2.0 * hom(k,1,2,0) )           &
                   / ( v_comp(k) - gv + 1.0E-20 )                               &
@@ -2718,5 +2723,5 @@
                   *   weight_substep(intermediate_timestep_count)
             ENDDO
-            sums_vs2_ws_l(nzb_v_inner(j,i)) = sums_vs2_ws_l(nzb_v_inner(j,i)+1)
+            sums_vs2_ws_l(nzb_v_inner(j,i),tn) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,tn)
          ENDDO
        ENDDO
@@ -2808,9 +2813,9 @@
 !-          At last vertical momentum flux is accumulated.
             DO  k = nzb_v_inner(j,i), nzt
-               sums_wsvs_ws_l(k) = sums_wsvs_ws_l(k)                           &
+               sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                           &
                                + ( flux_t(k) + diss_t(k) )                     &
                                *   weight_substep(intermediate_timestep_count)
             ENDDO
-            sums_vs2_ws_l(nzb_v_inner(j,i)) = sums_vs2_ws_l(nzb_v_inner(j,i)+1)
+            sums_vs2_ws_l(nzb_v_inner(j,i),tn) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,tn)
           ENDDO
        ENDDO
@@ -2833,5 +2838,5 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
+       INTEGER ::  i, j, k, tn = 0
        REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
        REAL    ::  flux_t(nzb:nzt+1), diss_t(nzb:nzt+1)
@@ -3184,5 +3189,5 @@
 !--         at last vertical momentum flux is accumulated
             DO  k = nzb_w_inner(j,i), nzt
-               sums_ws2_ws_l(k)  = sums_ws2_ws_l(k)                            &
+               sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                            &
                                + ( flux_t(k) + diss_t(k) )                     &
                                *   weight_substep(intermediate_timestep_count)
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 800)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 801)
@@ -4,6 +4,5 @@
 ! Current revisions:
 ! -----------------
-! Calculation of turbulent fluxes with WS-scheme only for the whole model
-! domain, not for user-defined subregions. 
+! Calculation of turbulent fluxes in advec_ws is now thread-safe.
 !
 ! Former revisions:
@@ -179,6 +178,6 @@
 !--       the corresponding fluxes has to satisfiy the same bc.
           IF ( ocean )  THEN
-             sums_us2_ws_l(nzt+1) = sums_us2_ws_l(nzt)
-             sums_vs2_ws_l(nzt+1) = sums_vs2_ws_l(nzt)    
+             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
+             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)    
           ENDIF
 
@@ -186,17 +185,17 @@
 !          
 !--          Swap the turbulent quantities evaluated in advec_ws.
-             sums_l(:,13,i) = sums_wsus_ws_l(:)       ! w*u* 
-             sums_l(:,15,i) = sums_wsvs_ws_l(:)       ! w*v* 
-             sums_l(:,30,i) = sums_us2_ws_l(:)        ! u*2 
-             sums_l(:,31,i) = sums_vs2_ws_l(:)        ! v*2 
-             sums_l(:,32,i) = sums_ws2_ws_l(:)        ! w*2 
+             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u* 
+             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v* 
+             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2 
+             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2 
+             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2 
              sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        & 
-                              ( sums_us2_ws_l(:) + sums_vs2_ws_l(:) +       &
-                                sums_ws2_ws_l(:) )    ! e*
+                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
+                                sums_ws2_ws_l(:,i) )    ! e*
              DO  k = nzb, nzt
-                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * (  &
-                                                      sums_us2_ws_l(k) +     &
-                                                      sums_vs2_ws_l(k) +     &
-                                                      sums_ws2_ws_l(k) )
+                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
+                                                      sums_us2_ws_l(k,i) +  &
+                                                      sums_vs2_ws_l(k,i) +  &
+                                                      sums_ws2_ws_l(k,i) )
              ENDDO
           ENDDO
@@ -207,8 +206,8 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i) = sums_wspts_ws_l(:)      ! w*pt* from advec_s_ws
-             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:) ! w*sa*
+             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
+             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
              IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
-                                                   sums_wsqs_ws_l(:) !w*q*
+                                                   sums_wsqs_ws_l(:,i) !w*q*
           ENDDO
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 800)
+++ palm/trunk/SOURCE/modules.f90	(revision 801)
@@ -7,4 +7,8 @@
 ! particle arrays (particles, parrticles_temp) implemented as pointers,
 ! +particles_1, particles_2, sort_count
+!
+! Dimesion of sums_wsus_ws_l, ! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l,
+! sums_ws2_ws_l, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l increased.
+! for thread-safe summation in advec_ws.
 !
 ! Former revisions:
@@ -1439,11 +1443,11 @@
     REAL ::     u_max, v_max, w_max
     REAL, DIMENSION(:), ALLOCATABLE       ::  sums_divnew_l, sums_divold_l,   &
+                                              weight_substep, weight_pres
+    REAL, DIMENSION(:,:), ALLOCATABLE     ::  sums, sums_wsts_bc_l, ts_value, &
                                               sums_wsus_ws_l, sums_wsvs_ws_l, &
                                               sums_us2_ws_l, sums_vs2_ws_l,   &
                                               sums_ws2_ws_l,                  &
                                               sums_wspts_ws_l,                &
-                                              sums_wssas_ws_l,sums_wsqs_ws_l, &                 
-                                              weight_substep, weight_pres
-    REAL, DIMENSION(:,:), ALLOCATABLE     ::  sums, sums_wsts_bc_l, ts_value
+                                              sums_wssas_ws_l,sums_wsqs_ws_l
     REAL, DIMENSION(:,:,:), ALLOCATABLE   ::  hom_sum, rmask, spectrum_x, &
                                               spectrum_y, sums_l, sums_l_l, &
