Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 4467)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 4468)
@@ -20,9 +20,14 @@
 ! Current revisions:
 ! ------------------
-!
-!
+! 
+! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! - bugfix for last commit in openacc branch
+! - some loop bounds revised (only to be consistent with cache version)
+! - setting of nzb_max_l for advection of the w-component revised
+!
+! 4466 2020-03-20 16:14:41Z suehring
 ! - vector branch further optimized (linear dependencies along z removed and
 !   loops are splitted)
@@ -1822,5 +1827,5 @@
                                          )                                     &
                           ) * ddy                                              &
-                        + ( w(k,j,i) * rho_air_zw(k) *                         &
+                        + ( w(k,j,i)   * rho_air_zw(k) *                       &
                                          ( ibit6 + ibit7 + ibit8 )             &
                           - w(k-1,j,i) * rho_air_zw(k-1) *                     &
@@ -1860,5 +1865,5 @@
           div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                      &
                         + ( v(k,j+1,i) - v(k,j,i) ) * ddy                      &
-                        + ( w(k,j,i) * rho_air_zw(k)                           &
+                        + ( w(k,j,i)   * rho_air_zw(k)                         &
                           - w(k-1,j,i) * rho_air_zw(k-1)                       &
                                                   ) * drho_air(k) * ddzw(k)
@@ -3188,5 +3193,5 @@
            ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
            ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
-          nzb_max_l = nzt
+          nzb_max_l = nzt - 1
        ELSE
           nzb_max_l = nzb_max
@@ -3705,5 +3710,5 @@
        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
        INTEGER(iwp) ::  sk_num    !< integer identifier, used for assign fluxes to the correct dimension in the analysis array
-       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
+       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread (is always zero here)
 
        INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::   &
@@ -3755,10 +3760,17 @@
        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side of the grid box
        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top of the grid box
-
-       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
-       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
-
-       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
-       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
+#ifndef _OPENACC
+       REAL(wp), DIMENSION(nzb+1:nzt)         ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
+       REAL(wp), DIMENSION(nzb+1:nzt)         ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
+#endif
+#ifdef _OPENACC
+       REAL(wp), DIMENSION(nzb+1:nzt)         ::  diss_l            !< discretized artificial dissipation at leftward-side
+       REAL(wp), DIMENSION(nzb+1:nzt)         ::  diss_s            !< discretized artificial dissipation at southward-side
+       REAL(wp), DIMENSION(nzb+1:nzt)         ::  flux_l            !< discretized 6th-order flux at leftward-side
+       REAL(wp), DIMENSION(nzb+1:nzt)         ::  flux_s            !< discretized 6th-order flux at southward-side
+#endif
+
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
+       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
 
        CALL cpu_log( log_point_s(49), 'advec_s_ws', 'start' )
@@ -3795,7 +3807,7 @@
        !$ACC PRIVATE(ibit3_s, ibit4_s, ibit5_s) &
        !$ACC PRIVATE(ibit6, ibit7, ibit8) &
-       !$ACC PRIVATE(flux_r, diss_r) &
-       !$ACC PRIVATE(flux_n, diss_n) &
-       !$ACC PRIVATE(swap_diss_y_local, swap_flux_y_local) &
+       !$ACC PRIVATE(nzb_max_l) &
+       !$ACC PRIVATE(diss_l, diss_r, flux_l, flux_r) &
+       !$ACC PRIVATE(diss_n, diss_s, flux_s, flux_n) &
        !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
        !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s) &
@@ -3838,5 +3850,5 @@
 
                    u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
-                   swap_flux_x_local(k,j,tn) = u_comp * (                      &
+                   swap_flux_x_local(k,j)    = u_comp * (                      &
                                                ( 37.0_wp * ibit2 * adv_sca_5   &
                                             +     7.0_wp * ibit1 * adv_sca_3   &
@@ -3853,5 +3865,5 @@
                                                         )
 
-                   swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (              &
+                   swap_diss_x_local(k,j)    = -ABS( u_comp ) * (              &
                                                 ( 10.0_wp * ibit2 * adv_sca_5  &
                                              +     3.0_wp * ibit1 * adv_sca_3  &
@@ -3873,5 +3885,5 @@
 
                    u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
-                   swap_flux_x_local(k,j,tn) = u_comp * (                      &
+                   swap_flux_x_local(k,j)    = u_comp * (                      &
                                         37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) )&
                                       -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )&
@@ -3879,5 +3891,5 @@
                                                         ) * adv_sca_5
 
-                   swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (              &
+                   swap_diss_x_local(k,j)    = -ABS( u_comp ) * (              &
                                         10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )&
                                       -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )&
@@ -3900,5 +3912,5 @@
 
                    v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
-                   swap_flux_y_local(k,tn) = v_comp *         (                &
+                   swap_flux_y_local(k)    = v_comp *         (                &
                                                 ( 37.0_wp * ibit5 * adv_sca_5  &
                                              +     7.0_wp * ibit4 * adv_sca_3  &
@@ -3915,5 +3927,5 @@
                                                               )
 
-                   swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                &
+                   swap_diss_y_local(k)    = -ABS( v_comp ) * (                &
                                                 ( 10.0_wp * ibit5 * adv_sca_5  &
                                              +     3.0_wp * ibit4 * adv_sca_3  &
@@ -3936,10 +3948,10 @@
 
                    v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
-                   swap_flux_y_local(k,tn) = v_comp * (                        &
+                   swap_flux_y_local(k)    = v_comp * (                        &
                                      37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )   &
                                    -  8.0_wp * ( 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,tn) = -ABS( v_comp ) * (                &
+                   swap_diss_y_local(k)    = -ABS( v_comp ) * (                &
                                      10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )   &
                                    -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )   &
@@ -4000,6 +4012,6 @@
                 ibit0_l = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
 
-                u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k)
-                swap_flux_x_local(k,j,tn) = u_comp_l * (                       &
+                u_comp_l                  = u(k,j,i) - u_gtrans + u_stokes_zu(k)
+                flux_l(k)                 = u_comp_l * (                       &
                                              ( 37.0_wp * ibit2_l * adv_sca_5   &
                                           +     7.0_wp * ibit1_l * adv_sca_3   &
@@ -4016,5 +4028,5 @@
                                                        )
 
-                swap_diss_x_local(k,j,tn) = -ABS( u_comp_l ) * (               &
+                diss_l(k)                 = -ABS( u_comp_l ) * (               &
                                             ( 10.0_wp * ibit2_l * adv_sca_5    &
                                          +     3.0_wp * ibit1_l * adv_sca_3    &
@@ -4025,9 +4037,12 @@
                                          +              ibit1_l * adv_sca_3    &
                                             ) *                                &
-                                        ( sk(k,j,i+1) - sk(k,j,i-2)  )         &
+                                         ( sk(k,j,i+1) - sk(k,j,i-2) )         &
                                      +      (           ibit2_l * adv_sca_5    &
                                             ) *                                &
                                          ( sk(k,j,i+2) - sk(k,j,i-3) )         &
                                                                )
+#else
+                flux_l(k) = swap_flux_x_local(k,j)
+                diss_l(k) = swap_diss_x_local(k,j)
 #endif
                 ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
@@ -4073,10 +4088,10 @@
 
                 v_comp_s                = v(k,j,i) - v_gtrans + v_stokes_zu(k)
-                swap_flux_y_local(k,tn) = v_comp_s * (                         &
+                flux_s(k)               = v_comp_s * (                         &
                                              ( 37.0_wp * ibit5_s * adv_sca_5   &
                                           +     7.0_wp * ibit4_s * adv_sca_3   &
                                           +              ibit3_s * adv_sca_1   &
                                              ) *                               &
-                                         ( sk(k,j,i)  + sk(k,j-1,i)     )      &
+                                         ( sk(k,j,i)   + sk(k,j-1,i)     )     &
                                        -     (  8.0_wp * ibit5_s * adv_sca_5   &
                                           +              ibit4_s * adv_sca_3   &
@@ -4085,8 +4100,8 @@
                                       +      (           ibit5_s * adv_sca_5   &
                                              ) *                               &
-                                        ( sk(k,j+2,i) + sk(k,j-3,i)     )      &
+                                         ( sk(k,j+2,i) + sk(k,j-3,i)    )      &
                                                      )
 
-                swap_flux_y_local(k,tn) = -ABS( v_comp_s ) * (                 &
+                diss_s(k)               = -ABS( v_comp_s ) * (                 &
                                              ( 10.0_wp * ibit5_s * adv_sca_5   &
                                           +     3.0_wp * ibit4_s * adv_sca_3   &
@@ -4102,4 +4117,7 @@
                                           ( sk(k,j+2,i) - sk(k,j-3,i)    )     &
                                                              )
+#else
+                flux_s(k) = swap_flux_y_local(k)
+                diss_s(k) = swap_diss_y_local(k)
 #endif
              ENDDO
@@ -4123,5 +4141,5 @@
 !--             Recompute the left fluxes.
                 u_comp_l                  = u(k,j,i) - u_gtrans + u_stokes_zu(k)
-                swap_flux_x_local(k,j,tn) = u_comp_l * (                       &
+                flux_l(k)                 = u_comp_l * (                       &
                                       37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) )  &
                                     -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
@@ -4129,9 +4147,13 @@
                                                        ) * adv_sca_5
 
-                swap_diss_x_local(k,j,tn) = -ABS( u_comp_l ) * (               &
+                diss_l(k)                 = -ABS( u_comp_l ) * (               &
                                       10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )  &
                                     -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
                                     +           ( sk(k,j,i+2) - sk(k,j,i-3) )  &
                                                                ) * adv_sca_5
+#else
+                flux_l(k) = swap_flux_x_local(k,j)
+                diss_l(k) = swap_diss_x_local(k,j)
+
 #endif
 
@@ -4149,14 +4171,17 @@
 !--             Recompute the south fluxes.
                 v_comp_s                = v(k,j,i) - v_gtrans + v_stokes_zu(k)
-                swap_flux_y_local(k,tn) = v_comp_s * (                         &
+                flux_s(k)               = v_comp_s * (                         &
                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )    &
                                   -  8.0_wp * ( 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,tn) = -ABS( v_comp_s ) * (                 &
+                diss_s(k)               = -ABS( v_comp_s ) * (                 &
                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )    &
                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
                                   +             sk(k,j+2,i) - sk(k,j-3,i)      &
                                                              ) * adv_sca_5
+#else
+                flux_s(k) = swap_flux_y_local(k)
+                diss_s(k) = swap_diss_y_local(k)
 #endif
              ENDDO
@@ -4304,6 +4329,6 @@
              DO  k = nzb+1, nzb_max_l
 
-                flux_d    = flux_t(k-1)
-                diss_d    = diss_t(k-1)
+                flux_d = flux_t(k-1)
+                diss_d = diss_t(k-1)
 
                 ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
@@ -4336,5 +4361,5 @@
                                                )                               &
                                 ) * ddy                                        &
-                              + ( w(k,j,i) * rho_air_zw(k) *                   &
+                              + ( w(k,j,i)   * rho_air_zw(k)   *               &
                                                ( ibit6 + ibit7 + ibit8 )       &
                                 - w(k-1,j,i) * rho_air_zw(k-1) *               &
@@ -4347,19 +4372,17 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                  &
-                         ( 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        &
-                       + ( ( flux_t(k) + diss_t(k) ) -                         &
-                           ( flux_d    + diss_d    )                           &
-                                                                )              &
-                          * drho_air(k) * ddzw(k) &
+                         ( flux_r(k) + diss_r(k) -                             &
+                           flux_l(k) - diss_l(k) ) * ddx                       &
+                       + ( flux_n(k) + diss_n(k) -                             &
+                           flux_s(k) - diss_s(k) ) * ddy                       &
+                       + ( flux_t(k) + diss_t(k) -                             &
+                           flux_d    - diss_d    ) * drho_air(k) * ddzw(k)     &
                                             ) + sk(k,j,i) * div
 
 #ifndef _OPENACC
-                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)
+                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)
 #endif
 
@@ -4368,6 +4391,6 @@
              DO  k = nzb_max_l+1, nzt
 
-                flux_d    = flux_t(k-1)
-                diss_d    = diss_t(k-1)
+                flux_d = flux_t(k-1)
+                diss_d = diss_t(k-1)
 !
 !--             Calculate the divergence of the velocity field. A respective
@@ -4376,23 +4399,22 @@
                 div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                &
                               + ( v(k,j+1,i) - v(k,j,i) ) * ddy                &
-                              + ( w(k,j,i) * rho_air_zw(k)                     &
+                              + ( w(k,j,i)   * rho_air_zw(k)                   &
                                 - w(k-1,j,i) * rho_air_zw(k-1)                 &
                                                         ) * drho_air(k) * ddzw(k)
 
                 tend(k,j,i) = tend(k,j,i) - (                                  &
-                         ( 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        &
-                       + ( ( flux_t(k) + diss_t(k) ) -                         &
-                           ( flux_d    + diss_d    )            )              &
-                           * drho_air(k) * ddzw(k)                             &
+                         ( flux_r(k) + diss_r(k) -                             &
+                           flux_l(k) - diss_l(k) ) * ddx                       &
+                       + ( flux_n(k) + diss_n(k) -                             &
+                           flux_s(k) - diss_s(k) ) * ddy                       &
+                       + ( flux_t(k) + diss_t(k) -                             &
+                           flux_d    - diss_d    ) * drho_air(k) * ddzw(k)     &
                                             ) + sk(k,j,i) * div
 
 #ifndef _OPENACC
-                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)
+                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)
 #endif
             ENDDO
@@ -4569,4 +4591,5 @@
        !$ACC PRIVATE(ibit0_l, ibit1_l, ibit2_l) &
        !$ACC PRIVATE(ibit3_s, ibit4_s, ibit5_s) &
+       !$ACC PRIVATE(nzb_max_l) &
        !$ACC PRIVATE(ibit6, ibit7, ibit8) &
        !$ACC PRIVATE(flux_r, diss_r) &
@@ -4931,5 +4954,5 @@
              diss_t(nzb) = 0.0_wp
              w_comp(nzb) = 0.0_wp
-             DO  k = nzb+1, nzb+2
+             DO  k = nzb+1, nzb+1
 !
 !--             k index has to be modified near bottom and top, else array
@@ -4975,5 +4998,5 @@
              ENDDO
 
-             DO  k = nzb+3, nzt-2
+             DO  k = nzb+2, nzt-2
 
                 ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
@@ -5282,4 +5305,5 @@
        !$ACC PRIVATE(ibit9_l, ibit10_l, ibit11_l) &
        !$ACC PRIVATE(ibit12_s, ibit13_s, ibit14_s) &
+       !$ACC PRIVATE(nzb_max_l) &
        !$ACC PRIVATE(flux_r, diss_r) &
        !$ACC PRIVATE(flux_n, diss_n) &
@@ -5637,5 +5661,5 @@
              diss_t(nzb) = 0.0_wp
              w_comp(nzb) = 0.0_wp
-             DO  k = nzb+1, nzb+2
+             DO  k = nzb+1, nzb+1
 !
 !--             k index has to be modified near bottom and top, else array
@@ -5681,5 +5705,5 @@
              ENDDO
 
-             DO  k = nzb+3, nzt-2
+             DO  k = nzb+2, nzt-2
                 ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
                 ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
@@ -6004,4 +6028,5 @@
        !$ACC PRIVATE(ibit18_l, ibit19_l, ibit20_l) &
        !$ACC PRIVATE(ibit21_s, ibit22_s, ibit23_s) &
+       !$ACC PRIVATE(nzb_max_l) &
        !$ACC PRIVATE(flux_r, diss_r) &
        !$ACC PRIVATE(flux_n, diss_n) &
@@ -6027,5 +6052,5 @@
                  ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
                  ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
-                nzb_max_l = nzt
+                nzb_max_l = nzt - 1
              ELSE
                 nzb_max_l = nzb_max
