Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 732)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 736)
@@ -3,5 +3,9 @@
 !-----------------------------------------------------------------------------!
 ! Current revisions:
-! -----------------
+! ------------------
+! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
+! index where fluxes on left side have to be calculated explicitly is 
+! different on each thread. Furthermore the swapping of fluxes is now 
+! thread-safe by adding an additional dimension.
 !
 ! Former revisions:
@@ -43,7 +47,4 @@
 ! A further routine local_diss_ij is available (next weeks) and is used if a 
 ! control of dissipative fluxes is desired.
-! In case of vertical grid stretching the order of advection scheme is
-! degraded. This is also realized for the ocean where the stretching is 
-! downwards.
 !
 ! OUTSTANDING: - Dirichlet and Radiation boundary conditions for 
@@ -102,4 +103,5 @@
        USE control_parameters
        USE indices
+       USE pegrid
        USE statistics
 
@@ -166,13 +168,16 @@
           IF ( ws_scheme_mom )  THEN
 
-             ALLOCATE( flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),  &
-                       flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),  &
-                       diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt) )
-             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn),               &
-                       flux_l_v(nzb+1:nzt,nys:nyn),               &
-                       flux_l_w(nzb+1:nzt,nys:nyn),               &
-                       diss_l_u(nzb+1:nzt,nys:nyn),               &
-                       diss_l_v(nzb+1:nzt,nys:nyn),               &
-                       diss_l_w(nzb+1:nzt,nys:nyn) )
+             ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),            &
+                       flux_s_v(nzb+1:nzt,0:threads_per_task-1),            &
+                       flux_s_w(nzb+1:nzt,0:threads_per_task-1),            &
+                       diss_s_u(nzb+1:nzt,0:threads_per_task-1),            &
+                       diss_s_v(nzb+1:nzt,0:threads_per_task-1),            &
+                       diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
+             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                       flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                       flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                       diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                       diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                       diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
 
           ENDIF
@@ -180,21 +185,25 @@
           IF ( ws_scheme_sca )  THEN
 
-             ALLOCATE( flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), &
-                       diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt) ) 
-             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn),              &
-                       flux_l_e(nzb+1:nzt,nys:nyn),               &
-                       diss_l_pt(nzb+1:nzt,nys:nyn),              &
-                       diss_l_e(nzb+1:nzt,nys:nyn) )
+             ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
+                       flux_s_e(nzb+1:nzt,0:threads_per_task-1),            &
+                       diss_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
+                       diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 
+             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                       flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                       diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                       diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
 
              IF ( humidity .OR. passive_scalar )  THEN
-                ALLOCATE( flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt) )
-                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn),            &
-                          diss_l_q(nzb+1:nzt,nys:nyn) )
+                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),         &
+                          diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
 
              IF ( ocean )  THEN
-                ALLOCATE( flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt) )
-                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn),           &
-                          diss_l_sa(nzb+1:nzt,nys:nyn) )
+                ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1),         &
+                          diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
 
@@ -267,5 +276,5 @@
     SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
                               swap_diss_y_local, swap_flux_x_local, &
-                              swap_diss_x_local  )
+                              swap_diss_x_local, i_omp, tn )
 
        USE arrays_3d
@@ -274,9 +283,10 @@
        USE grid_variables
        USE indices
+       USE pegrid
        USE statistics
 
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
+       INTEGER ::  i, i_omp, j, k, tn
        LOGICAL ::  degraded_l, degraded_s
        REAL    ::  flux_d, diss_d, u_comp, v_comp
@@ -284,8 +294,9 @@
        REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
                                              flux_n, diss_n
-       REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &
-                                             swap_diss_y_local
-       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
-                                             swap_diss_x_local
+       REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local,  &
+                                                          swap_diss_y_local
+       REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::             &
+                                                          swap_flux_x_local,  &
+                                                          swap_diss_x_local
        CHARACTER (LEN = *), INTENT(IN)    :: sk_char
 
@@ -368,7 +379,7 @@
 !--                Compute leftside fluxes for the left boundary of PE domain
                    u_comp                 = u(k,j,i) - u_gtrans
