Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3775)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3776)
@@ -25,42 +25,45 @@
 ! -----------------
 ! $Id$
+! write out preprocessor directives; remove tailing whitespaces
+!
+! 3775 2019-03-04 12:40:20Z gronemeier
 ! removed unused variables
-! 
+!
 ! 3724 2019-02-06 16:28:23Z kanani
 ! Correct double-used log_point_s units
-! 
+!
 ! 3719 2019-02-06 13:10:18Z kanani
-! Changed log_point to log_point_s, otherwise this overlaps with 
+! Changed log_point to log_point_s, otherwise this overlaps with
 ! 'all progn.equations' cpu measurement.
-! 
+!
 ! 3684 2019-01-20 20:20:58Z knoop
 ! Remove unused variable simulated_time
-! 
+!
 ! 3636 2018-12-19 13:48:34Z raasch
 ! nopointer option removed
-! 
+!
 ! 3634 2018-12-18 12:31:28Z knoop
 ! 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
-! 
+!
 ! 3545 2018-11-21 11:19:41Z gronemeier
 ! - Set rans_mode according to value of turbulence_closure
 ! - removed debug output
-! 
+!
 ! 3430 2018-10-25 13:36:23Z maronga
 ! Added support for buildings in the dynamic SGS model
-! 
+!
 ! 3398 2018-10-22 19:30:24Z knoop
 ! Refactored production_e and production_e_ij (removed code redundancy)
-! 
+!
 ! 3386 2018-10-19 16:28:22Z gronemeier
 ! Renamed tcm_prognostic to tcm_prognostic_equations
-! 
+!
 ! 3385 2018-10-19 14:52:29Z knoop
 ! Restructured loops and ifs in production_e to ease vectorization on GPUs
-! 
+!
 ! 3300 2018-10-02 14:16:54Z gronemeier
 ! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall
@@ -68,25 +71,25 @@
 ! - removed maxloc call as this produced different results for different
 !   compiler options
-! 
+!
 ! 3294 2018-10-01 02:37:10Z raasch
 ! changes concerning modularization of ocean option
-! 
+!
 ! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
-! 
+!
 ! 3245 2018-09-13 14:08:16Z knoop
 ! unused variables removed, shortest_distance has wp now
-! 
+!
 ! 3183 2018-07-27 14:25:55Z suehring
 ! Rename variables and remove unused variable from USE statement
-! 
+!
 ! 3182 2018-07-27 13:36:03Z suehring
 ! Use MOST for km only in RANS mode
-! 
+!
 ! 3130 2018-07-16 11:08:55Z gronemeier
 ! - move boundary condition of km and kh to tcm_diffusivities
 ! - calculate km at boundaries according to MOST
 ! - move phi_m to surface_layer_fluxes_mod
-! 
+!
 ! 3129 2018-07-16 07:45:13Z gronemeier
 ! - move limitation of diss to boundary_conds
@@ -96,18 +99,18 @@
 ! - resort output after case select -> reduced code duplication
 ! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d
-! 
+!
 ! 3121 2018-07-11 18:46:49Z gronemeier
 ! - created the function phi_m
 ! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
 !   surfaces
-! 
+!
 ! 3120 2018-07-11 18:30:57Z gronemeier
 ! - changed tcm_diffusivities to tcm_diffusivities_default
 ! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
 !   and tcm_diffusivities_dynamic
-! 
+!
 ! 3086 2018-06-25 09:08:04Z gronemeier
 ! bugfix: set rans_const_sigma(1) = 1.3
-! 
+!
 ! 3083 2018-06-19 14:03:12Z gronemeier
 ! - set limits of diss at the end of prognostic equations
@@ -130,17 +133,17 @@
 !    - calculate km based on l_wall
 !    - initialize diss if 1D model is not used
-! 
+!
 ! 3045 2018-05-28 07:55:41Z Giersch
 ! Error message revised
-! 
+!
 ! 3014 2018-05-09 08:42:38Z maronga
 ! Bugfix: nzb_do and nzt_do were not used for 3d data output
-! 
+!
 ! 3004 2018-04-27 12:33:25Z Giersch
 ! Further allocation checks implemented
-! 
+!
 ! 2938 2018-03-27 15:52:42Z suehring
 ! Further todo's
-! 
+!
 ! 2936 2018-03-27 14:49:27Z gronemeier
 ! - defined l_grid only within this module
@@ -153,17 +156,17 @@
 ! - Moved init_mixing_length from init_grid.f90 and
 !   renamed it to tcm_init_mixing_length
-! 
+!
 ! 2764 2018-01-22 09:25:36Z gronemeier
 ! Bugfix: remove duplicate SAVE statements
