Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 1132)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 1171)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! use_reference-case activated in accelerator version
 !
 ! Former revisions:
@@ -317,74 +317,72 @@
 !--    runtime errors occur with -C hopt on NEC
        IF ( use_reference )  THEN
-          STOP '+++ use_reference not implemented in diffusion_e_acc'
-!          DO  i = nxl, nxr
-!             DO  j = nys, nyn
-!                DO  k = nzb_s_inner(j,i)+1, nzt
-!
-!--                Calculate the mixing length (for dissipation)
-!                   dvar_dz = atmos_ocean_sign * &
-!                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-!                   IF ( dvar_dz > 0.0 ) THEN
-!                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-!                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
-!                   ELSE
-!                      l_stable = l_grid(k)
-!                   ENDIF
-!
-!--                Adjustment of the mixing length
-!                   IF ( wall_adjustment )  THEN
-!                      l(k,j)  = MIN( wall_adjustment_factor *          &
-!                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
-!                                     l_grid(k), l_stable )
-!                      ll(k,j) = MIN( wall_adjustment_factor *          &
-!                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
-!                                     l_grid(k) )
-!                   ELSE
-!                      l(k,j)  = MIN( l_grid(k), l_stable )
-!                      ll(k,j) = l_grid(k)
-!                   ENDIF
-!
-!                ENDDO
-!             ENDDO
-!
-!
-!--          Calculate the tendency terms
-!             DO  j = nys, nyn
-!                DO  k = nzb_s_inner(j,i)+1, nzt
-!
-!                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
-!                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
-!
-!                    tend(k,j,i) = tend(k,j,i)                                  &
-!                                        + (                                    &
-!                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
-!                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
-!                                          ) * ddx2                             &
-!                                        + (                                    &
-!                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
-!                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
-!                                          ) * ddy2                             &
-!                                        + (                                    &
-!               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
-!             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-!                                          ) * ddzw(k)                          &
-!                             - dissipation(k,j)
-!
-!                ENDDO
-!             ENDDO
-!
-!
-!--          Store dissipation if needed for calculating the sgs particle
-!--          velocities
-!             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
-!                DO  j = nys, nyn
-!                   DO  k = nzb_s_inner(j,i)+1, nzt
-!                      diss(k,j,i) = dissipation(k,j)
-!                   ENDDO
-!                ENDDO
-!             ENDIF
-!
-!          ENDDO
-!
+
+          !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
+          !$acc         present( nzb_s_inner, rif, tend, var, zu, zw )
+          !$acc loop
+          DO  i = i_left, i_right
+             DO  j = j_south, j_north
+                !$acc loop vector( 32 )
+                DO  k = 1, nzt
+
+                   IF ( k > nzb_s_inner(j,i) )  THEN
+!
+!--                   Calculate the mixing length (for dissipation)
+                      dvar_dz = atmos_ocean_sign * &
+                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                      IF ( dvar_dz > 0.0 ) THEN
+                         l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                    SQRT( g / var_reference * dvar_dz ) + 1E-5
+                      ELSE
+                         l_stable = l_grid(k)
+                      ENDIF
+!
+!--                   Adjustment of the mixing length
+                      IF ( wall_adjustment )  THEN
+                         l  = MIN( wall_adjustment_factor *          &
+                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                   l_grid(k), l_stable )
+                         ll = MIN( wall_adjustment_factor *          &
+                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                   l_grid(k) )
+                      ELSE
+                         l  = MIN( l_grid(k), l_stable )
+                         ll = l_grid(k)
+                      ENDIF
+!
+!--                   Calculate the tendency terms
+                      dissipation = ( 0.19 + 0.74 * l / ll ) * &
+                                    e(k,j,i) * SQRT( e(k,j,i) ) / l
+
+                      tend(k,j,i) = tend(k,j,i)                                  &
+                                        + (                                    &
+                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
+                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
+                                          ) * ddx2                             &
+                                        + (                                    &
+                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
+                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
+                                          ) * ddy2                             &
+                                        + (                                    &
+               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
+                                          ) * ddzw(k)                          &
+                                  - dissipation
+
+!
+!--                   Store dissipation if needed for calculating the sgs particle
+!--                   velocities
+                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
+                           turbulence )  THEN
+                         diss(k,j,i) = dissipation
+                      ENDIF
+
+                   ENDIF
+
+                ENDDO
+             ENDDO
+          ENDDO
+          !$acc end kernels
+
        ELSE
 
@@ -439,5 +437,5 @@
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
                                           ) * ddzw(k)                          &
-                             - dissipation
+                                  - dissipation
 
 !
