Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 1002)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -56,5 +56,5 @@
 
     PRIVATE
-    PUBLIC diffusion_v
+    PUBLIC diffusion_v, diffusion_v_acc
 
     INTERFACE diffusion_v
@@ -62,4 +62,8 @@
        MODULE PROCEDURE diffusion_v_ij
     END INTERFACE diffusion_v
+
+    INTERFACE diffusion_v_acc
+       MODULE PROCEDURE diffusion_v_acc
+    END INTERFACE diffusion_v_acc
 
  CONTAINS
@@ -218,4 +222,188 @@
 
 !------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_v_acc
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  kmxm, kmxp, kmzm, kmzp
+
+       !$acc declare create ( vsus )
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
+
+!
+!--    First calculate horizontal momentum flux v'u' at vertical walls,
+!--    if neccessary
+       IF ( topography /= 'flat' )  THEN
+          CALL wall_fluxes_acc( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
+                                nzb_v_outer, wall_v )
+       ENDIF
+
+       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )   &
+       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )           &
+       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nysv, nyn
+!
+!--          Compute horizontal diffusion
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF ( k > nzb_v_outer(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
+                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
+                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
+                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
+                         &   ) * ddx                                           &
+                         & + 2.0 * (                                           &
+                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
+                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
+                         &         ) * ddy2
+                ENDIF
+             ENDDO
+
+!
+!--          Wall functions at the left and right walls, respectively
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND. &
+                    wall_v(j,i) /= 0.0 )  THEN
+
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
+                   
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                 + 2.0 * (                                     &
+                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
+                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
+                                         ) * ddy2                              &
+                                 + (   fxp(j,i) * (                            &
+                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
+                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
+                                                  )                            &
+                                     - fxm(j,i) * (                            &
+                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
+                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
+                                                  )                            &
+                                     + wall_v(j,i) * vsus(k,j,i)               &
+                                   ) * ddx
+                ENDIF
+             ENDDO
+
+!
+!--          Compute vertical diffusion. In case of simulating a Prandtl
+!--          layer, index k starts at nzb_v_inner+2.
+             !$acc loop vector(32)
+             DO  k = 1, nzt_diff
+                IF ( k >= nzb_diff_v(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmzp = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                   kmzm = 0.25 * &
+                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
+                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
+                         &            )                                        &
+                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
+                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
+                         &            )                                        &
+                         &   ) * ddzw(k)
+                ENDIF
+             ENDDO
+
+          ENDDO
+       ENDDO
+
+!
+!--    Vertical diffusion at the first grid point above the surface,
+!--    if the momentum flux at the bottom is given by the Prandtl law
+!--    or if it is prescribed by the user.
+!--    Difference quotient of the momentum flux is not formed over
+!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
+!--    comparison with other (LES) modell showed that the values of
+!--    the momentum flux becomes too large in this case.
+!--    The term containing w(k-1,..) (see above equation) is removed here
+!--    because the vertical velocity is assumed to be zero at the surface.
+       IF ( use_surface_fluxes )  THEN
+
+          !$acc loop
+          DO  i = nxl, nxr
+             !$acc loop vector(32)
+             DO  j = nysv, nyn
+          
+                k = nzb_v_inner(j,i)+1
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
+                      &   ) * ddzw(k)                                        &
+                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
+                      &   + vsws(j,i)                                        &
+                      &   ) * ddzw(k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+!
+!--    Vertical diffusion at the first gridpoint below the top boundary,
+!--    if the momentum flux at the top is prescribed by the user
+       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
+
+          k = nzt
+
+          !$acc loop
+          DO  i = nxl, nxr
+             !$acc loop vector(32)
+             DO  j = nysv, nyn
+
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -vswst(j,i)                                      &
+                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
+                      &   ) * ddzw(k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+       !$acc end kernels
+
+    END SUBROUTINE diffusion_v_acc
+
+
+!------------------------------------------------------------------------------!
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