-! 
+!
 ! 2746 2018-01-15 12:06:04Z suehring
 ! Move flag plant canopy to modules
-! 
+!
 ! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
-! 
+!
 ! 2701 2017-12-15 15:40:50Z suehring
 ! Changes from last commit documented
-! 
+!
 ! 2698 2017-12-14 18:46:24Z suehring
 ! Bugfix in get_topography_top_index
@@ -172,6 +175,6 @@
 ! Initial revision
 !
-! 
-! 
+!
+!
 !
 ! Authors:
@@ -191,5 +194,5 @@
 !------------------------------------------------------------------------------!
  MODULE turbulence_closure_mod
- 
+
 
     USE arrays_3d,                                                             &
@@ -264,5 +267,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid     !< geometric mean of grid sizes dx, dy, dz 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid     !< geometric mean of grid sizes dx, dy, dz
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
@@ -309,5 +312,5 @@
 
 !
-!-- Initialization actions  
+!-- Initialization actions
     INTERFACE tcm_init
        MODULE PROCEDURE tcm_init
@@ -476,5 +479,5 @@
        c_0 = 0.1_wp    !according to Lilly (1967) and Deardorff (1980)
 
-       dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 
+       dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
                        !of K_e which is used in RANS mode
 
@@ -489,5 +492,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_check_data_output( var, unit )
- 
+
     IMPLICIT NONE
 
@@ -519,5 +522,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
-    
+
     IMPLICIT NONE
 
@@ -567,5 +570,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_3d_data_averaging( mode, variable )
- 
+
 
     USE averaging,                                                             &
@@ -616,5 +619,5 @@
 
           CASE ( 'diss' )
-             IF ( ALLOCATED( diss_av ) ) THEN             
+             IF ( ALLOCATED( diss_av ) ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -662,5 +665,5 @@
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         diss_av(k,j,i) = diss_av(k,j,i)                       & 
+                         diss_av(k,j,i) = diss_av(k,j,i)                       &
                                         / REAL( average_count_3d, KIND=wp )
                       ENDDO
@@ -674,5 +677,5 @@
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         kh_av(k,j,i) = kh_av(k,j,i)                           & 
+                         kh_av(k,j,i) = kh_av(k,j,i)                           &
                                         / REAL( average_count_3d, KIND=wp )
                       ENDDO
@@ -686,5 +689,5 @@
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         km_av(k,j,i) = km_av(k,j,i)                           & 
+                         km_av(k,j,i) = km_av(k,j,i)                           &
                                         / REAL( average_count_3d, KIND=wp )
                       ENDDO
@@ -707,5 +710,5 @@
  SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
                                 nzb_do, nzt_do )
- 
+
     USE averaging,                                                             &
         ONLY:  diss_av, kh_av, km_av
@@ -734,5 +737,5 @@
 
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
-    
+
     found = .TRUE.
     resorted = .FALSE.
@@ -791,13 +794,13 @@
                 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
                                          REAL( fill_value, KIND = wp ),        &
-                                         BTEST( wall_flags_0(k,j,i), flag_nr ) ) 
+                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
        ENDDO
     ENDIF
- 
+
  END SUBROUTINE tcm_data_output_2d
 
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -806,5 +809,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
- 
+
 
     USE averaging,                                                             &
@@ -873,5 +876,5 @@
              to_be_resorted => km_av
           ENDIF
-          
+
        CASE DEFAULT
           found = .FALSE.
@@ -923,8 +926,8 @@
 
 !
-!-- Allocate arrays required for dissipation. 
+!-- Allocate arrays required for dissipation.
 !-- Please note, if it is a nested run, arrays need to be allocated even if
-!-- they do not necessarily need to be transferred, which is attributed to 
-!-- the design of the model coupler which allocates memory for each variable. 
+!-- they do not necessarily need to be transferred, which is attributed to
+!-- the design of the model coupler which allocates memory for each variable.
     IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.       &
          collision_turbulence  .OR.  nested_run )  THEN
@@ -1025,5 +1028,5 @@
              ELSE
                 IF ( .NOT. ocean_mode )  THEN
-                   kh   = 0.01_wp   ! there must exist an initial diffusion, because 
+                   kh   = 0.01_wp   ! there must exist an initial diffusion, because
                    km   = 0.01_wp   ! otherwise no TKE would be produced by the
                                     ! production terms, as long as not yet
@@ -1069,5 +1072,5 @@
           ELSE
              IF ( .NOT. ocean_mode )  THEN
