SUBROUTINE diffusivities( var, var_reference )

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! OpenACC statements added + code changes required for GPU optimization,
! adjustment of mixing length to the Prandtl mixing length at first grid point
! above ground removed
! 
! Former revisions:
! -----------------
! $Id: diffusivities.f90 1015 2012-09-27 09:23:24Z raasch $
!
! 667 2010-12-23 12:06:00Z suehring/gryschka
! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
!
! 137 2007-11-28 08:50:10Z letzel
! Bugfix for summation of sums_l_l for flow_statistics
! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
!
! 97 2007-06-21 08:23:15Z raasch
! Adjustment of mixing length calculation for the ocean version.
! This is also a bugfix, because the height above the topography is now
! used instead of the height above level k=0.
! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
! use_pt_reference renamed use_reference
!
! 57 2007-03-09 12:05:41Z raasch
! Reference temperature pt_reference can be used in buoyancy term
!
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.24  2006/04/26 12:16:26  raasch
! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
!
! Revision 1.1  1997/09/19 07:41:10  raasch
! Initial revision
!
!
! Description:
! ------------
! Computation of the turbulent diffusion coefficients for momentum and heat 
! according to Prandtl-Kolmogorov
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE control_parameters
    USE grid_variables
    USE indices
    USE pegrid
    USE statistics

    IMPLICIT NONE

    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn

    REAL    ::  dvar_dz, l, ll, l_stable, sqrt_e, var_reference

    REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)


!
!-- Default thread number in case of one thread
    tn = 0

!
!-- Initialization for calculation of the mixing length profile
    sums_l_l = 0.0

!
!-- Compute the turbulent diffusion coefficient for momentum
    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn)
!$  tn = omp_get_thread_num()

!
!-- Data declerations for accelerators
    !$acc data present( dd2zu, e, km, kh, l_grid, l_wall, nzb_s_inner, rif, var )
    !$acc kernels

!
!-- Introduce an optional minimum tke
    IF ( e_min > 0.0 )  THEN
       !$OMP DO
       !$acc loop
       DO  i = nxlg, nxrg
          DO  j = nysg, nyng
             !$acc loop vector( 32 )
             DO  k = 1, nzt
                IF ( k > nzb_s_inner(j,i) )  THEN
                   e(k,j,i) = MAX( e(k,j,i), e_min )
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ENDIF

    !$OMP DO
    !$acc loop
    DO  i = nxlg, nxrg
       DO  j = nysg, nyng
          !$acc loop vector( 32 )
          DO  k = 1, nzt

             IF ( k > nzb_s_inner(j,i) )  THEN

                sqrt_e = SQRT( e(k,j,i) )
!
!--             Determine the mixing length
                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
                          ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
                IF ( dvar_dz > 0.0 ) THEN
                   IF ( use_reference )  THEN
                      l_stable = 0.76 * sqrt_e / &
                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
                   ELSE
                      l_stable = 0.76 * sqrt_e / &
                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
                   ENDIF
                ELSE
                   l_stable = l_grid(k)
                ENDIF
!
!--             Adjustment of the mixing length
                IF ( wall_adjustment )  THEN
                   l  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
                   ll = MIN( l_wall(k,j,i), l_grid(k) )
                ELSE
                   l  = MIN( l_grid(k), l_stable )
                   ll = l_grid(k)
                ENDIF

      !
      !--       Compute diffusion coefficients for momentum and heat
                km(k,j,i) = 0.1 * l * sqrt_e
                kh(k,j,i) = ( 1.0 + 2.0 * l / ll ) * km(k,j,i)

             ENDIF

#if ! defined( __openacc )
!
!++          Statistics still have to be realized for accelerators
!--          Summation for averaged profile (cf. flow_statistics)
             DO  sr = 0, statistic_regions
                sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)
             ENDDO
#endif

          ENDDO
       ENDDO
    ENDDO

#if ! defined( __openacc )
!
!++ Statistics still have to be realized for accelerators
    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
                                                  ! data output
#endif
    !$OMP END PARALLEL

!
!-- Set vertical boundary values (Neumann conditions both at bottom and top).
!-- Horizontal boundary conditions at vertical walls are not set because
!-- so far vertical walls require usage of a Prandtl-layer where the boundary
!-- values of the diffusivities are not needed
    !$OMP PARALLEL DO
    !$acc loop
    DO  i = nxlg, nxrg
       DO  j = nysg, nyng
          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
          km(nzt+1,j,i)            = km(nzt,j,i)
          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
          kh(nzt+1,j,i)            = kh(nzt,j,i)
       ENDDO
    ENDDO

!
!-- Set Neumann boundary conditions at the outflow boundaries in case of
!-- non-cyclic lateral boundaries
    IF ( outflow_l )  THEN
       km(:,:,nxl-1) = km(:,:,nxl)
       kh(:,:,nxl-1) = kh(:,:,nxl)
    ENDIF
    IF ( outflow_r )  THEN
       km(:,:,nxr+1) = km(:,:,nxr)
       kh(:,:,nxr+1) = kh(:,:,nxr)
    ENDIF
    IF ( outflow_s )  THEN
       km(:,nys-1,:) = km(:,nys,:)
       kh(:,nys-1,:) = kh(:,nys,:)
    ENDIF
    IF ( outflow_n )  THEN
       km(:,nyn+1,:) = km(:,nyn,:)
       kh(:,nyn+1,:) = kh(:,nyn,:)
    ENDIF

    !$acc end kernels
    !$acc end data

 END SUBROUTINE diffusivities