Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3550)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3550 2018-11-21 16:01:01Z gronemeier
 ! - calculate diss production same in vector and cache optimization
 ! - move boundary condition for e and diss to boundary_conds
@@ -1619,4 +1622,6 @@
        ENDDO  !k loop
 
+       !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
+
     ENDIF  !LES or RANS mode
 
@@ -1624,4 +1629,7 @@
 !-- Set lateral boundary conditions for l_wall
     CALL exchange_horiz( l_wall, nbgp )
+
+    !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
+    !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
 
     CONTAINS
@@ -1777,4 +1785,6 @@
 !--    Default surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu)
        DO  m = 1, surf_def_h(0)%ns
 
@@ -1811,4 +1821,6 @@
 !--    Default surfaces, downward-facing surfaces
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km)
        DO  m = 1, surf_def_h(1)%ns
 
@@ -1842,4 +1854,6 @@
 !--    Natural surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu)
        DO  m = 1, surf_lsm_h%ns
 
@@ -1876,4 +1890,6 @@
 !--    Urban surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu)
        DO  m = 1, surf_usm_h%ns
 
@@ -1970,5 +1986,7 @@
              CALL advec_s_up( e )
           ELSE
+             !$ACC KERNELS PRESENT(tend)
              tend = 0.0_wp
+             !$ACC END KERNELS
              IF ( timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( ws_scheme_sca )  THEN
@@ -2006,4 +2024,8 @@
 !--    reasons in the course of the integration. In such cases the old TKE
 !--    value is reduced by 90%.
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(e, tend, te_m, wall_flags_0) &
+       !$ACC PRESENT(tsc(3:3)) &
+       !$ACC PRESENT(e_p)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -2024,4 +2046,6 @@
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( intermediate_timestep_count == 1 )  THEN
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+             !$ACC PRESENT(tend, te_m)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -2033,4 +2057,6 @@
           ELSEIF ( intermediate_timestep_count < &
                    intermediate_timestep_count_max )  THEN
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+             !$ACC PRESENT(tend, te_m)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -2380,6 +2406,15 @@
 !-- points first, gradients at surface-bounded grid points will be
 !-- overwritten further below.
+    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) &
+    !$ACC PRIVATE(surf_s, surf_e) &
+    !$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) &
+    !$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_0) &
+    !$ACC PRESENT(tend) &
+    !$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) &
+    !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
+    !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3))
     DO  i = nxl, nxr
        DO  j = nys, nyn
+          !$ACC LOOP PRIVATE(k)
           DO  k = nzb+1, nzt
 
@@ -2423,4 +2458,5 @@
                 surf_s = surf_def_v(l)%start_index(j,i)
                 surf_e = surf_def_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k           = surf_def_v(l)%k(m)
@@ -2441,4 +2477,5 @@
                 surf_s = surf_lsm_v(l)%start_index(j,i)
                 surf_e = surf_lsm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k           = surf_lsm_v(l)%k(m)
@@ -2459,4 +2496,5 @@
                 surf_s = surf_usm_v(l)%start_index(j,i)
                 surf_e = surf_usm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k           = surf_usm_v(l)%k(m)
@@ -2479,4 +2517,5 @@
                 surf_s = surf_def_v(l)%start_index(j,i)
                 surf_e = surf_def_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k     = surf_def_v(l)%k(m)
@@ -2497,4 +2536,5 @@
                 surf_s = surf_lsm_v(l)%start_index(j,i)
                 surf_e = surf_lsm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k     = surf_lsm_v(l)%k(m)
@@ -2515,4 +2555,5 @@
                 surf_s = surf_usm_v(l)%start_index(j,i)
                 surf_e = surf_usm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k     = surf_usm_v(l)%k(m)
@@ -2534,4 +2575,5 @@
              surf_s = surf_def_h(0)%start_index(j,i)
              surf_e = surf_def_h(0)%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_def_h(0)%k(m)
@@ -2550,4 +2592,5 @@
              surf_s = surf_lsm_h%start_index(j,i)
              surf_e = surf_lsm_h%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_lsm_h%k(m)
@@ -2561,4 +2604,5 @@
              surf_s = surf_usm_h%start_index(j,i)
              surf_e = surf_usm_h%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_usm_h%k(m)
@@ -2573,4 +2617,5 @@
              surf_s = surf_def_h(1)%start_index(j,i)
              surf_e = surf_def_h(1)%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_def_h(1)%k(m)
@@ -2585,4 +2630,5 @@
 
 
+          !$ACC LOOP PRIVATE(k, def, flag)
           DO  k = nzb+1, nzt
 
@@ -2684,7 +2730,16 @@
           ELSE ! or IF ( .NOT. ocean_mode )  THEN
 
+             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+             !$ACC PRIVATE(surf_s, surf_e) &
+             !$ACC PRIVATE(tmp_flux(nzb+1:nzt)) &
+             !$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_0) &
+             !$ACC PRESENT(tend) &
+             !$ACC PRESENT(surf_def_h(0:2)) &
+             !$ACC PRESENT(surf_lsm_h) &
+             !$ACC PRESENT(surf_usm_h)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
 