-                kh   = 0.01_wp   ! there must exist an initial diffusion, because 
+                kh   = 0.01_wp   ! there must exist an initial diffusion, because
                 km   = 0.01_wp   ! otherwise no TKE would be produced by the
                                  ! production terms, as long as not yet
@@ -1112,5 +1115,5 @@
 !
 !--    In case of complex terrain and cyclic fill method as initialization,
-!--    shift initial data in the vertical direction for each point in the 
+!--    shift initial data in the vertical direction for each point in the
 !--    x-y-plane depending on local surface height
        IF ( complex_terrain  .AND.                                             &
@@ -1142,5 +1145,5 @@
           mean_inflow_profiles(:,5) = hom_sum(:,8,0)   ! e
 !
-!--       In case of complex terrain, determine vertical displacement at inflow 
+!--       In case of complex terrain, determine vertical displacement at inflow
 !--       boundary and adjust mean inflow profiles
           IF ( complex_terrain )  THEN
@@ -1344,5 +1347,5 @@
                                            SQRT( 0.25_wp * dy**2 +             &
                                           ( zw(k) - zu(k) )**2 ) )
-!   
+!
 !--                 xz-distance (vertical edges, down)
                     IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 )  .OR.       &
@@ -1384,5 +1387,5 @@
                                            SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
                                                  +  ( zw(k) - zu(k) )**2  ) )
-                     
+
                 ENDIF
              ENDDO
@@ -1426,5 +1429,5 @@
 !
 !--    Limit mixing length to either nearest wall or Blackadar mixing length.
-!--    For that, analyze each grid point (i/j/k) ("analysed grid point") and 
+!--    For that, analyze each grid point (i/j/k) ("analysed grid point") and
 !--    search within its vicinity for the shortest distance to a wall by cal-
 !--    culating the distance between the analysed grid point and the "viewed
@@ -1459,5 +1462,5 @@
 !
 !--       When analysed grid point lies above maximum topography, set search
-!--       radius to 0 if the distance between the analysed grid point and max 
+!--       radius to 0 if the distance between the analysed grid point and max
 !--       topography height is larger than the maximum search radius
           IF ( zu(k-rad_k_b) > zu(k_max_topo) )  rad_k_b = 0
@@ -1628,5 +1631,5 @@
 ! Description:
 ! ------------
-!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and 
+!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and
 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
 !> closest to the origin (0/0) of 'array'.
@@ -1694,5 +1697,5 @@
           ENDDO
        ENDIF
-       
+
     END FUNCTION
 
@@ -1700,5 +1703,5 @@
 ! Description:
 ! ------------
-!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point 
+!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point
 !> (kp,jp,ip) containing the first bit of wall_flags_0 into the array
 !> 'vicinity'. Only copy first bit as this indicates the presence of topography.
@@ -1753,5 +1756,5 @@
 
     USE surface_mod,                                                           &
-        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h 
+        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -1761,5 +1764,5 @@
     INTEGER(iwp) ::  k      !< grid index z-direction
     INTEGER(iwp) ::  m      !< running index surface elements
-    
+
     REAL(wp) ::  km_sfc     !< diffusion coefficient, used to compute virtual velocities
 
@@ -1780,10 +1783,10 @@
        DO  m = 1, surf_def_h(0)%ns
 
-          i = surf_def_h(0)%i(m)            
+          i = surf_def_h(0)%i(m)
           j = surf_def_h(0)%j(m)
           k = surf_def_h(0)%k(m)
 !
-!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
-!--       and km are not on the same grid. Actually, a further 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
+!--       and km are not on the same grid. Actually, a further
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
@@ -1798,5 +1801,5 @@
                                      drho_air_zw(k-1)               *          &
                                      ( zu(k+1) - zu(k-1)    )       /          &
-                                     ( km_sfc  + 1.0E-20_wp )   
+                                     ( km_sfc  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
@@ -1816,10 +1819,10 @@
        DO  m = 1, surf_def_h(1)%ns
 
-          i = surf_def_h(1)%i(m)            
+          i = surf_def_h(1)%i(m)
           j = surf_def_h(1)%j(m)
           k = surf_def_h(1)%k(m)
 !
-!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
-!--       and km are not on the same grid. Actually, a further 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
+!--       and km are not on the same grid. Actually, a further
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
@@ -1827,9 +1830,9 @@
                                      drho_air_zw(k-1) *                        &
                                      ( zu(k+1)    - zu(k-1)    )  /            &
-                                     ( km(k,j,i)  + 1.0E-20_wp )  
+                                     ( km(k,j,i)  + 1.0E-20_wp )
           surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *          &
                                      drho_air_zw(k-1) *                        &
                                      ( zu(k+1)    - zu(k-1)    )  /            &
-                                     ( km(k,j,i)  + 1.0E-20_wp )  
+                                     ( km(k,j,i)  + 1.0E-20_wp )
 
           IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                     &
@@ -1853,6 +1856,6 @@
           k = surf_lsm_h%k(m)
 !
-!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
-!--       and km are not on the same grid. Actually, a further 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
+!--       and km are not on the same grid. Actually, a further
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
@@ -1863,5 +1866,5 @@
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km_sfc  + 1.0E-20_wp )  
+                                        ( km_sfc  + 1.0E-20_wp )
           surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
                                         drho_air_zw(k-1)         *             &
