Index: palm/trunk/SOURCE/wind_turbine_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 3049)
+++ palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 3065)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! dz was replaced by dz(1), error message concerning grid stretching has been
+! introduced
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -132,4 +136,5 @@
 !> automatically).
 !>
+!> @todo Replace dz(1) appropriatly to account for grid stretching
 !> @todo Revise code according to PALM Coding Standard
 !> @todo Implement ADM and ALM turbine models
@@ -740,5 +745,5 @@
        delta_r_factor = segment_width
        delta_t_factor = segment_length
-       delta_r_init   = delta_r_factor * MIN( dx, dy, dz)
+       delta_r_init   = delta_r_factor * MIN( dx, dy, dz(1))
 
        DO inot = 1, nturbines
@@ -937,4 +942,7 @@
 
     
+       USE control_parameters,                                                 &
+           ONLY:  dz_stretch_level_start
+    
        IMPLICIT NONE
 
@@ -992,4 +1000,17 @@
                  pi**( 1.0_wp / 6.0_wp ) * eps_kernel
 !
+!--    Stretching (non-uniform grid spacing) is not considered in the wind 
+!--    turbine model. Therefore, vertical stretching has to be applied above
+!--    the area where the wtm is active. ABS (...) is required because the 
+!--    default value of dz_stretch_level_start is -9999999.9_wp (negative). 
+       IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL(rcz) + MAXVAL(rr) +     &
+       eps_min)  THEN
+          WRITE( message_string, * ) 'The lowest level where vertical ',       &
+                                     'stretching is applied &have to be ',     &
+                                     'greater than ', MAXVAL(rcz) +            &
+                                      MAXVAL(rr) + eps_min
+          CALL message( 'init_grid', 'PA0484', 1, 2, 0, 6, 0 )
+       ENDIF 
+!
 !--    Square of eps_min:
        eps_min2 = eps_min**2
@@ -1022,5 +1043,5 @@
           i_hub(inot) = INT(   rcx(inot)                 / dx )
           j_hub(inot) = INT( ( rcy(inot) + 0.5_wp * dy ) / dy )
-          k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz ) / dz )
+          k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz(1) ) / dz(1) )
 
 !
@@ -1031,5 +1052,5 @@
           i_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dx )
           j_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dy )
-          k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz )
+          k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz(1) )
        
        ENDDO
@@ -1121,6 +1142,6 @@
 !--       surface. All points between z=0 and z=dz/s would already be contained
 !--       in grid box 1.
-          index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz ) + 1
-          index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz ) + 1
+          index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz(1) ) + 1
+          index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz(1) ) + 1
 
 !