+                   !$ACC LOOP PRIVATE(k)
                    DO  k = nzb+1, nzt
                       tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
@@ -2697,4 +2752,5 @@
                          surf_s = surf_def_h(l)%start_index(j,i)
                          surf_e = surf_def_h(l)%end_index(j,i)
+                         !$ACC LOOP PRIVATE(m, k)
                          DO  m = surf_s, surf_e
                             k = surf_def_h(l)%k(m)
@@ -2706,4 +2762,5 @@
                       surf_s = surf_lsm_h%start_index(j,i)
                       surf_e = surf_lsm_h%end_index(j,i)
+                      !$ACC LOOP PRIVATE(m, k)
                       DO  m = surf_s, surf_e
                          k = surf_lsm_h%k(m)
@@ -2714,4 +2771,5 @@
                       surf_s = surf_usm_h%start_index(j,i)
                       surf_e = surf_usm_h%end_index(j,i)
+                      !$ACC LOOP PRIVATE(m, k)
                       DO  m = surf_s, surf_e
                          k = surf_usm_h%k(m)
@@ -2723,4 +2781,5 @@
                       surf_s = surf_def_h(2)%start_index(j,i)
                       surf_e = surf_def_h(2)%end_index(j,i)
+                      !$ACC LOOP PRIVATE(m, k)
                       DO  m = surf_s, surf_e
                          k = surf_def_h(2)%k(m)
@@ -2732,4 +2791,5 @@
 
 !--                   Compute tendency for TKE-production from shear
+                     !$ACC LOOP PRIVATE(k, flag)
                       DO  k = nzb+1, nzt
                          flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
@@ -3628,5 +3688,9 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, drho_air, rho_air_zw
+        ONLY:  ddzu, dd2zu, ddzw, drho_air, rho_air_zw
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
+               wall_adjustment, wall_adjustment_factor
 
     USE grid_variables,                                                        &
@@ -3653,4 +3717,17 @@
     REAL(wp)     ::  ll             !< adjusted l
     REAL(wp)     ::  var_reference  !< reference temperature
+#ifdef _OPENACC
+!
+!-- From mixing_length_les:
+    REAL(wp)     :: l_stable        !< mixing length according to stratification
+!
+!-- From mixing_length_rans:
+    REAL(wp)     :: duv2_dz2        !< squared vertical gradient of wind vector
+    REAL(wp)     :: l_diss          !< mixing length for dissipation
+    REAL(wp)     :: rif             !< Richardson flux number
+!
+!-- From both:
+    REAL(wp)     :: dvar_dz         !< vertical gradient of var
+#endif
 
 #if defined( __nopointer )
@@ -3659,9 +3736,16 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
 #endif
-    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !< TKE dissipation
+    REAL(wp)     ::  dissipation  !< TKE dissipation
 
 
 !
 !-- Calculate the tendency terms
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+    !$ACC PRIVATE(flag, l, ll, dissipation) &
+    !$ACC PRIVATE(l_stable, duv2_dz2, l_diss, rif, dvar_dz) &
+    !$ACC PRESENT(e, u, v, km, var, wall_flags_0) &
+    !$ACC PRESENT(ddzu, dd2zu, ddzw, rho_air_zw, drho_air) &
+    !$ACC PRESENT(l_black, l_grid, l_wall) &
+    !$ACC PRESENT(diss, tend)
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -3675,20 +3759,80 @@
              IF ( les_dynamic .OR. les_mw )  THEN
 
+#ifdef _OPENACC
+!
+!--             Cannot call subroutine mixing_length_les because after adding all required
+!--             OpenACC directives the execution crashed reliably when inside the called
+!--             subroutine. I'm not sure why that is...
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+                   IF ( use_single_reference_value )  THEN
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                   ELSE
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+                ELSE
+                   l_stable = l_grid(k)
+                ENDIF
+!
+!--             Adjustment of the mixing length
+                IF ( wall_adjustment )  THEN
+                   l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
+                   ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
+                ELSE
+                   l  = MIN( l_grid(k), l_stable )
+                   ll = l_grid(k)
+                ENDIF
+#else
                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
-
-                dissipation(k,j) = ( 0.19_wp + 0.74_wp * l / ll )              &
-                                   * e(k,j,i) * SQRT( e(k,j,i) ) / l
+#endif
+
+                dissipation = ( 0.19_wp + 0.74_wp * l / ll )                   &
+                              * e(k,j,i) * SQRT( e(k,j,i) ) / l
 
              ELSEIF ( rans_tke_l )  THEN
 
+#ifdef _OPENACC
+!
+!--             Same reason as above...
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+
+                duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2                 &
+                           + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2                 &
+                           + 1E-30_wp
+
+                IF ( use_single_reference_value )  THEN
+                   rif = g / var_reference * dvar_dz / duv2_dz2
+                ELSE
+                   rif = g / var(k,j,i) * dvar_dz / duv2_dz2
+                ENDIF
+
+                rif = MAX( rif, -5.0_wp )
+                rif = MIN( rif,  1.0_wp )
+
+!
+!--             Calculate diabatic mixing length using Dyer-profile functions
+                IF ( rif >= 0.0_wp )  THEN
+                   l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
+                   l_diss = l
+                ELSE
+!
+!--                In case of unstable stratification, use mixing length of neutral case
+!--                for l, but consider profile functions for l_diss
+                   l      = l_wall(k,j,i)
+                   l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
+                ENDIF
+#else
                 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