@@ -1889,6 +1892,6 @@
           k = surf_usm_h%k(m)
 !
-!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v 
-!--       and km are not on the same grid. Actually, a further 
+!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
+!--       and km are not on the same grid. Actually, a further
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
@@ -2365,5 +2368,5 @@
     REAL(wp)     ::  theta       !< virtual potential temperature
     REAL(wp)     ::  temp        !< theta * Exner-function
-    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation 
+    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
     REAL(wp)     ::  usvs        !< momentum flux u"v"
     REAL(wp)     ::  vsus        !< momentum flux v"u"
@@ -3669,5 +3672,5 @@
  SUBROUTINE diffusion_e( var, var_reference )
 
-#ifdef _OPENACC
+#if defined( _OPENACC )
     USE arrays_3d,                                                             &
         ONLY:  ddzu, dd2zu, ddzw, drho_air, rho_air_zw
@@ -3703,5 +3706,5 @@
     REAL(wp)     ::  ll             !< adjusted l
     REAL(wp)     ::  var_reference  !< reference temperature
-#ifdef _OPENACC
+#if defined( _OPENACC )
 !
 !-- From mixing_length_les:
@@ -3742,5 +3745,5 @@
              IF ( les_dynamic .OR. les_mw )  THEN
 
-#ifdef _OPENACC
+#if defined( _OPENACC )
 !
 !--             Cannot call subroutine mixing_length_les because after adding all required
@@ -3777,5 +3780,5 @@
              ELSEIF ( rans_tke_l )  THEN
 
-#ifdef _OPENACC
+#if defined( _OPENACC )
 !
 !--             Same reason as above...
@@ -3857,5 +3860,5 @@
 !--    Upward facing surfaces
        DO  m = 1, bc_h(0)%ns
-          i = bc_h(0)%i(m)            
+          i = bc_h(0)%i(m)
           j = bc_h(0)%j(m)
           k = bc_h(0)%k(m)
@@ -3865,5 +3868,5 @@
 !--    Downward facing surfaces
        DO  m = 1, bc_h(1)%ns
-          i = bc_h(1)%i(m)            
+          i = bc_h(1)%i(m)
           j = bc_h(1)%j(m)
           k = bc_h(1)%k(m)
@@ -3889,5 +3892,5 @@
     USE grid_variables,                                                        &
         ONLY:  ddx2, ddy2