@@ -1148,5 +1169,5 @@
                             IF ( j == tower_s )  THEN
                                tow_cd_surf(k,j,i) = ( rcz(inot) -              &
-                                    ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    & ! extension in z-direction
+                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction
                                   ( ( tower_s + 1.0_wp + 0.5_wp ) * dy    -    &
                                     ( rcy(inot) - 0.5_wp * dtow(inot) ) ) *    & ! extension in y-direction
@@ -1154,5 +1175,5 @@
                             ELSEIF ( j == tower_n )  THEN
                                tow_cd_surf(k,j,i) = ( rcz(inot)            -   &
-                                    ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    & ! extension in z-direction
+                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction
                                   ( ( rcy(inot) + 0.5_wp * dtow(inot) )   -    &
                                     ( tower_n + 0.5_wp ) * dy )           *    & ! extension in y-direction
@@ -1163,5 +1184,5 @@
                             ELSE
                                tow_cd_surf(k,j,i) = ( rcz(inot) -              &
-                                    ( k_hub(inot) * dz - 0.5_wp * dz ) )  *    &
+                                    ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*&
                                     dy * turb_cd_tower(inot)
                             ENDIF
@@ -1169,6 +1190,6 @@
 !--                      tower lies completely within one grid box:
                          ELSE
-                            tow_cd_surf(k,j,i) = ( rcz(inot)                 - &
-                                       ( k_hub(inot) * dz - 0.5_wp * dz ) ) *  &
+                            tow_cd_surf(k,j,i) = ( rcz(inot) - ( k_hub(inot) * &
+                                       dz(1) - 0.5_wp * dz(1) ) ) *            &
                                        dtow(inot) * turb_cd_tower(inot)
                          ENDIF
@@ -1182,10 +1203,10 @@
 !--                         leftmost and rightmost grid box:
                             IF ( j == tower_s )  THEN                         
-                               tow_cd_surf(k,j,i) = dz * (                     &
+                               tow_cd_surf(k,j,i) = dz(1) * (                  &
                                       ( tower_s + 1 + 0.5_wp ) * dy         -  &
                                       ( rcy(inot) - 0.5_wp * dtow(inot) )      &
                                                         ) * turb_cd_tower(inot)
                             ELSEIF ( j == tower_n )  THEN
-                               tow_cd_surf(k,j,i) = dz * (                     &
+                               tow_cd_surf(k,j,i) = dz(1) * (                  &
                                       ( rcy(inot) + 0.5_wp * dtow(inot) )   -  &
                                       ( tower_n + 0.5_wp ) * dy                &
@@ -1195,10 +1216,11 @@
 !--                         (where tow_cd_surf = grid box area):
                             ELSE
-                               tow_cd_surf(k,j,i) = dz * dy * turb_cd_tower(inot)
+                               tow_cd_surf(k,j,i) = dz(1) * dy *               &
+                                                    turb_cd_tower(inot)
                             ENDIF
 !
 !--                         tower lies completely within one grid box:
                          ELSE
-                            tow_cd_surf(k,j,i) = dz * dtow(inot) *             &
+                            tow_cd_surf(k,j,i) = dz(1) * dtow(inot) *          &
                                                 turb_cd_tower(inot)
                          ENDIF ! end if larger than grid box
@@ -1302,6 +1324,6 @@
        ENDDO   ! end of loop over turbines
 
-       tow_cd_surf   = tow_cd_surf   / ( dx * dy * dz )      ! Normalize tower drag
-       nac_cd_surf = nac_cd_surf / ( dx * dy * dz )      ! Normalize nacelle drag
+       tow_cd_surf   = tow_cd_surf   / ( dx * dy * dz(1) )  ! Normalize tower drag
+       nac_cd_surf = nac_cd_surf / ( dx * dy * dz(1) )      ! Normalize nacelle drag
 
        CALL wtm_read_blade_tables
@@ -1746,5 +1768,5 @@
                    ii =   rbx(ring,rseg) * ddx
                    jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
-                   kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz
+                   kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1)
 !
 !--                Interpolate only if all required information is available on
@@ -1776,5 +1798,5 @@
 
                          u_int_1_l(inot,ring,rseg) = u_int_l          +        &
-                                     ( rbz(ring,rseg) - zu(kk) ) / dz *        &
+                                     ( rbz(ring,rseg) - zu(kk) ) / dz(1) *     &
                                      ( u_int_u - u_int_l )
 
@@ -1791,5 +1813,5 @@
                    ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx
                    jj =   rby(ring,rseg)                 * ddy
-                   kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz 
+                   kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 
 !
 !--                Interpolate only if all required information is available on
@@ -1821,5 +1843,5 @@
 
                          v_int_1_l(inot,ring,rseg) = v_int_l +                 &
-                                     ( rbz(ring,rseg) - zu(kk) ) / dz *        &
+                                     ( rbz(ring,rseg) - zu(kk) ) / dz(1) *     &
                                      ( v_int_u - v_int_l )
 
@@ -1836,5 +1858,5 @@
                    ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx
                    jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy
-                   kk =   rbz(ring,rseg)                 / dz
+                   kk =   rbz(ring,rseg)                 / dz(1)
 !
 !--                Interpolate only if all required information is available on
@@ -1866,5 +1888,5 @@
 
                          w_int_1_l(inot,ring,rseg) = w_int_l +                 &
-                                     ( rbz(ring,rseg) - zw(kk) ) / dz *        &
+                                     ( rbz(ring,rseg) - zw(kk) ) / dz(1) *     &
                                      ( w_int_u - w_int_l )
                       ELSE
@@ -2273,11 +2295,11 @@
                             dist_u_3d = ( i * dx               - rbx(ring,rseg) )**2 + &
                                         ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + &
-                                        ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2
+                                        ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2
                             dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + &
                                         ( j * dy               - rby(ring,rseg) )**2 + &
-                                        ( k * dz - 0.5_wp * dz - rbz(ring,rseg) )**2
+                                        ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2
                             dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + &
                                         ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + &
-                                        ( k * dz               - rbz(ring,rseg) )**2
+                                        ( k * dz(1)               - rbz(ring,rseg) )**2
 
 !