-                   swap_flux_x_local(k,j) = u_comp * (                         &
+                   swap_flux_x_local(k,j,tn) = u_comp * (                         &
                                             sk(k,j,i) + sk(k,j,i-1) ) * 0.5
-                   swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
+                   swap_diss_x_local(k,j,tn) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
                                                       sk(k,j,i), sk(k,j,i-1),  &
                                                       sk(k,j,i-1), u_comp,     &
@@ -402,7 +413,7 @@
 !--                Compute southside fluxes for the south boundary of PE domain
                    v_comp               = v(k,j,i) - v_gtrans
-                   swap_flux_y_local(k) = v_comp *                             &
+                   swap_flux_y_local(k,tn) = v_comp *                             &
                                           ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 
-                   swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
+                   swap_diss_y_local(k,tn) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
                                                     sk(k,j,i), sk(k,j-1,i),    &
                                                     sk(k,j-1,i), v_comp,       &
@@ -477,14 +488,14 @@
 !       
 !--    Compute left- and southside fluxes of the respective PE bounds.        
-       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
+       IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
        
            DO  k = nzb_s_inner(j,i)+1, nzt
               v_comp               = v(k,j,i) - v_gtrans
-              swap_flux_y_local(k) = v_comp * (                               &
+              swap_flux_y_local(k,tn) = v_comp * (                               &
                                      37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
                                    -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
                                    +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
                                               ) * adv_sca_5
-              swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
+              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                       &
                                      10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
                                    -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
@@ -495,14 +506,14 @@
         ENDIF
         
-        IF ( i == nxl  .AND.  .NOT. degraded_l )  THEN
+        IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
         
            DO  k = nzb_s_inner(j,i)+1, nzt
               u_comp                 = u(k,j,i) - u_gtrans
-              swap_flux_x_local(k,j) = u_comp * (                             &
+              swap_flux_x_local(k,j,tn) = u_comp * (                             &
                                        37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
                                      -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
                                      +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
                                                 ) * adv_sca_5
-              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
+              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                     &
                                        10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
                                      -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
@@ -518,14 +529,14 @@
           
           tend(k,j,i) = tend(k,j,i) - (                                      &
-                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
-                            swap_diss_x_local(k,j) ) * ddx                   &
-                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   - &
-                            swap_diss_y_local(k)   ) * ddy                   &
+                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
+                            swap_diss_x_local(k,j,tn) ) * ddx                   &
+                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
+                            swap_diss_y_local(k,tn)   ) * ddy                   &
                                       )
 
-          swap_flux_y_local(k)   = flux_n(k)
-          swap_diss_y_local(k)   = diss_n(k)
-          swap_flux_x_local(k,j) = flux_r(k)
-          swap_diss_x_local(k,j) = diss_r(k)
+          swap_flux_y_local(k,tn)   = flux_n(k)
+          swap_diss_y_local(k,tn)   = diss_n(k)
+          swap_flux_x_local(k,j,tn) = flux_r(k)
+          swap_diss_x_local(k,j,tn) = diss_r(k)
           
        ENDDO
@@ -661,5 +672,5 @@
 ! Advection of u-component - Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE advec_u_ws_ij( i, j )
+    SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
 
        USE arrays_3d
@@ -672,9 +683,9 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
+       INTEGER ::  i, i_omp, j, k, tn
        LOGICAL ::  degraded_l, degraded_s
        REAL    ::  gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp
-       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_r, diss_r,        &
-                                       flux_n, diss_n, u_comp
+       REAL, DIMENSION(nzb:nzt+1) ::  flux_t, diss_t, flux_r, diss_r,        &
+                                      flux_n, diss_n, u_comp
                                        
        degraded_l = .FALSE.