-        
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  collision_turbulence
@@ -3984,6 +3987,6 @@
 !--    For each surface type determine start and end index (in case of elevated
 !--    topography several up/downward facing surfaces may exist.
-       surf_s = bc_h(0)%start_index(j,i)   
-       surf_e = bc_h(0)%end_index(j,i)   
+       surf_s = bc_h(0)%start_index(j,i)
+       surf_e = bc_h(0)%end_index(j,i)
        DO  m = surf_s, surf_e
           k             = bc_h(0)%k(m)
@@ -3992,6 +3995,6 @@
 !
 !--    Downward facing surfaces
-       surf_s = bc_h(1)%start_index(j,i)   
-       surf_e = bc_h(1)%end_index(j,i)   
+       surf_s = bc_h(1)%start_index(j,i)
+       surf_e = bc_h(1)%end_index(j,i)
        DO  m = surf_s, surf_e
           k             = bc_h(1)%k(m)
@@ -4298,5 +4301,5 @@
           kh(k,j,i) = 1.35_wp * km(k,j,i)
        ENDDO
-!      
+!
 !--    Natural surfaces
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
@@ -4309,5 +4312,5 @@
           kh(k,j,i) = 1.35_wp * km(k,j,i)
        ENDDO
-!      
+!
 !--    Urban surfaces
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
@@ -4320,10 +4323,10 @@
           kh(k,j,i) = 1.35_wp * km(k,j,i)
        ENDDO
-       
-!      
+
+!
 !--    North-, south-, west and eastward facing surfaces
 !--    Do not consider stratification at these surfaces.
        DO  n = 0, 3
-!      
+!
 !--       Default surfaces
           !$OMP PARALLEL DO PRIVATE(i,j,k,m)
@@ -4335,5 +4338,5 @@
              kh(k,j,i) = 1.35_wp * km(k,j,i)
           ENDDO
-!      
+!
 !--       Natural surfaces
           !$OMP PARALLEL DO PRIVATE(i,j,k,m)
@@ -4345,5 +4348,5 @@
              kh(k,j,i) = 1.35_wp * km(k,j,i)
           ENDDO
-!      
+!
 !--       Urban surfaces
           !$OMP PARALLEL DO PRIVATE(i,j,k,m)
@@ -4356,5 +4359,5 @@
           ENDDO
        ENDDO
-       
+
        CALL exchange_horiz( km, nbgp )
        CALL exchange_horiz( kh, nbgp )
@@ -4424,10 +4427,10 @@
 ! Description:
 ! ------------
-!> Computation of the turbulent diffusion coefficients for momentum and heat 
+!> Computation of the turbulent diffusion coefficients for momentum and heat
 !> according to Prandtl-Kolmogorov.
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_diffusivities_default( var, var_reference )
- 
-#ifdef _OPENACC
+
+#if defined( _OPENACC )
     USE arrays_3d,                                                             &
         ONLY:  dd2zu
@@ -4454,5 +4457,5 @@
     REAL(wp)     ::  ll                  !< adjusted mixing length
     REAL(wp)     ::  var_reference       !< reference temperature
-#ifdef _OPENACC
+#if defined( _OPENACC )
 !
 !-- From mixing_length_les:
@@ -4493,5 +4496,5 @@
 !
 !--             Determine the mixing length for LES closure
-#ifdef _OPENACC
+#if defined( _OPENACC )
 !
 !--             Cannot call subroutine mixing_length_les because after adding all required
@@ -4602,13 +4605,13 @@
 ! Description:
 ! ------------
-!> Calculates the eddy viscosity dynamically using the linear dynamic model 
-!> according to 
-!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via 
-!> stochastic analysis." 
+!> Calculates the eddy viscosity dynamically using the linear dynamic model
+!> according to
+!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via
+!> stochastic analysis."
 !> Monte Carlo Methods and Applications 14.4 (2008): 311-329.
 !>
-!> Furthermore dynamic bounds are used to limit the absolute value of c* as 
+!> Furthermore dynamic bounds are used to limit the absolute value of c* as
 !> described in
-!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation: 
+!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation:
 !> Stability via realizability." Physics of Fluids 29.10 (2017): 105104.
 !>
@@ -4617,8 +4620,8 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_diffusivities_dynamic
-        
+
     USE arrays_3d,                                                             &
         ONLY:  ddzw, dzw, dd2zu, w, ug, vg
-    
+
     USE grid_variables,                                                        &
         ONLY : ddx, ddy, dx, dy
@@ -4643,5 +4646,5 @@
 
     REAL(wp)     ::  flag        !< topography flag
-    
+
     REAL(wp)     ::  uc(-1:1,-1:1)  !< u on grid center
     REAL(wp)     ::  vc(-1:1,-1:1)  !< v on grid center
@@ -4706,7 +4709,7 @@
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-            
+
              flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-             
+
 !
 !--          Compute the deviatoric shear tensor s_ij on grid centers:
@@ -4829,5 +4832,5 @@
              cst_max = fac_cmax * SQRT( e(k,j,i) ) /                           &
                        ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
-             
+
              IF ( ABS( cst ) > cst_max )  THEN
                 cst = cst_max * cst / ABS( cst )
@@ -4835,5 +4838,5 @@
 
              km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag
-             
+
           ENDDO
        ENDDO
@@ -4850,5 +4853,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_box_filter_2d_single( var, var_fil )
- 
+
     IMPLICIT NONE
 
@@ -4856,5 +4859,5 @@
     REAL(wp)     ::  var_fil             !< filtered variable
 !
-!-- It is assumed that a box with a side length of 2 * dx and centered at the 
+!-- It is assumed that a box with a side length of 2 * dx and centered at the
 !-- variable*s position contains one half of the four closest neigbours and one
 !-- forth of the diagonally closest neighbours.
@@ -4884,5 +4887,5 @@
     REAL(wp)     ::  var_fil(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< filtered variable
 !
-!-- It is assumed that a box with a side length of 2 * dx and centered at the 
+!-- It is assumed that a box with a side length of 2 * dx and centered at the
 !-- variable's position contains one half of the four closest neigbours and one
 !-- forth of the diagonally closest neighbours.