-
-                dissipation(k,j) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
-
-                diss(k,j,i) = dissipation(k,j) * flag
+#endif
+
+                dissipation = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
+
+                diss(k,j,i) = dissipation * flag
 
              ELSEIF ( rans_tke_e )  THEN
 
-                dissipation(k,j) = diss(k,j,i)
+                dissipation = diss(k,j,i)
 
              ENDIF
@@ -3710,22 +3854,15 @@
                                            ) * ddzw(k) * drho_air(k)           &
                                          ) * flag * dsig_e                     &
-                          - dissipation(k,j) * flag
+                          - dissipation * flag
+!
+!--          Store dissipation if needed for calculating the sgs particle
+!--          velocities
+             IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.         &
+                  wang_kernel  .OR.  collision_turbulence  ) )  THEN
+                diss(k,j,i) = dissipation * flag
+             ENDIF
 
           ENDDO
        ENDDO
-
-!
-!--    Store dissipation if needed for calculating the sgs particle
-!--    velocities
-       IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.               &
-            wang_kernel  .OR.  collision_turbulence  ) )  THEN
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                diss(k,j,i) = dissipation(k,j) * MERGE( 1.0_wp, 0.0_wp,        &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-             ENDDO
-          ENDDO
-       ENDIF
-
     ENDDO
 
@@ -4257,4 +4394,6 @@
 !-- Downward facing surfaces
     !$OMP PARALLEL DO PRIVATE(i,j,k)
+    !$ACC PARALLEL LOOP PRIVATE(i,j,k) &
+    !$ACC PRESENT(bc_h(1), kh, km)
     DO  m = 1, bc_h(1)%ns
        i = bc_h(1)%i(m)
@@ -4267,4 +4406,6 @@
 !-- Downward facing surfaces
     !$OMP PARALLEL DO PRIVATE(i,j,k)
+    !$ACC PARALLEL LOOP PRIVATE(i,j,k) &
+    !$ACC PRESENT(bc_h(0), kh, km)
     DO  m = 1, bc_h(0)%ns
        i = bc_h(0)%i(m)
@@ -4277,4 +4418,6 @@
 !-- Model top
     !$OMP PARALLEL DO
+    !$ACC PARALLEL LOOP COLLAPSE(2) &
+    !$ACC PRESENT(kh, km)
     DO  i = nxlg, nxrg
        DO  j = nysg, nyng
@@ -4315,4 +4458,11 @@
  SUBROUTINE tcm_diffusivities_default( var, var_reference )
  
+    USE arrays_3d,                                                             &
+        ONLY:  dd2zu
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
+               wall_adjustment, wall_adjustment_factor
+
     USE statistics,                                                            &
         ONLY :  rmask, sums_l_l
@@ -4331,4 +4481,10 @@
     REAL(wp)     ::  ll                  !< adjusted mixing length
     REAL(wp)     ::  var_reference       !< reference temperature
+#ifdef _OPENACC
+!
+!-- From mixing_length_les:
+    REAL(wp)     :: l_stable        !< mixing length according to stratification
+    REAL(wp)     :: dvar_dz         !< vertical gradient of var
+#endif
 
 #if defined( __nopointer )
@@ -4344,5 +4500,7 @@
 !
 !-- Initialization for calculation of the mixing length profile
+    !$ACC KERNELS PRESENT(sums_l_l)
     sums_l_l = 0.0_wp
+    !$ACC END KERNELS
 
 !
@@ -4353,4 +4511,8 @@
     IF ( les_dynamic .OR. les_mw )  THEN
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(sr) &
+       !$ACC PRIVATE(flag, dvar_dz, l_stable, l, ll) &
+       !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) &
+       !$ACC PRESENT(kh, km, sums_l_l)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -4361,5 +4523,34 @@
 !
 !--             Determine the mixing length for LES closure
+#ifdef _OPENACC
+!
+!--             Cannot call subroutine mixing_length_les because after adding all required
+!--             OpenACC directives the execution crashed reliably when inside the called
+!--             subroutine. I'm not sure why that is...
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+                   IF ( use_single_reference_value )  THEN
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                   ELSE
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+                ELSE
+                   l_stable = l_grid(k)
+                ENDIF
+!
+!--             Adjustment of the mixing length
+                IF ( wall_adjustment )  THEN
+                   l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
+                   ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
+                ELSE
+                   l  = MIN( l_grid(k), l_stable )
+                   ll = l_grid(k)
+                ENDIF
+#else
                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
+#endif
+
 !
 !--             Compute diffusion coefficients for momentum and heat
@@ -4429,6 +4620,8 @@
     ENDIF
 
+    !$ACC KERNELS PRESENT(sums_l_l)
     sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
                                                 ! data output
+    !$ACC END KERNELS
 !$OMP END PARALLEL
 
