Changeset 2553

Timestamp:

Oct 18, 2017 8:03:45 AM (7 years ago)

Author:

Giersch

Message:

Changed programming style and bugfixes

File:

• palm/trunk/SOURCE/wind_turbine_model_mod.f90 (modified) (18 diffs)

palm/trunk/SOURCE/wind_turbine_model_mod.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! Bugfix of vertical loop in wtm_tendencies to account for different turbine 
+! heights, bugfix of the interpolation of the u-component concerning the 
+! vertical index and further small adjustments of the programming style 
+! 
+! 2410 2017-09-06 08:16:30Z Giersch
 ! Revise error message PA0462
 ! 
@@ -93 +98 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dd2zu, tend, u, v, w, zu, zw
-
-    USE constants
+        ONLY:  tend, u, v, w, zu, zw
+
+    USE constants,                                                             &
+        ONLY:  pi
 
     USE control_parameters,                                                    &
@@ -140 +146 @@
 
     REAL(wp), DIMENSION(1:100) ::  dtow             = 0.0_wp  !< tower diameter [m]
-    REAL(wp), DIMENSION(1:100) ::  omega_rot        = 0.0_wp  !< inital or constant rotor speed [rad/s]
+    REAL(wp), DIMENSION(1:100) ::  omega_rot        = 0.9_wp  !< inital or constant rotor speed [rad/s]
     REAL(wp), DIMENSION(1:100) ::  phi_yaw          = 0.0_wp  !< yaw angle [degree] ( clockwise, 0 = facing west )
     REAL(wp), DIMENSION(1:100) ::  pitch_add        = 0.0_wp  !< constant pitch angle
@@ -240 +246 @@
      
     REAL(wp) ::  cur_r                       !< 
-    REAL(wp) ::  delta_t                     !<  tangential segment length
     REAL(wp) ::  phi_rotor                   !< 
     REAL(wp) ::  pre_factor                  !<  
@@ -383 +388 @@
     INTERFACE wtm_forces
        MODULE PROCEDURE wtm_forces
+    END INTERFACE wtm_forces
+
+    INTERFACE wtm_yawcontrol
        MODULE PROCEDURE wtm_yawcontrol
-    END INTERFACE wtm_forces
+    END INTERFACE wtm_yawcontrol
     
     INTERFACE wtm_rotate_rotor
@@ -402 +410 @@
     
     PUBLIC wtm_check_parameters, wtm_forces, wtm_init, wtm_init_arrays,        &
-           wtm_parin, wtm_tendencies, wtm_tendencies_ij, wind_turbine 
+           wtm_parin, wtm_tendencies, wind_turbine 
 
  CONTAINS
@@ -524 +532 @@
        delta_t_factor = segment_length
        delta_r_init   = delta_r_factor * MIN( dx, dy, dz)
-       delta_t        = delta_t_factor * MIN( dx, dy, dz)
 
        DO inot = 1, nturbines
@@ -1484 +1491 @@
                    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 ) / dz
 !
 !--                Interpolate only if all required information is available on
@@ -1710 +1717 @@
 !
 !--                In the following the 3D-velocity field is projected its
-!--                components perpedicular and parallel to the rotor area
+!--                components perpendicular and parallel to the rotor area
 !--                The calculation of forces will be done in the rotor-
 !--                coordinates y' and z.
@@ -2230 +2237 @@
     SUBROUTINE wtm_yawcontrol( inot )
     
-       USE constants
        USE kinds
                 
@@ -2461 +2467 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   DO  k = nzb+1, k_hub(1) + k_smear(1)
+                   DO  k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear)
 !
 !--                   Calculate the thrust generated by the nacelle and the tower
                       tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) *               &
-                                         SIGN( u(k,j,i)**2 , u(k,j,i) )     
+                                      SIGN( u(k,j,i)**2 , u(k,j,i) )     
                       tend_tow_x   = 0.5_wp * tow_cd_surf(k,j,i) *             &
                                          SIGN( u(k,j,i)**2 , u(k,j,i) ) 
-                                                   
+                                                
                       tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i)        &
                                   - tend_nac_x - tend_tow_x )                  &
@@ -2484 +2490 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   DO  k = nzb+1, k_hub(1) + k_smear(1)
+                   DO  k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear)
                       tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) *               &
                                          SIGN( v(k,j,i)**2 , v(k,j,i) )     
@@ -2504 +2510 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   DO  k = nzb+1,  k_hub(1) + k_smear(1)
+                   DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
                       tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i)            &
                                       * MERGE( 1.0_wp, 0.0_wp,                 &
@@ -2548 +2554 @@
 !--       Apply the x-component of the force to the u-component of the flow:
           IF ( simulated_time >= time_turbine_on )  THEN
-
-             DO  k = nzb+1,  k_hub(1) + k_smear(1)
+             DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
 !
 !--             Calculate the thrust generated by the nacelle and the tower 
@@ -2559 +2564 @@
                             - tend_nac_x - tend_tow_x )                        &
                                       * MERGE( 1.0_wp, 0.0_wp,                 &
-                                               BTEST( wall_flags_0(k,j,i), 1 ) )
+                                            BTEST( wall_flags_0(k,j,i), 1 ) )
              ENDDO
           ENDIF
@@ -2567 +2572 @@
 !--       Apply the y-component of the force to the v-component of the flow:
           IF ( simulated_time >= time_turbine_on )  THEN
-             DO  k = nzb+1,  k_hub(1) + k_smear(1)
+             DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
                 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) *                     &
                                    SIGN( v(k,j,i)**2 , v(k,j,i) )     
@@ -2576 +2581 @@
                                       * MERGE( 1.0_wp, 0.0_wp,                 &
                                                BTEST( wall_flags_0(k,j,i), 2 ) )
-             ENDDO
+                ENDDO
           ENDIF
 
@@ -2583 +2588 @@
 !--       Apply the z-component of the force to the w-component of the flow:
           IF ( simulated_time >= time_turbine_on )  THEN
-             DO  k = nzb+1,  k_hub(1) + k_smear(1)
+             DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
                 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i)                  &
                                       * MERGE( 1.0_wp, 0.0_wp,                 &