@@ -696,5 +707,5 @@
                             -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
                   diss_r(k) = - ABS( u_comp(k) - gu ) * (                     &
-                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               & 
+                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
                             -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
               ENDDO
@@ -717,5 +728,5 @@
                             -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
                   diss_n(k) = - ABS( v_comp ) * (                             &
-                              3.0 * ( u(k,j+1,i) - u(k,j,i)     )             & 
+                              3.0 * ( u(k,j+1,i) - u(k,j,i)     )             &
                             -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
                ENDDO
@@ -735,6 +746,6 @@
 !--               Compute leftside fluxes for the left boundary of PE domain
                   u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
-                  flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
-                  diss_l_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
+                  flux_l_u(k,j,tn) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
+                  diss_l_u(k,j,tn) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
                                             u(k,j,i-1), u(k,j,i-1), u_comp(k),&
                                             0.25, ddx )
@@ -745,5 +756,5 @@
                             -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
                   diss_r(k) = - ABS( u_comp(k) -gu ) * (                      &
-                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               & 
+                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
                             -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
                ENDDO
@@ -757,5 +768,5 @@
                             -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
                   diss_n(k) = - ABS( v_comp ) * (                             &
-                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               & 
+                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
                             -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
                ENDDO 
@@ -764,8 +775,8 @@
                DO  k = nzb_u_inner(j,i)+1, nzt
 !               
-!--              Compute southside fluxes for the south boundary of PE domain           
+!--              Compute southside fluxes for the south boundary of PE domain
                   v_comp      = v(k,j,i) + v(k,j,i-1) - gv
-                  flux_s_u(k) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25 
-                  diss_s_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
+                  flux_s_u(k,tn) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25 
+                  diss_s_u(k,tn) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
                                           u(k,j-1,i), u(k,j-1,i), v_comp,     &
                                           0.25, ddy )
@@ -776,5 +787,5 @@
                             -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
                   diss_n(k) = - ABS( v_comp ) * (                             &
-                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               & 
+                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
                             -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
                ENDDO  
@@ -819,5 +830,5 @@
        ELSE
 !       
-!--       Compute the fifth order fluxes for the interior of PE domain.                  
+!--       Compute the fifth order fluxes for the interior of PE domain.
           DO  k = nzb_u_inner(j,i)+1, nzt
              u_comp(k) = u(k,j,i+1) + u(k,j,i)
@@ -844,14 +855,14 @@
        ENDIF
 !       
-!--    Compute left- and southside fluxes for the respective boundary of PE      
+!--    Compute left- and southside fluxes for the respective boundary of PE
        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
        
           DO  k = nzb_u_inner(j,i)+1, nzt
              v_comp      = v(k,j,i) + v(k,j,i-1) - gv
-             flux_s_u(k) = v_comp * (                                         &
+             flux_s_u(k,tn) = v_comp * (                                         &
                            37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
                          -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
                          +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
-             diss_s_u(k) = - ABS(v_comp) * (                                  &
+             diss_s_u(k,tn) = - ABS(v_comp) * (                                  &
                            10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
                          -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
@@ -861,13 +872,13 @@
        ENDIF
        
-       IF ( i == nxlu .AND. ( .NOT. degraded_l ) )  THEN
+       IF ( i == i_omp .AND. ( .NOT. degraded_l ) )  THEN
        
           DO  k = nzb_u_inner(j,i)+1, nzt
              u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
-             flux_l_u(k,j) = u_comp_l * (                                     &
+             flux_l_u(k,j,tn) = u_comp_l * (                                     &
                              37.0 * ( u(k,j,i) + u(k,j,i-1)   )               &
                            -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )               &
                            +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
-             diss_l_u(k,j) = - ABS(u_comp_l) * (                              &
+             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                              &
                              10.0 * ( u(k,j,i) - u(k,j,i-1)   )               &
                            -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )               &
@@ -881,12 +892,12 @@
           tend(k,j,i) = tend(k,j,i) - (                                       &
                             ( flux_r(k) + diss_r(k)                           &
-                          -   flux_l_u(k,j) - diss_l_u(k,j) ) * ddx         &
+                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx           &
                           + ( flux_n(k) + diss_n(k)                           &
-                          -   flux_s_u(k) - diss_s_u(k)     ) * ddy  )
-
-           flux_l_u(k,j) = flux_r(k)
-           diss_l_u(k,j) = diss_r(k)
-           flux_s_u(k)   = flux_n(k)
-           diss_s_u(k)   = diss_n(k)
+                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy  )
+
+           flux_l_u(k,j,tn) = flux_r(k)
+           diss_l_u(k,j,tn) = diss_r(k)
+           flux_s_u(k,tn)   = flux_n(k)
+           diss_s_u(k,tn)   = diss_n(k)
 !
 !--        Statistical Evaluation of u'u'. The factor has to be applied for 
@@ -897,5 +908,5 @@
              / ( u_comp(k) - gu + 1.0E-20      )                              &
              +   diss_r(k) *                                                  &
-                 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                        & 
+                 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                        &
              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )                          &
              *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
@@ -998,8 +1009,8 @@
 
 
-!------------------------------------------------------------------------------!
+!-----------------------------------------------------------------------------!
 ! Advection of v-component - Call for grid point i,j
-!------------------------------------------------------------------------------!
-   SUBROUTINE advec_v_ws_ij( i, j )
+!-----------------------------------------------------------------------------!
+   SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
 
        USE arrays_3d
@@ -1012,5 +1023,5 @@
        IMPLICIT NONE
 
-       INTEGER  ::  i, j, k
+       INTEGER  ::  i, i_omp, j, k, tn
        LOGICAL  ::  degraded_l, degraded_s
        REAL     ::  gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp
@@ -1036,5 +1047,5 @@
                              -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
                    diss_r(k) = - ABS( u_comp ) * (                            &
-                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              & 
+                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
                              -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
                 ENDDO
@@ -1056,5 +1067,5 @@
                              -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
                    diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
-                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              & 
+                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
                              -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
                 ENDDO
@@ -1077,5 +1088,5 @@
                              -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
                    diss_r(k) = - ABS( u_comp ) * (                            &
-                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )              & 
+                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
                              -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
                 ENDDO
@@ -1084,6 +1095,6 @@
                 DO  k = nzb_v_inner(j,i)+1, nzt
                    u_comp        = u(k,j-1,i) + u(k,j,i) - gu
-                   flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
-                   diss_l_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
+                   flux_l_v(k,j,tn) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
+                   diss_l_v(k,j,tn) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
                                              v(k,j,i-1), v(k,j,i-1), u_comp,  &
                                              0.25, ddx )
@@ -1094,5 +1105,5 @@
                              -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
                    diss_r(k) = - ABS( u_comp ) * (                            &
-                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              & 
+                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
                              -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 
                 ENDDO
@@ -1102,6 +1113,6 @@
                 DO  k = nzb_v_inner(j,i)+1, nzt
                    v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
-                   flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
-                   diss_s_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
+                   flux_s_v(k,tn) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
+                   diss_s_v(k,tn) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
                                            v(k,j-1,i), v(k,j-1,i), v_comp(k), &
                                            0.25, ddy )
@@ -1112,5 +1123,5 @@
                              -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
                    diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
-                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              & 
+                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
                              -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
                 ENDDO
@@ -1181,12 +1192,13 @@
 !       
 !--    Compute left- and southside fluxes for the respective boundary       
-       IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
+       IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
+
           DO  k = nzb_v_inner(j,i)+1, nzt
              u_comp        = u(k,j-1,i) + u(k,j,i) - gu
-             flux_l_v(k,j) = u_comp * (                                       &
+             flux_l_v(k,j,tn) = u_comp * (                                       &
                              37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
                            -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
                            +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
-             diss_l_v(k,j) = - ABS( u_comp ) * (                              &
+             diss_l_v(k,j,tn) = - ABS( u_comp ) * (                              &
                              10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
                            -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
@@ -1196,13 +1208,13 @@
        ENDIF
        
-       IF ( j == nysv .AND. ( .NOT. degraded_s ) )  THEN
+       IF ( j == nysv  .AND.  .NOT. degraded_s )  THEN
        
           DO  k = nzb_v_inner(j,i)+1, nzt
              v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
-             flux_s_v(k) = v_comp_l * (                                       &
+             flux_s_v(k,tn) = v_comp_l * (                                       &
                            37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
                          -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
                          +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
-             diss_s_v(k) = - ABS( v_comp_l ) * (                              &
+             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                              &
                            10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
                          -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
@@ -1216,12 +1228,12 @@
           tend(k,j,i) = tend(k,j,i) - (                                       &
                          ( flux_r(k) + diss_r(k)                              &
-                       -   flux_l_v(k,j) - diss_l_v(k,j)   ) * ddx            &
+                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx            &
                        + ( flux_n(k) + diss_n(k)                              &
-                       -   flux_s_v(k) - diss_s_v(k)       ) * ddy     )
+                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy     )
         
-           flux_l_v(k,j) = flux_r(k)
-           diss_l_v(k,j) = diss_r(k)
-           flux_s_v(k)   = flux_n(k)
-           diss_s_v(k)   = diss_n(k)
+           flux_l_v(k,j,tn) = flux_r(k)
+           diss_l_v(k,j,tn) = diss_r(k)
+           flux_s_v(k,tn)   = flux_n(k)
+           diss_s_v(k,tn)   = diss_n(k)
 
 !
@@ -1333,5 +1345,5 @@
 ! Advection of w-component - Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE advec_w_ws_ij( i, j )
+    SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
 
        USE arrays_3d
@@ -1344,5 +1356,5 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k
+       INTEGER ::  i, i_omp, j, k, tn
        LOGICAL ::  degraded_l, degraded_s
        REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
@@ -1355,5 +1367,6 @@
        gu = 2.0 * u_gtrans
        gv = 2.0 * v_gtrans
-       
+
+              
        IF ( boundary_flags(j,i) /= 0 )  THEN
 !       
@@ -1425,6 +1438,6 @@
                   
                   u_comp        = u(k+1,j,i) + u(k,j,i) - gu
-                  flux_l_w(k,j) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
-                  diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
+                  flux_l_w(k,j,tn) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
+                  diss_l_w(k,j,tn) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
                                             w(k,j,i-1), w(k,j,i-1), u_comp,   &
                                             0.25, ddx )
@@ -1455,6 +1468,6 @@
 !--              Compute southside fluxes for the south boundary of PE domain           
                   v_comp      = v(k+1,j,i) + v(k,j,i) - gv
-                  flux_s_w(k) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
-                  diss_s_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
+                  flux_s_w(k,tn) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
+                  diss_s_w(k,tn) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
                                           w(k,j-1,i), w(k,j-1,i), v_comp,     &
                                           0.25, ddy )                 
@@ -1526,13 +1539,13 @@
 !       
 !--    Compute left- and southside fluxes for the respective boundary      
-       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
+       IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
        
           DO  k = nzb_w_inner(j,i)+1, nzt
              v_comp      = v(k+1,j,i) + v(k,j,i) - gv
-             flux_s_w(k) = v_comp * (                                         &
+             flux_s_w(k,tn) = v_comp * (                                         &
                            37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
                          -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
                          +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
-             diss_s_w(k) = - ABS( v_comp ) * (                                &
+             diss_s_w(k,tn) = - ABS( v_comp ) * (                                &
                            10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
                          -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
@@ -1542,13 +1555,13 @@
        ENDIF
        
-       IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN
-       
+       IF ( i == i_omp  .AND.   .NOT. degraded_l ) THEN
+          
           DO  k = nzb_w_inner(j,i)+1, nzt
             u_comp        = u(k+1,j,i) + u(k,j,i) - gu
-            flux_l_w(k,j) = u_comp * (                                        &
+            flux_l_w(k,j,tn) = u_comp * (                                        &
                             37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
                           -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
                           +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
-            diss_l_w(k,j) = - ABS( u_comp ) * (                               &
+            diss_l_w(k,j,tn) = - ABS( u_comp ) * (                               &
                             10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
                           -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
@@ -1557,4 +1570,5 @@
           
        ENDIF
+
 !       
 !--    Now compute the tendency terms for the horizontal parts.         
@@ -1562,12 +1576,12 @@
           tend(k,j,i) = tend(k,j,i) - (                                       &
                       ( flux_r(k) + diss_r(k)                                 &
-                    -   flux_l_w(k,j) - diss_l_w(k,j)   ) * ddx               &
+                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx               &
                     + ( flux_n(k) + diss_n(k)                                 &
-                    -   flux_s_w(k) - diss_s_w(k)       ) * ddy  )
-
-          flux_l_w(k,j) = flux_r(k)
-          diss_l_w(k,j) = diss_r(k)
-          flux_s_w(k) = flux_n(k)
-          diss_s_w(k) = diss_n(k)
+                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy  )
+
+          flux_l_w(k,j,tn) = flux_r(k)
+          diss_l_w(k,j,tn) = diss_r(k)
+          flux_s_w(k,tn) = flux_n(k)
+          diss_s_w(k,tn) = diss_n(k)
        ENDDO
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 732)
+++ palm/trunk/SOURCE/modules.f90	(revision 736)
@@ -5,5 +5,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Dimension of fluxes needed for WS-scheme increased. 
 !
 ! Former revisions:
@@ -246,14 +246,12 @@
           km_damp_x, km_damp_y, lad, l_grid, pt_init, q_init, rdf, sa_init,    &
           ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs,  &
-          zu, zw, flux_s_u, flux_s_v, flux_s_w, diss_s_u, diss_s_v, diss_s_w,  &
-          flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q,        &
-          flux_s_sa, diss_s_sa
+          zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
           c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg,                  &
           mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, ts, us, z0,    &
-          flux_l_u, flux_l_v, flux_l_w, diss_l_u, diss_l_v, diss_l_w,          &
-          flux_l_pt, diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q,        &
-          flux_l_sa, diss_l_sa, total_2d_o, total_2d_a
+          flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q,        &
+          flux_s_sa, diss_s_sa, flux_s_u, flux_s_v, flux_s_w, diss_s_u,        &
+          diss_s_v, diss_s_w, total_2d_o, total_2d_a
 
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
@@ -270,5 +268,7 @@
           canopy_heat_flux, cdc, d, diss, lad_s, lad_u, lad_v, lad_w, lai,     &
           l_wall, p_loc, sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,    &
-          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s, flux_l_pt,          &
+          diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q, flux_l_sa,        &
+          diss_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_l_u, diss_l_v, diss_l_w
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 732)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 736)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: determination of first thread index i for WS-scheme
 !
 ! Former revisions:
@@ -147,5 +147,5 @@
 
     CHARACTER (LEN=9) ::  time_to_string
-    INTEGER ::  i, j, k
+    INTEGER ::  i, i_omp_start, j, k, tn = 0
     REAL    ::  sat, sbt
 
@@ -172,4 +172,5 @@
 !
 !-- u-tendency terms with no communication
+    i_omp_start = nxlu
     DO  i = nxlu, nxr
        DO  j = nys, nyn
@@ -179,5 +180,5 @@
              tend(:,j,i) = 0.0
              IF ( ws_scheme_mom )  THEN
-                 CALL advec_u_ws( i, j )
+                 CALL advec_u_ws( i, j, i_omp_start, tn )
              ELSE
                  CALL advec_u_pw( i, j )
@@ -257,4 +258,5 @@
 !
 !-- v-tendency terms with no communication
+    i_omp_start = nxl
     DO  i = nxl, nxr
        DO  j = nysv, nyn
@@ -264,5 +266,5 @@
              tend(:,j,i) = 0.0
              IF ( ws_scheme_mom )  THEN
-                 CALL advec_v_ws( i, j )
+                 CALL advec_v_ws( i, j, i_omp_start, tn )
              ELSE 
                  CALL advec_v_pw( i, j )
@@ -348,5 +350,5 @@
              tend(:,j,i) = 0.0
              IF ( ws_scheme_mom )  THEN
-                 CALL advec_w_ws( i, j )
+                 CALL advec_w_ws( i, j, i_omp_start, tn )
              ELSE  
                  CALL advec_w_pw( i, j )
@@ -453,5 +455,5 @@
                 IF ( ws_scheme_sca )  THEN
                    CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
-                                 diss_s_pt, flux_l_pt, diss_l_pt )
+                                 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
                 ELSE
                     CALL advec_s_pw( i, j, pt )
@@ -574,5 +576,5 @@
                    IF ( ws_scheme_sca )  THEN
                        CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
-                                    diss_s_sa, flux_l_sa, diss_l_sa )
+                                    diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
                    ELSE
                        CALL advec_s_pw( i, j, sa )
@@ -673,5 +675,5 @@
                    IF ( ws_scheme_sca )  THEN
                        CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
-                                   diss_s_q, flux_l_q, diss_l_q )
+                                   diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
                    ELSE
                        CALL advec_s_pw( i, j, q )
@@ -813,5 +815,5 @@
                       IF ( ws_scheme_sca )  THEN
                           CALL advec_s_ws( i, j, e, 'e', flux_s_e, & 
-                                      diss_s_e, flux_l_e, diss_l_e )
+                                      diss_s_e, flux_l_e, diss_l_e, i_omp_start, tn )
                       ELSE
                           CALL advec_s_pw( i, j, e )
@@ -916,5 +918,6 @@
 
     CHARACTER (LEN=9) ::  time_to_string
-    INTEGER ::  i, j, k
+    INTEGER ::  i, i_omp_start, j, k, omp_get_thread_num, tn = 0
+    LOGICAL ::  loop_start
 
 
@@ -934,10 +937,21 @@
          intermediate_timestep_count == 1 )  CALL ws_statistics
 
-
 !
 !-- Loop over all prognostic equations
-!$OMP PARALLEL private (i,j,k)
+!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
+
+!$  tn = omp_get_thread_num() 
+    loop_start = .TRUE.
 !$OMP DO
     DO  i = nxl, nxr
+
+!
+!--    Store the first loop index. It differs for each thread and is required
+!--    later in advec_ws
+       IF ( loop_start )  THEN
+          loop_start  = .FALSE.
+          i_omp_start = i 
+       ENDIF
+ 
        DO  j = nys, nyn
 !
@@ -948,10 +962,14 @@
              IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( ws_scheme_mom )  THEN
-    !               CALL local_diss( i, j, u)    ! dissipation control
-                   CALL advec_u_ws( i, j )
+                   IF ( outflow_l .AND. i_omp_start == nxl )  THEN 
+!                     CALL local_diss( i, j, u)    ! dissipation control
+                      CALL advec_u_ws( i, j, i_omp_start + 1, tn )
+                   ELSE
+                      CALL advec_u_ws( i, j, i_omp_start, tn )
+                   ENDIF
                 ELSE 
                    CALL advec_u_pw( i, j )
                 ENDIF 
-            ELSE
+             ELSE
                 CALL advec_u_up( i, j )
              ENDIF
@@ -1017,5 +1035,5 @@
                 IF ( ws_scheme_mom )  THEN
                  !   CALL local_diss( i, j, v)
-                    CALL advec_v_ws( i, j )
+                    CALL advec_v_ws( i, j, i_omp_start, tn )
                 ELSE 
                     CALL advec_v_pw( i, j )
@@ -1081,5 +1099,5 @@
              IF ( ws_scheme_mom )  THEN
              !   CALL local_diss( i, j, w)
-                CALL advec_w_ws( i, j )
+                CALL advec_w_ws( i, j, i_omp_start, tn )
              ELSE 
                 CALL advec_w_pw( i, j )
@@ -1145,5 +1163,5 @@
        !            CALL local_diss( i, j, pt )
                    CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
-                             diss_s_pt, flux_l_pt, diss_l_pt )
+                             diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
                 ELSE
                    CALL advec_s_pw( i, j, pt )
@@ -1227,5 +1245,5 @@
             !        CALL local_diss( i, j, sa )
                     CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
-                                diss_s_sa, flux_l_sa, diss_l_sa  )
+                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
                 ELSE 
                     CALL advec_s_pw( i, j, sa )
@@ -1285,5 +1303,5 @@
           !         CALL local_diss( i, j, q )
                    CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 
-                                diss_s_q, flux_l_q, diss_l_q )
+                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
                 ELSE 
                    CALL advec_s_pw( i, j, q )
@@ -1361,5 +1379,5 @@
                  !    CALL local_diss( i, j, e )
                      CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
-                                diss_s_e, flux_l_e, diss_l_e  )
+                                diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn )
                  ELSE
                      CALL advec_s_pw( i, j, e )
