Changeset 4497 for palm/trunk/SOURCE/wind_turbine_model_mod.f90
- Timestamp:
- Apr 15, 2020 10:20:51 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r4495 r4497 1 1 !> @file wind_turbine_model_mod.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the 6 ! terms of the GNU General Public License as published by the Free Software 7 ! Foundation, either version 3 of the License, or (at your option) any later 8 ! version. 9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANYr 11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 12 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 13 ! 14 ! You should have received a copy of the GNU General Public License along with 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 2009-2020 Carl von Ossietzky Universitaet Oldenburg 18 17 ! Copyright 1997-2020 Leibniz Universitaet Hannover 19 !------------------------------------------------------------------------------! 18 !--------------------------------------------------------------------------------------------------! 19 ! 20 20 ! 21 21 ! Current revisions: … … 26 26 ! ----------------- 27 27 ! $Id$ 28 ! file re-formatted to follow the PALM coding standard 29 ! 30 ! 31 ! 4495 2020-04-13 20:11:20Z raasch 28 32 ! restart data handling with MPI-IO added 29 33 ! 30 34 ! 4481 2020-03-31 18:55:54Z maronga 31 35 ! ASCII output cleanup 32 ! 36 ! 33 37 ! 4465 2020-03-20 11:35:48Z maronga 34 ! Removed old ASCII outputm, some syntax layout adjustments, added output for 35 ! rotor and tower diameters. Added warning message in case of NetCDF 3 (no 36 ! WTM output file will be produced). 37 ! 38 ! Removed old ASCII outputm, some syntax layout adjustments, added output for rotor and tower 39 ! diameters. Added warning message in case of NetCDF 3 (no WTM output file will be produced). 40 ! 38 41 ! 4460 2020-03-12 16:47:30Z oliver.maas 39 42 ! allow for simulating up to 10 000 wind turbines 40 ! 43 ! 41 44 ! 4459 2020-03-12 09:35:23Z oliver.maas 42 45 ! avoid division by zero in tip loss correction factor calculation 43 ! 46 ! 44 47 ! 4457 2020-03-11 14:20:43Z raasch 45 48 ! use statement for exchange horiz added 46 ! 49 ! 47 50 ! 4447 2020-03-06 11:05:30Z oliver.maas 48 51 ! renamed wind_turbine_parameters namelist variables 49 ! 52 ! 50 53 ! 4438 2020-03-03 20:49:28Z suehring 51 54 ! Bugfix: shifted netcdf preprocessor directive to correct position 52 ! 55 ! 53 56 ! 4436 2020-03-03 12:47:02Z oliver.maas 54 57 ! added optional netcdf data input for wtm array input parameters 55 ! 58 ! 56 59 ! 4426 2020-02-27 10:02:19Z oliver.maas 57 ! define time as unlimited dimension so that no maximum number 58 ! of time steps has to be given forwtm_data_output59 ! 60 ! define time as unlimited dimension so that no maximum number of time steps has to be given for 61 ! wtm_data_output 62 ! 60 63 ! 4423 2020-02-25 07:17:11Z maronga 61 64 ! Switched to serial output as data is aggerated before anyway. 62 ! 65 ! 63 66 ! 4420 2020-02-24 14:13:56Z maronga 64 67 ! Added output control for wind turbine model 65 ! 68 ! 66 69 ! 4412 2020-02-19 14:53:13Z maronga 67 70 ! Bugfix: corrected character length in dimension_names 68 ! 71 ! 69 72 ! 4411 2020-02-18 14:28:02Z maronga 70 73 ! Added output in NetCDF format using DOM (only netcdf4-parallel is supported). 71 74 ! Old ASCII output is still available at the moment. 72 ! 75 ! 73 76 ! 4360 2020-01-07 11:25:50Z suehring 74 ! Introduction of wall_flags_total_0, which currently sets bits based on static 75 ! topographyinformation used in wall_flags_static_076 ! 77 ! Introduction of wall_flags_total_0, which currently sets bits based on static topography 78 ! information used in wall_flags_static_0 79 ! 77 80 ! 4343 2019-12-17 12:26:12Z oliver.maas 78 ! replaced <= by < in line 1464 to ensure that ialpha will not be 79 ! greater than dlen 80 ! 81 ! replaced <= by < in line 1464 to ensure that ialpha will not be greater than dlen 82 ! 81 83 ! 4329 2019-12-10 15:46:36Z motisi 82 84 ! Renamed wall_flags_0 to wall_flags_static_0 83 ! 85 ! 84 86 ! 4326 2019-12-06 14:16:14Z oliver.maas 85 ! changed format of turbine control output to allow for higher 86 ! torque and power values 87 ! 87 ! changed format of turbine control output to allow for higher torque and power values 88 ! 88 89 ! 4182 2019-08-22 15:20:23Z scharf 89 90 ! Corrected "Former revisions" section 90 ! 91 ! 91 92 ! 4144 2019-08-06 09:11:47Z raasch 92 93 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc. 93 ! 94 ! 94 95 ! 4056 2019-06-27 13:53:16Z Giersch 95 ! CASE DEFAULT action in wtm_actions needs to be CONTINUE. Otherwise an abort 96 ! will happen for location values that are not implemented as CASE statements97 ! but are already realized in the code (e.g. pt-tendency)98 ! 96 ! CASE DEFAULT action in wtm_actions needs to be CONTINUE. Otherwise an abort will happen for 97 ! location values that are not implemented as CASE statements but are already realized in the code 98 ! (e.g. pt-tendency) 99 ! 99 100 ! 3885 2019-04-11 11:29:34Z kanani 100 ! Changes related to global restructuring of location messages and introduction 101 ! of additional debugmessages102 ! 101 ! Changes related to global restructuring of location messages and introduction of additional debug 102 ! messages 103 ! 103 104 ! 3875 2019-04-08 17:35:12Z knoop 104 105 ! Addaped wtm_tendency to fit the module actions interface 105 ! 106 ! 106 107 ! 3832 2019-03-28 13:16:58Z raasch 107 108 ! instrumented with openmp directives 108 ! 109 ! 109 110 ! 3725 2019-02-07 10:11:02Z raasch 110 111 ! unused variables removed 111 ! 112 ! 112 113 ! 3685 2019-01-21 01:02:11Z knoop 113 114 ! Some interface calls moved to module_interface + cleanup 114 ! 115 ! 115 116 ! 3655 2019-01-07 16:51:22Z knoop 116 117 ! Replace degree symbol by 'degrees' 117 ! 118 ! 118 119 ! 1914 2016-05-26 14:44:07Z witha 119 120 ! Initial revision 120 ! 121 ! 121 122 ! 122 123 ! Description: 123 124 ! ------------ 124 !> This module calculates the effect of wind turbines on the flow fields. The 125 !> initial version contains only the advanced actuator disk with rotation method 126 !> (ADM-R). 125 !> This module calculates the effect of wind turbines on the flow fields. The initial version 126 !> contains only the advanced actuator disk with rotation method (ADM-R). 127 127 !> The wind turbines include the tower effect, can be yawed and tilted. 128 !> The wind turbine model includes controllers for rotational speed, pitch and 129 !> yaw. 130 !> Currently some specifications of the NREL 5 MW reference turbine 131 !> are hardcoded whereas most input data comes from separate files (currently 132 !> external, planned to be included as namelist which will be read in 133 !> automatically). 128 !> The wind turbine model includes controllers for rotational speed, pitch and yaw. 129 !> Currently some specifications of the NREL 5 MW reference turbine are hardcoded whereas most input 130 !> data comes from separate files (currently external, planned to be included as namelist which will 131 !> be read in automatically). 134 132 !> 135 133 !> @todo Replace dz(1) appropriatly to account for grid stretching … … 145 143 !> @todo Allow different turbine types in one simulation 146 144 ! 147 !------------------------------------------------------------------------------ !145 !--------------------------------------------------------------------------------------------------! 148 146 MODULE wind_turbine_model_mod 149 147 150 USE arrays_3d, &148 USE arrays_3d, & 151 149 ONLY: tend, u, v, w, zu, zw 152 150 153 USE basic_constants_and_equations_mod, &151 USE basic_constants_and_equations_mod, & 154 152 ONLY: pi 155 153 156 USE control_parameters, & 157 ONLY: coupling_char, & 158 debug_output, & 159 dt_3d, dz, end_time, message_string, restart_data_format_output, time_since_reference_point,& 160 wind_turbine, initializing_actions, origin_date_time 161 162 USE cpulog, & 154 USE control_parameters, & 155 ONLY: coupling_char, & 156 debug_output, & 157 dt_3d, dz, end_time, initializing_actions, message_string, & 158 origin_date_time, restart_data_format_output, time_since_reference_point, & 159 wind_turbine 160 161 USE cpulog, & 163 162 ONLY: cpu_log, log_point_s 164 163 165 164 USE data_output_module 166 167 USE grid_variables, &165 166 USE grid_variables, & 168 167 ONLY: ddx, dx, ddy, dy 169 168 170 USE indices, &171 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, &169 USE indices, & 170 ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, & 172 171 nzb, nzt, wall_flags_total_0 173 172 174 173 USE kinds 175 174 176 USE netcdf_data_input_mod, &177 ONLY: check_existence, &178 close_input_file, &179 get_variable, &180 input_pids_wtm, &181 inquire_num_variables, &182 inquire_variable_names, &183 input_file_wtm, &184 num_var_pids, &185 open_read_file, &186 pids_id, &175 USE netcdf_data_input_mod, & 176 ONLY: check_existence, & 177 close_input_file, & 178 get_variable, & 179 input_pids_wtm, & 180 inquire_num_variables, & 181 inquire_variable_names, & 182 input_file_wtm, & 183 num_var_pids, & 184 open_read_file, & 185 pids_id, & 187 186 vars_pids 188 187 189 188 USE pegrid 190 189 … … 192 191 ONLY: rrd_mpi_io_global_array, wrd_mpi_io_global_array 193 192 194 193 195 194 IMPLICIT NONE 196 195 197 196 PRIVATE 198 197 199 INTEGER(iwp), PARAMETER :: n_turbines_max = 10000 !< maximum number of turbines (for array allocation)200 201 198 CHARACTER(LEN=100) :: variable_name !< name of output variable 202 CHARACTER(LEN=30) :: nc_filename 203 199 CHARACTER(LEN=30) :: nc_filename !< 200 201 202 INTEGER(iwp), PARAMETER :: n_turbines_max = 1000 !< maximum number of turbines (for array allocation) 203 204 204 205 ! 205 206 !-- Variables specified in the namelist wind_turbine_par 206 207 INTEGER(iwp) :: n_airfoils = 8 !< number of airfoils of the used turbine model (for ADM-R and ALM) 208 INTEGER(iwp) :: n_turbines = 1 !< number of turbines 209 210 211 REAL(wp), DIMENSION(:), POINTER :: output_values_1d_pointer !< pointer for 2d output array 212 REAL(wp), POINTER :: output_values_0d_pointer !< pointer for 2d output array 213 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: output_values_1d_target !< pointer for 2d output array 214 REAL(wp), TARGET :: output_values_0d_target !< pointer for 2d output array 215 216 LOGICAL :: pitch_control = .FALSE. !< switch for use of pitch controller 217 LOGICAL :: speed_control = .FALSE. !< switch for use of speed controller 218 LOGICAL :: yaw_control = .FALSE. !< switch for use of yaw controller 219 LOGICAL :: tip_loss_correction = .FALSE. !< switch for use of tip loss correct. 220 221 LOGICAL :: initial_write_coordinates = .FALSE. 222 223 REAL(wp) :: dt_data_output_wtm = 0.0_wp !< data output interval 224 REAL(wp) :: time_wtm = 0.0_wp !< time since last data output 225 226 207 INTEGER(iwp) :: n_airfoils = 8 !< number of airfoils of the used turbine model (for ADM-R and ALM) 208 INTEGER(iwp) :: n_turbines = 1 !< number of turbines 209 210 LOGICAL :: pitch_control = .FALSE. !< switch for use of pitch controller 211 LOGICAL :: speed_control = .FALSE. !< switch for use of speed controller 212 LOGICAL :: tip_loss_correction = .FALSE. !< switch for use of tip loss correct. 213 LOGICAL :: yaw_control = .FALSE. !< switch for use of yaw controller 214 215 216 LOGICAL :: initial_write_coordinates = .FALSE. !< 217 218 219 REAL(wp), DIMENSION(:), POINTER :: output_values_1d_pointer !< pointer for 2d output array 220 REAL(wp), POINTER :: output_values_0d_pointer !< pointer for 2d output array 221 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: output_values_1d_target !< pointer for 2d output array 222 REAL(wp), TARGET :: output_values_0d_target !< pointer for 2d output array 223 224 225 REAL(wp) :: dt_data_output_wtm = 0.0_wp !< data output interval 226 REAL(wp) :: time_wtm = 0.0_wp !< time since last data output 227 228 227 229 REAL(wp) :: segment_length_tangential = 1.0_wp !< length of the segments, the rotor area is divided into 228 230 !< (in tangential direction, as factor of MIN(dx,dy,dz)) … … 230 232 !< (in radial direction, as factor of MIN(dx,dy,dz)) 231 233 REAL(wp) :: time_turbine_on = 0.0_wp !< time at which turbines are started 232 REAL(wp) :: tilt_angle = 0.0_wp !< vertical tilt_angle of the rotor [degree] ( positive = backwards ) 233 234 235 REAL(wp), DIMENSION(1:n_turbines_max) :: tower_diameter = 0.0_wp !< tower diameter [m] 236 REAL(wp), DIMENSION(1:n_turbines_max), TARGET :: rotor_speed = 0.9_wp !< inital or constant rotor speed [rad/s] 237 REAL(wp), DIMENSION(1:n_turbines_max) :: yaw_angle = 0.0_wp !< yaw angle [degree] ( clockwise, 0 = facing west ) 238 REAL(wp), DIMENSION(1:n_turbines_max) :: pitch_angle = 0.0_wp !< constant pitch angle 234 REAL(wp) :: tilt_angle = 0.0_wp !< vertical tilt_angle of the rotor [degree] ( positive = backwards ) 235 236 !< ( clockwise, 0 = facing west ) 239 237 REAL(wp), DIMENSION(1:n_turbines_max) :: hub_x = 9999999.9_wp !< position of hub in x-direction 240 238 REAL(wp), DIMENSION(1:n_turbines_max) :: hub_y = 9999999.9_wp !< position of hub in y-direction 241 239 REAL(wp), DIMENSION(1:n_turbines_max) :: hub_z = 9999999.9_wp !< position of hub in z-direction 242 240 REAL(wp), DIMENSION(1:n_turbines_max) :: nacelle_radius = 0.0_wp !< nacelle diameter [m] 243 REAL(wp), DIMENSION(1:n_turbines_max) :: rotor_radius = 63.0_wp !< rotor radius [m]244 241 ! REAL(wp), DIMENSION(1:n_turbines_max) :: nacelle_cd = 0.85_wp !< drag coefficient for nacelle 245 242 REAL(wp), DIMENSION(1:n_turbines_max) :: tower_cd = 1.2_wp !< drag coefficient for tower 243 REAL(wp), DIMENSION(1:n_turbines_max) :: pitch_angle = 0.0_wp !< constant pitch angle 244 REAL(wp), DIMENSION(1:n_turbines_max) :: rotor_radius = 63.0_wp !< rotor radius [m] 245 REAL(wp), DIMENSION(1:n_turbines_max), TARGET :: rotor_speed = 0.9_wp !< inital or constant rotor speed [rad/s] 246 REAL(wp), DIMENSION(1:n_turbines_max) :: tower_diameter = 0.0_wp !< tower diameter [m] 247 REAL(wp), DIMENSION(1:n_turbines_max) :: yaw_angle = 0.0_wp !< yaw angle [degree] 248 246 249 247 250 ! 248 251 !-- Variables specified in the namelist for speed controller 249 252 !-- Default values are from the NREL 5MW research turbine (Jonkman, 2008) 250 251 REAL(wp) :: generator_power_rated = 5296610.0_wp !< rated turbine power [W] 252 REAL(wp) :: gear_ratio = 97.0_wp !< Gear ratio from rotor to generator 253 REAL(wp) :: rotor_inertia = 34784179.0_wp !< Inertia of the rotor [kg*m2] 254 REAL(wp) :: generator_inertia = 534.116_wp !< Inertia of the generator [kg*m2] 255 REAL(wp) :: generator_efficiency = 0.944_wp !< Electric efficiency of the generator 256 REAL(wp) :: gear_efficiency = 1.0_wp !< Loss between rotor and generator 257 REAL(wp) :: air_density = 1.225_wp !< Air density to convert to W [kg/m3] 258 REAL(wp) :: generator_speed_rated = 121.6805_wp !< Rated generator speed [rad/s] 259 REAL(wp) :: generator_torque_max = 47402.91_wp !< Maximum of the generator torque [Nm] 260 REAL(wp) :: region_2_slope = 2.332287_wp !< Slope constant for region 2 261 REAL(wp) :: region_2_min = 91.21091_wp !< Lower generator speed boundary of region 2 [rad/s] 262 REAL(wp) :: region_15_min = 70.16224_wp !< Lower generator speed boundary of region 1.5 [rad/s] 263 REAL(wp) :: generator_torque_rate_max = 15000.0_wp !< Max generator torque increase [Nm/s] 264 REAL(wp) :: pitch_rate = 8.0_wp !< Max pitch rate [degree/s] 265 266 267 ! 268 !-- Variables specified in the namelist for yaw control 269 270 REAL(wp) :: yaw_speed = 0.005236_wp !< speed of the yaw actuator [rad/s] 271 REAL(wp) :: yaw_misalignment_max = 0.08726_wp !< maximum tolerated yaw missalignment [rad] 272 REAL(wp) :: yaw_misalignment_min = 0.008726_wp !< minimum yaw missalignment for which the actuator stops [rad] 273 274 275 ! 276 !-- Variables for initialization of the turbine model 277 278 INTEGER(iwp) :: inot !< turbine loop index (turbine id) 279 INTEGER(iwp) :: nsegs_max !< maximum number of segments (all turbines, required for allocation of arrays) 280 INTEGER(iwp) :: nrings_max !< maximum number of rings (all turbines, required for allocation of arrays) 281 INTEGER(iwp) :: ring !< ring loop index (ring number) 282 INTEGER(iwp) :: upper_end !< 283 284 INTEGER(iwp), DIMENSION(1) :: lct !< 285 286 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_hub !< index belonging to x-position of the turbine 287 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_smear !< index defining the area for the smearing of the forces (x-direction) 288 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_hub !< index belonging to y-position of the turbine 289 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_smear !< index defining the area for the smearing of the forces (y-direction) 290 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_hub !< index belonging to hub height 291 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_smear !< index defining the area for the smearing of the forces (z-direction) 292 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nrings !< number of rings per turbine 293 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nsegs_total !< total number of segments per turbine 294 295 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nsegs !< number of segments per ring and turbine 296 297 ! 298 !- parameters for the smearing from the rotor to the cartesian grid 299 REAL(wp) :: pol_a !< parameter for the polynomial smearing fct 300 REAL(wp) :: pol_b !< parameter for the polynomial smearing fct 301 REAL(wp) :: delta_t_factor !< 302 REAL(wp) :: eps_factor !< 303 REAL(wp) :: eps_min !< 304 REAL(wp) :: eps_min2 !< 305 306 ! 307 !-- Variables for the calculation of lift and drag coefficients 308 REAL(wp), DIMENSION(:), ALLOCATABLE :: ard !< 309 REAL(wp), DIMENSION(:), ALLOCATABLE :: crd !< 310 REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_r !< radial segment length 311 REAL(wp), DIMENSION(:), ALLOCATABLE :: lrd !< 312 313 REAL(wp) :: accu_cl_cd_tab = 0.1_wp !< Accuracy of the interpolation of 314 !< the lift and drag coeff [deg] 315 316 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cd_tab !< table of the blade drag coefficient 317 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_tab !< table of the blade lift coefficient 318 319 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: nac_cd_surf !< 3d field of the tower drag coefficient 253 REAL(wp) :: air_density = 1.225_wp !< Air density to convert to W [kg/m3] 254 REAL(wp) :: gear_efficiency = 1.0_wp !< Loss between rotor and generator 255 REAL(wp) :: gear_ratio = 97.0_wp !< Gear ratio from rotor to generator 256 REAL(wp) :: generator_power_rated = 5296610.0_wp !< rated turbine power [W] 257 REAL(wp) :: generator_inertia = 534.116_wp !< Inertia of the generator [kg*m2] 258 REAL(wp) :: generator_efficiency = 0.944_wp !< Electric efficiency of the generator 259 REAL(wp) :: generator_speed_rated = 121.6805_wp !< Rated generator speed [rad/s] 260 REAL(wp) :: generator_torque_max = 47402.91_wp !< Maximum of the generator torque [Nm] 261 REAL(wp) :: generator_torque_rate_max = 15000.0_wp !< Max generator torque increase [Nm/s] 262 REAL(wp) :: pitch_rate = 8.0_wp !< Max pitch rate [degree/s] 263 REAL(wp) :: region_2_slope = 2.332287_wp !< Slope constant for region 2 264 REAL(wp) :: region_2_min = 91.21091_wp !< Lower generator speed boundary of region 2 [rad/s] 265 REAL(wp) :: region_15_min = 70.16224_wp !< Lower generator speed boundary of region 1.5 [rad/s] 266 REAL(wp) :: rotor_inertia = 34784179.0_wp !< Inertia of the rotor [kg*m2] 267 268 269 270 ! 271 !-- Variables specified in the namelist for yaw control: 272 REAL(wp) :: yaw_misalignment_max = 0.08726_wp !< maximum tolerated yaw missalignment [rad] 273 REAL(wp) :: yaw_misalignment_min = 0.008726_wp !< minimum yaw missalignment for which the actuator stops [rad] 274 REAL(wp) :: yaw_speed = 0.005236_wp !< speed of the yaw actuator [rad/s] 275 276 ! 277 !-- Variables for initialization of the turbine model: 278 INTEGER(iwp) :: inot !< turbine loop index (turbine id) 279 INTEGER(iwp) :: nsegs_max !< maximum number of segments (all turbines, required for allocation of arrays) 280 INTEGER(iwp) :: nrings_max !< maximum number of rings (all turbines, required for allocation of arrays) 281 INTEGER(iwp) :: ring !< ring loop index (ring number) 282 INTEGER(iwp) :: upper_end !< 283 284 INTEGER(iwp), DIMENSION(1) :: lct !< 285 286 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_hub !< index belonging to x-position of the turbine 287 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_smear !< index defining the area for the smearing of the forces (x-direction) 288 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_hub !< index belonging to y-position of the turbine 289 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_smear !< index defining the area for the smearing of the forces (y-direction) 290 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_hub !< index belonging to hub height 291 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_smear !< index defining the area for the smearing of the forces (z-direction) 292 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nrings !< number of rings per turbine 293 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nsegs_total !< total number of segments per turbine 294 295 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nsegs !< number of segments per ring and turbine 296 297 ! 298 !- Parameters for the smearing from the rotor to the cartesian grid: 299 REAL(wp) :: delta_t_factor !< 300 REAL(wp) :: eps_factor !< 301 REAL(wp) :: eps_min !< 302 REAL(wp) :: eps_min2 !< 303 REAL(wp) :: pol_a !< parameter for the polynomial smearing fct 304 REAL(wp) :: pol_b !< parameter for the polynomial smearing fct 305 ! 306 !-- Variables for the calculation of lift and drag coefficients: 307 REAL(wp), DIMENSION(:), ALLOCATABLE :: ard !< 308 REAL(wp), DIMENSION(:), ALLOCATABLE :: crd !< 309 REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_r !< radial segment length 310 REAL(wp), DIMENSION(:), ALLOCATABLE :: lrd !< 311 312 REAL(wp) :: accu_cl_cd_tab = 0.1_wp !< Accuracy of the interpolation of the lift and drag coeff [deg] 313 314 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cd_tab !< table of the blade drag coefficient 315 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_tab !< table of the blade lift coefficient 316 317 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: nac_cd_surf !< 3d field of the tower drag coefficient 320 318 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tow_cd_surf !< 3d field of the nacelle drag coefficient 321 319 322 320 ! 323 !-- Variables for the calculation of the forces 324 325 REAL(wp) :: cur_r !< 326 REAL(wp) :: phi_rotor !< 327 REAL(wp) :: pre_factor !< 328 REAL(wp) :: torque_seg !< 329 REAL(wp) :: u_int_l !< 330 REAL(wp) :: u_int_u !< 331 REAL(wp) :: u_rot !< 332 REAL(wp) :: v_int_l !< 333 REAL(wp) :: v_int_u !< 334 REAL(wp) :: w_int_l !< 335 REAL(wp) :: w_int_u !< 321 !-- Variables for the calculation of the forces: 322 REAL(wp) :: cur_r !< 323 REAL(wp) :: phi_rotor !< 324 REAL(wp) :: pre_factor !< 325 REAL(wp) :: torque_seg !< 326 REAL(wp) :: u_int_l !< 327 REAL(wp) :: u_int_u !< 328 REAL(wp) :: u_rot !< 329 REAL(wp) :: v_int_l !< 330 REAL(wp) :: v_int_u !< 331 REAL(wp) :: w_int_l !< 332 REAL(wp) :: w_int_u !< 336 333 !$OMP THREADPRIVATE (cur_r, phi_rotor, pre_factor, torque_seg, u_int_l, u_int_u, u_rot, & 337 334 !$OMP& v_int_l, v_int_u, w_int_l, w_int_u) 338 335 ! 339 !- Tendencies from the nacelle and tower thrust 340 REAL(wp) :: tend_nac_x = 0.0_wp !< 341 REAL(wp) :: tend_ tow_x = 0.0_wp !<342 REAL(wp) :: tend_ nac_y = 0.0_wp !<336 !- Tendencies from the nacelle and tower thrust: 337 REAL(wp) :: tend_nac_x = 0.0_wp !< 338 REAL(wp) :: tend_nac_y = 0.0_wp !< 339 REAL(wp) :: tend_tow_x = 0.0_wp !< 343 340 REAL(wp) :: tend_tow_y = 0.0_wp !< 344 341 !$OMP THREADPRIVATE (tend_nac_x, tend_tow_x, tend_nac_y, tend_tow_y) 345 342 346 REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_attack !<347 REAL(wp), DIMENSION(:), ALLOCATABLE :: chord !<348 REAL(wp), DIMENSION(:), ALLOCATABLE :: phi_rel !<349 REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_total !<350 REAL(wp), DIMENSION(:), ALLOCATABLE :: thrust_rotor !<351 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_c l !<352 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_c d !<353 REAL(wp), DIMENSION(:), ALLOCATABLE :: vrel !<354 REAL(wp), DIMENSION(:), ALLOCATABLE :: vtheta !<343 REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_attack !< 344 REAL(wp), DIMENSION(:), ALLOCATABLE :: chord !< 345 REAL(wp), DIMENSION(:), ALLOCATABLE :: phi_rel !< 346 REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_total !< 347 REAL(wp), DIMENSION(:), ALLOCATABLE :: thrust_rotor !< 348 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd !< 349 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl !< 350 REAL(wp), DIMENSION(:), ALLOCATABLE :: vrel !< 351 REAL(wp), DIMENSION(:), ALLOCATABLE :: vtheta !< 355 352 356 353 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rbx, rby, rbz !< coordinates of the blade elements … … 358 355 359 356 ! 360 !- Fields for the interpolation of velocities on the rotor grid 361 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_int !< 362 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_int_1_l !< 363 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_int !< 364 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_int_1_l !< 365 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_int !< 366 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_int_1_l !< 367 368 ! 369 !- rotor tendencies on the segments 370 REAL(wp), DIMENSION(:), ALLOCATABLE :: thrust_seg !< 371 REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_seg_y !< 372 REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_seg_z !< 373 374 ! 375 !- rotor tendencies on the rings 376 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: thrust_ring !< 377 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: torque_ring_y !< 378 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: torque_ring_z !< 379 380 ! 381 !- rotor tendencies on rotor grids for all turbines 382 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: thrust !< 383 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: torque_y !< 384 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: torque_z !< 385 386 ! 387 !- rotor tendencies on coordinate grid 388 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_x !< 389 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_y !< 390 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_z !< 391 ! 392 !- variables for the rotation of the rotor coordinates 357 !- Fields for the interpolation of velocities on the rotor grid: 358 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_int !< 359 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_int_1_l !< 360 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_int !< 361 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_int_1_l !< 362 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_int !< 363 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_int_1_l !< 364 365 ! 366 !- Rotor tendencies on the segments: 367 REAL(wp), DIMENSION(:), ALLOCATABLE :: thrust_seg !< 368 REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_seg_y !< 369 REAL(wp), DIMENSION(:), ALLOCATABLE :: torque_seg_z !< 370 371 ! 372 !- Rotor tendencies on the rings: 373 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: thrust_ring !< 374 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: torque_ring_y !< 375 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: torque_ring_z !< 376 377 ! 378 !- Rotor tendencies on rotor grids for all turbines: 379 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: thrust !< 380 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: torque_y !< 381 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: torque_z !< 382 383 ! 384 !- Rotor tendencies on coordinate grid: 385 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_x !< 386 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_y !< 387 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rot_tend_z !< 388 389 ! 390 !- Variables for the rotation of the rotor coordinates: 393 391 REAL(wp), DIMENSION(1:n_turbines_max,1:3,1:3) :: rot_coord_trans !< matrix for rotation of rotor coordinates 394 395 REAL(wp), DIMENSION(1:3) :: rot_eigen_rad !< 396 REAL(wp), DIMENSION(1:3) :: rot_eigen_azi !< 397 REAL(wp), DIMENSION(1:3) :: rot_eigen_nor !< 398 REAL(wp), DIMENSION(1:3) :: re !< 399 REAL(wp), DIMENSION(1:3) :: rea !< 400 REAL(wp), DIMENSION(1:3) :: ren !< 401 REAL(wp), DIMENSION(1:3) :: rote !< 402 REAL(wp), DIMENSION(1:3) :: rota !< 403 REAL(wp), DIMENSION(1:3) :: rotn !< 404 405 ! 406 !-- Fixed variables for the speed controller 407 408 LOGICAL :: start_up = .TRUE. !< 409 410 REAL(wp) :: Fcorner !< corner freq for the controller low pass filter 411 REAL(wp) :: region_25_min !< min region 2.5 412 REAL(wp) :: om_rate !< rotor speed change 413 REAL(wp) :: slope15 !< slope in region 1.5 414 REAL(wp) :: region_25_slope !< slope in region 2.5 415 REAL(wp) :: trq_rate !< torque change 416 REAL(wp) :: vs_sysp !< 417 REAL(wp) :: lp_coeff !< coeff for the controller low pass filter 418 419 REAL(wp), DIMENSION(n_turbines_max) :: rotor_speed_l = 0.0_wp !< local rot speed [rad/s] 420 421 ! 422 !-- Fixed variables for the yaw controller 423 392 393 REAL(wp), DIMENSION(1:3) :: rot_eigen_rad !< 394 REAL(wp), DIMENSION(1:3) :: rot_eigen_azi !< 395 REAL(wp), DIMENSION(1:3) :: rot_eigen_nor !< 396 REAL(wp), DIMENSION(1:3) :: re !< 397 REAL(wp), DIMENSION(1:3) :: rea !< 398 REAL(wp), DIMENSION(1:3) :: ren !< 399 REAL(wp), DIMENSION(1:3) :: rote !< 400 REAL(wp), DIMENSION(1:3) :: rota !< 401 REAL(wp), DIMENSION(1:3) :: rotn !< 402 403 ! 404 !-- Fixed variables for the speed controller: 405 LOGICAL :: start_up = .TRUE. !< 406 407 REAL(wp) :: fcorner !< corner freq for the controller low pass filter 408 REAL(wp) :: om_rate !< rotor speed change 409 REAL(wp) :: region_25_min !< min region 2.5 410 REAL(wp) :: region_25_slope !< slope in region 2.5 411 REAL(wp) :: slope15 !< slope in region 1.5 412 REAL(wp) :: trq_rate !< torque change 413 REAL(wp) :: vs_sysp !< 414 REAL(wp) :: lp_coeff !< coeff for the controller low pass filter 415 416 REAL(wp), DIMENSION(n_turbines_max) :: rotor_speed_l = 0.0_wp !< local rot speed [rad/s] 417 418 ! 419 !-- Fixed variables for the yaw controller: 420 INTEGER(iwp) :: wdlon !< 421 INTEGER(iwp) :: wdsho !< 422 423 LOGICAL, DIMENSION(1:n_turbines_max) :: doyaw = .FALSE. !< 424 425 REAL(wp), DIMENSION(:) , ALLOCATABLE :: u_inflow !< wind speed at hub 426 REAL(wp), DIMENSION(:) , ALLOCATABLE :: u_inflow_l !< 427 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wdir !< wind direction at hub 428 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wdir_l !< 429 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd2_l !< local (cpu) short running avg of the wd 430 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd30_l !< local (cpu) long running avg of the wd 424 431 REAL(wp), DIMENSION(:) , ALLOCATABLE :: yawdir !< direction to yaw 425 432 REAL(wp), DIMENSION(:) , ALLOCATABLE :: yaw_angle_l !< local (cpu) yaw angle 426 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd30_l !< local (cpu) long running avg of the wd 427 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wd2_l !< local (cpu) short running avg of the wd 428 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wdir !< wind direction at hub 429 REAL(wp), DIMENSION(:) , ALLOCATABLE :: u_inflow !< wind speed at hub 430 REAL(wp), DIMENSION(:) , ALLOCATABLE :: wdir_l !< 431 REAL(wp), DIMENSION(:) , ALLOCATABLE :: u_inflow_l !< 433 434 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd2 !< 432 435 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd30 !< 433 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd2 !< 434 LOGICAL, DIMENSION(1:n_turbines_max) :: doyaw = .FALSE. !< 435 INTEGER(iwp) :: WDLON !< 436 INTEGER(iwp) :: WDSHO !< 437 438 ! 439 !-- Variables that have to be saved in the binary file for restarts 440 REAL(wp), DIMENSION(1:n_turbines_max) :: pitch_angle_old = 0.0_wp !< old constant pitch angle 436 437 438 ! 439 !-- Variables that have to be saved in the binary file for restarts: 440 REAL(wp), DIMENSION(1:n_turbines_max) :: pitch_angle_old = 0.0_wp !< old constant pitch angle 441 441 REAL(wp), DIMENSION(1:n_turbines_max) :: generator_speed = 0.0_wp !< curr. generator speed 442 442 REAL(wp), DIMENSION(1:n_turbines_max) :: generator_speed_f = 0.0_wp !< filtered generator speed … … 461 461 MODULE PROCEDURE wtm_data_output 462 462 END INTERFACE wtm_data_output 463 463 464 464 INTERFACE wtm_init_arrays 465 465 MODULE PROCEDURE wtm_init_arrays … … 473 473 MODULE PROCEDURE wtm_init_output 474 474 END INTERFACE wtm_init_output 475 475 476 476 INTERFACE wtm_actions 477 477 MODULE PROCEDURE wtm_actions … … 489 489 490 490 491 PUBLIC &492 dt_data_output_wtm, &493 time_wtm, &491 PUBLIC & 492 dt_data_output_wtm, & 493 time_wtm, & 494 494 wind_turbine 495 495 496 PUBLIC &497 wtm_parin, &498 wtm_check_parameters, &499 wtm_data_output, &500 wtm_init_arrays, &501 wtm_init_output, &502 wtm_init, &503 wtm_actions, &504 wtm_rrd_global, &496 PUBLIC & 497 wtm_parin, & 498 wtm_check_parameters, & 499 wtm_data_output, & 500 wtm_init_arrays, & 501 wtm_init_output, & 502 wtm_init, & 503 wtm_actions, & 504 wtm_rrd_global, & 505 505 wtm_wrd_global 506 506 … … 509 509 510 510 511 !------------------------------------------------------------------------------ !511 !--------------------------------------------------------------------------------------------------! 512 512 ! Description: 513 513 ! ------------ 514 514 !> Parin for &wind_turbine_par for wind turbine model 515 !------------------------------------------------------------------------------ !516 517 518 519 520 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file521 522 NAMELIST /wind_turbine_parameters/&523 air_density, tower_diameter, dt_data_output_wtm,&524 gear_efficiency, gear_ratio, generator_efficiency,&525 generator_inertia, rotor_inertia, yaw_misalignment_max,&526 generator_torque_max, generator_torque_rate_max, yaw_misalignment_min,&527 region_15_min, region_2_min, n_airfoils, n_turbines,&528 rotor_speed, yaw_angle, pitch_angle, pitch_control,&529 generator_speed_rated, generator_power_rated, hub_x, hub_y, hub_z,&530 nacelle_radius, rotor_radius, segment_length_tangential, segment_width_radial,&531 region_2_slope, speed_control, tilt_angle, time_turbine_on,&532 tower_cd, pitch_rate,&533 534 ! 535 536 ! 537 !-- Try to find wind turbine model package538 539 540 DO WHILE ( INDEX( line, '&wind_turbine_parameters' ) == 0 )541 READ ( 11, '(A)', END=12 ) line542 543 544 545 ! 546 !-- Read user-defined namelist547 548 ! 549 !-- Set flag that indicates that the wind turbine model is switched on550 551 552 553 554 10 BACKSPACE( 11 )555 READ( 11 , '(A)') line556 557 558 12 CONTINUE! TBD Change from continue, mit ierrn machen559 560 561 562 563 !------------------------------------------------------------------------------ !515 !--------------------------------------------------------------------------------------------------! 516 SUBROUTINE wtm_parin 517 518 IMPLICIT NONE 519 520 CHARACTER(LEN=80) :: line !< dummy string that contains the current line of the parameter file 521 522 NAMELIST /wind_turbine_parameters/ & 523 air_density, tower_diameter, dt_data_output_wtm, & 524 gear_efficiency, gear_ratio, generator_efficiency, & 525 generator_inertia, rotor_inertia, yaw_misalignment_max, & 526 generator_torque_max, generator_torque_rate_max, yaw_misalignment_min, & 527 region_15_min, region_2_min, n_airfoils, n_turbines, & 528 rotor_speed, yaw_angle, pitch_angle, pitch_control, & 529 generator_speed_rated, generator_power_rated, hub_x, hub_y, hub_z, & 530 nacelle_radius, rotor_radius, segment_length_tangential, segment_width_radial, & 531 region_2_slope, speed_control, tilt_angle, time_turbine_on, & 532 tower_cd, pitch_rate, & 533 yaw_control, yaw_speed, tip_loss_correction 534 ! , nacelle_cd 535 536 ! 537 !-- Try to find wind turbine model package: 538 REWIND ( 11 ) 539 line = ' ' 540 DO WHILE ( INDEX ( line, '&wind_turbine_parameters' ) == 0 ) 541 READ ( 11, '(A)', END = 12 ) line 542 ENDDO 543 BACKSPACE ( 11 ) 544 545 ! 546 !-- Read user-defined namelist: 547 READ ( 11, wind_turbine_parameters, ERR = 10, END = 12 ) 548 ! 549 !-- Set flag that indicates that the wind turbine model is switched on: 550 wind_turbine = .TRUE. 551 552 GOTO 12 553 554 10 BACKSPACE ( 11 ) 555 READ ( 11, '(A)' ) line 556 CALL parin_fail_message( 'wind_turbine_parameters', line ) 557 558 12 CONTINUE ! TBD Change from continue, mit ierrn machen 559 560 END SUBROUTINE wtm_parin 561 562 563 !--------------------------------------------------------------------------------------------------! 564 564 ! Description: 565 565 ! ------------ 566 566 !> This routine writes the respective restart data. 567 !------------------------------------------------------------------------------ !568 SUBROUTINE wtm_wrd_global569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 END SUBROUTINE wtm_wrd_global622 623 624 !------------------------------------------------------------------------------ !567 !--------------------------------------------------------------------------------------------------! 568 SUBROUTINE wtm_wrd_global 569 570 571 IMPLICIT NONE 572 573 574 IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN 575 576 CALL wrd_write_string( 'generator_speed' ) 577 WRITE ( 14 ) generator_speed 578 579 CALL wrd_write_string( 'generator_speed_f' ) 580 WRITE ( 14 ) generator_speed_f 581 582 CALL wrd_write_string( 'generator_speed_f_old' ) 583 WRITE ( 14 ) generator_speed_f_old 584 585 CALL wrd_write_string( 'generator_speed_old' ) 586 WRITE ( 14 ) generator_speed_old 587 588 CALL wrd_write_string( 'rotor_speed' ) 589 WRITE ( 14 ) rotor_speed 590 591 CALL wrd_write_string( 'yaw_angle' ) 592 WRITE ( 14 ) yaw_angle 593 594 CALL wrd_write_string( 'pitch_angle' ) 595 WRITE ( 14 ) pitch_angle 596 597 CALL wrd_write_string( 'pitch_angle_old' ) 598 WRITE ( 14 ) pitch_angle_old 599 600 CALL wrd_write_string( 'torque_gen' ) 601 WRITE ( 14 ) torque_gen 602 603 CALL wrd_write_string( 'torque_gen_old' ) 604 WRITE ( 14 ) torque_gen_old 605 606 ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' ) THEN 607 608 CALL wrd_mpi_io_global_array( 'generator_speed', generator_speed ) 609 CALL wrd_mpi_io_global_array( 'generator_speed_f', generator_speed_f ) 610 CALL wrd_mpi_io_global_array( 'generator_speed_f_old', generator_speed_f_old ) 611 CALL wrd_mpi_io_global_array( 'generator_speed_old', generator_speed_old ) 612 CALL wrd_mpi_io_global_array( 'rotor_speed', rotor_speed ) 613 CALL wrd_mpi_io_global_array( 'yaw_angle', yaw_angle ) 614 CALL wrd_mpi_io_global_array( 'pitch_angle', pitch_angle ) 615 CALL wrd_mpi_io_global_array( 'pitch_angle_old', pitch_angle_old ) 616 CALL wrd_mpi_io_global_array( 'torque_gen', torque_gen ) 617 CALL wrd_mpi_io_global_array( 'torque_gen_old', torque_gen_old ) 618 619 ENDIF 620 621 END SUBROUTINE wtm_wrd_global 622 623 624 !--------------------------------------------------------------------------------------------------! 625 625 ! Description: 626 626 ! ------------ 627 627 !> Read module-specific global restart data (Fortran binary format). 628 !------------------------------------------------------------------------------ !628 !--------------------------------------------------------------------------------------------------! 629 629 SUBROUTINE wtm_rrd_global_ftn( found ) 630 630 631 631 632 USE control_parameters, &632 USE control_parameters, & 633 633 ONLY: length, restart_string 634 634 … … 636 636 IMPLICIT NONE 637 637 638 LOGICAL, INTENT(OUT) :: found638 LOGICAL, INTENT(OUT) :: found 639 639 640 640 … … 670 670 671 671 END SELECT 672 672 673 673 674 674 END SUBROUTINE wtm_rrd_global_ftn … … 696 696 697 697 698 !------------------------------------------------------------------------------ !698 !--------------------------------------------------------------------------------------------------! 699 699 ! Description: 700 700 ! ------------ 701 701 !> Check namelist parameter 702 !------------------------------------------------------------------------------! 703 SUBROUTINE wtm_check_parameters 704 705 IMPLICIT NONE 706 707 IF ( .NOT. input_pids_wtm ) THEN 708 IF ( ( .NOT.speed_control ) .AND. pitch_control ) THEN 709 message_string = 'pitch_control = .TRUE. requires '// & 710 'speed_control = .TRUE.' 711 CALL message( 'wtm_check_parameters', 'PA0461', 1, 2, 0, 6, 0 ) 712 ENDIF 713 714 IF ( ANY( rotor_speed(1:n_turbines) < 0.0 ) ) THEN 715 message_string = 'rotor_speed < 0.0, Please set rotor_speed to ' // & 716 'a value equal or larger than zero' 717 CALL message( 'wtm_check_parameters', 'PA0462', 1, 2, 0, 6, 0 ) 718 ENDIF 719 720 721 IF ( ANY( hub_x(1:n_turbines) == 9999999.9_wp ) .OR. & 722 ANY( hub_y(1:n_turbines) == 9999999.9_wp ) .OR. & 723 ANY( hub_z(1:n_turbines) == 9999999.9_wp ) ) THEN 724 725 message_string = 'hub_x, hub_y, hub_z ' // & 726 'have to be given for each turbine.' 727 CALL message( 'wtm_check_parameters', 'PA0463', 1, 2, 0, 6, 0 ) 728 ENDIF 702 !--------------------------------------------------------------------------------------------------! 703 SUBROUTINE wtm_check_parameters 704 705 IMPLICIT NONE 706 707 IF ( .NOT. input_pids_wtm ) THEN 708 IF ( ( .NOT.speed_control ) .AND. pitch_control ) THEN 709 message_string = 'pitch_control = .TRUE. requires speed_control = .TRUE.' 710 CALL message( 'wtm_check_parameters', 'PA0461', 1, 2, 0, 6, 0 ) 729 711 ENDIF 730 731 END SUBROUTINE wtm_check_parameters 732 ! 733 734 !------------------------------------------------------------------------------! 712 713 IF ( ANY( rotor_speed(1:n_turbines) < 0.0 ) ) THEN 714 message_string = 'rotor_speed < 0.0, Please set rotor_speed to ' // & 715 'a value equal or larger than zero' 716 CALL message( 'wtm_check_parameters', 'PA0462', 1, 2, 0, 6, 0 ) 717 ENDIF 718 719 IF ( ANY( hub_x(1:n_turbines) == 9999999.9_wp ) .OR. & 720 ANY( hub_y(1:n_turbines) == 9999999.9_wp ) .OR. & 721 ANY( hub_z(1:n_turbines) == 9999999.9_wp ) ) THEN 722 723 message_string = 'hub_x, hub_y, hub_z have to be given for each turbine.' 724 CALL message( 'wtm_check_parameters', 'PA0463', 1, 2, 0, 6, 0 ) 725 ENDIF 726 ENDIF 727 728 END SUBROUTINE wtm_check_parameters 729 ! 730 !--------------------------------------------------------------------------------------------------! 735 731 ! Description: 736 732 ! ------------ 737 733 !> Allocate wind turbine model arrays 738 !------------------------------------------------------------------------------ !739 740 741 742 743 REAL(wp) :: delta_r_factor !<744 REAL(wp) :: delta_r_init !<734 !--------------------------------------------------------------------------------------------------! 735 SUBROUTINE wtm_init_arrays 736 737 IMPLICIT NONE 738 739 REAL(wp) :: delta_r_factor !< 740 REAL(wp) :: delta_r_init !< 745 741 746 742 #if defined( __netcdf ) 747 743 ! 748 ! Read wtm input file (netcdf) if it exists 749 IF ( input_pids_wtm ) THEN 750 751 ! 752 !-- Open the wtm input file 753 CALL open_read_file( TRIM( input_file_wtm ) // & 754 TRIM( coupling_char ), pids_id ) 755 756 CALL inquire_num_variables( pids_id, num_var_pids ) 757 758 ! 759 !-- Allocate memory to store variable names and read them 760 ALLOCATE( vars_pids(1:num_var_pids) ) 761 CALL inquire_variable_names( pids_id, vars_pids ) 762 763 ! 764 !-- Input of all wtm parameters 765 IF ( check_existence( vars_pids, 'tower_diameter' ) ) THEN 766 CALL get_variable( pids_id, 'tower_diameter', tower_diameter(1:n_turbines)) 767 ENDIF 768 769 IF ( check_existence( vars_pids, 'rotor_speed' ) ) THEN 770 CALL get_variable( pids_id, 'rotor_speed', rotor_speed(1:n_turbines)) 771 ENDIF 772 773 IF ( check_existence( vars_pids, 'pitch_angle' ) ) THEN 774 CALL get_variable( pids_id, 'pitch_angle', pitch_angle(1:n_turbines)) 775 ENDIF 776 777 IF ( check_existence( vars_pids, 'yaw_angle' ) ) THEN 778 CALL get_variable( pids_id, 'yaw_angle', yaw_angle(1:n_turbines)) 779 ENDIF 780 781 IF ( check_existence( vars_pids, 'hub_x' ) ) THEN 782 CALL get_variable( pids_id, 'hub_x', hub_x(1:n_turbines)) 783 ENDIF 784 785 IF ( check_existence( vars_pids, 'hub_y' ) ) THEN 786 CALL get_variable( pids_id, 'hub_y', hub_y(1:n_turbines)) 787 ENDIF 788 789 IF ( check_existence( vars_pids, 'hub_z' ) ) THEN 790 CALL get_variable( pids_id, 'hub_z', hub_z(1:n_turbines)) 791 ENDIF 792 793 IF ( check_existence( vars_pids, 'nacelle_radius' ) ) THEN 794 CALL get_variable( pids_id, 'nacelle_radius', nacelle_radius(1:n_turbines)) 795 ENDIF 796 797 IF ( check_existence( vars_pids, 'rotor_radius' ) ) THEN 798 CALL get_variable( pids_id, 'rotor_radius', rotor_radius(1:n_turbines)) 799 ENDIF 800 ! 801 ! IF ( check_existence( vars_pids, 'nacelle_cd' ) ) THEN 802 ! CALL get_variable( pids_id, 'nacelle_cd', nacelle_cd(1:n_turbines)) 803 ! ENDIF 804 805 IF ( check_existence( vars_pids, 'tower_cd' ) ) THEN 806 CALL get_variable( pids_id, 'tower_cd', tower_cd(1:n_turbines)) 807 ENDIF 808 ! 809 !-- Close wtm input file 810 CALL close_input_file( pids_id ) 811 744 ! Read wtm input file (netcdf) if it exists: 745 IF ( input_pids_wtm ) THEN 746 747 ! 748 !-- Open the wtm input file: 749 CALL open_read_file( TRIM( input_file_wtm ) // & 750 TRIM( coupling_char ), pids_id ) 751 752 CALL inquire_num_variables( pids_id, num_var_pids ) 753 754 ! 755 !-- Allocate memory to store variable names and read them: 756 ALLOCATE( vars_pids(1:num_var_pids) ) 757 CALL inquire_variable_names( pids_id, vars_pids ) 758 759 ! 760 !-- Input of all wtm parameters: 761 IF ( check_existence( vars_pids, 'tower_diameter' ) ) THEN 762 CALL get_variable( pids_id, 'tower_diameter', tower_diameter(1:n_turbines) ) 812 763 ENDIF 764 765 IF ( check_existence( vars_pids, 'rotor_speed' ) ) THEN 766 CALL get_variable( pids_id, 'rotor_speed', rotor_speed(1:n_turbines) ) 767 ENDIF 768 769 IF ( check_existence( vars_pids, 'pitch_angle' ) ) THEN 770 CALL get_variable( pids_id, 'pitch_angle', pitch_angle(1:n_turbines) ) 771 ENDIF 772 773 IF ( check_existence( vars_pids, 'yaw_angle' ) ) THEN 774 CALL get_variable( pids_id, 'yaw_angle', yaw_angle(1:n_turbines) ) 775 ENDIF 776 777 IF ( check_existence( vars_pids, 'hub_x' ) ) THEN 778 CALL get_variable( pids_id, 'hub_x', hub_x(1:n_turbines) ) 779 ENDIF 780 781 IF ( check_existence( vars_pids, 'hub_y' ) ) THEN 782 CALL get_variable( pids_id, 'hub_y', hub_y(1:n_turbines) ) 783 ENDIF 784 785 IF ( check_existence( vars_pids, 'hub_z' ) ) THEN 786 CALL get_variable( pids_id, 'hub_z', hub_z(1:n_turbines) ) 787 ENDIF 788 789 IF ( check_existence( vars_pids, 'nacelle_radius' ) ) THEN 790 CALL get_variable( pids_id, 'nacelle_radius', nacelle_radius(1:n_turbines) ) 791 ENDIF 792 793 IF ( check_existence( vars_pids, 'rotor_radius' ) ) THEN 794 CALL get_variable( pids_id, 'rotor_radius', rotor_radius(1:n_turbines) ) 795 ENDIF 796 ! 797 ! IF ( check_existence( vars_pids, 'nacelle_cd' ) ) THEN 798 ! CALL get_variable( pids_id, 'nacelle_cd', nacelle_cd(1:n_turbines) ) 799 ! ENDIF 800 801 IF ( check_existence( vars_pids, 'tower_cd' ) ) THEN 802 CALL get_variable( pids_id, 'tower_cd', tower_cd(1:n_turbines) ) 803 ENDIF 804 ! 805 !-- Close wtm input file: 806 CALL close_input_file( pids_id ) 807 808 ENDIF 813 809 #endif 814 810 815 811 ! 816 !-- To be able to allocate arrays with dimension of rotor rings and segments, 817 !-- the maximum possible numbers of rings and segments have to be calculated: 818 819 ALLOCATE( nrings(1:n_turbines) ) 820 ALLOCATE( delta_r(1:n_turbines) ) 821 822 nrings(:) = 0 823 delta_r(:) = 0.0_wp 824 825 ! 826 !-- Thickness (radial) of each ring and length (tangential) of each segment: 827 delta_r_factor = segment_width_radial 828 delta_t_factor = segment_length_tangential 829 delta_r_init = delta_r_factor * MIN( dx, dy, dz(1)) 830 831 DO inot = 1, n_turbines 832 ! 833 !-- Determine number of rings: 834 nrings(inot) = NINT( rotor_radius(inot) / delta_r_init ) 835 836 delta_r(inot) = rotor_radius(inot) / nrings(inot) 837 812 !-- To be able to allocate arrays with dimension of rotor rings and segments, 813 !-- the maximum possible numbers of rings and segments have to be calculated: 814 ALLOCATE( nrings(1:n_turbines) ) 815 ALLOCATE( delta_r(1:n_turbines) ) 816 817 nrings(:) = 0 818 delta_r(:) = 0.0_wp 819 820 ! 821 !-- Thickness (radial) of each ring and length (tangential) of each segment: 822 delta_r_factor = segment_width_radial 823 delta_t_factor = segment_length_tangential 824 delta_r_init = delta_r_factor * MIN( dx, dy, dz(1) ) 825 826 DO inot = 1, n_turbines 827 ! 828 !-- Determine number of rings: 829 nrings(inot) = NINT( rotor_radius(inot) / delta_r_init ) 830 831 delta_r(inot) = rotor_radius(inot) / nrings(inot) 832 833 ENDDO 834 835 nrings_max = MAXVAL( nrings ) 836 837 ALLOCATE( nsegs(1:nrings_max,1:n_turbines) ) 838 ALLOCATE( nsegs_total(1:n_turbines) ) 839 840 nsegs(:,:) = 0 841 nsegs_total(:) = 0 842 843 844 DO inot = 1, n_turbines 845 DO ring = 1, nrings(inot) 846 ! 847 !-- Determine number of segments for each ring: 848 nsegs(ring,inot) = MAX( 8, CEILING( delta_r_factor * pi * ( 2.0_wp * ring - 1.0_wp ) / & 849 delta_t_factor ) ) 838 850 ENDDO 839 840 nrings_max = MAXVAL(nrings) 841 842 ALLOCATE( nsegs(1:nrings_max,1:n_turbines) ) 843 ALLOCATE( nsegs_total(1:n_turbines) ) 844 845 nsegs(:,:) = 0 846 nsegs_total(:) = 0 847 848 849 DO inot = 1, n_turbines 850 DO ring = 1, nrings(inot) 851 ! 852 !-- Determine number of segments for each ring: 853 nsegs(ring,inot) = MAX( 8, CEILING( delta_r_factor * pi * & 854 ( 2.0_wp * ring - 1.0_wp ) / & 855 delta_t_factor ) ) 856 ENDDO 857 ! 858 !-- Total sum of all rotor segments: 859 nsegs_total(inot) = SUM( nsegs(:,inot) ) 860 861 ENDDO 862 863 ! 864 !-- Maximum number of segments per ring: 865 nsegs_max = MAXVAL(nsegs) 851 ! 852 !-- Total sum of all rotor segments: 853 nsegs_total(inot) = SUM( nsegs(:,inot) ) 854 ENDDO 855 856 ! 857 !-- Maximum number of segments per ring: 858 nsegs_max = MAXVAL( nsegs ) 866 859 867 860 !! 868 !!-- 869 ! 870 ! 871 ! PRINT*, '--------------------------------------------------'872 ! 873 ! PRINT*, '--------------------------------------------------'874 ! 875 ! 876 ! 877 ! 878 879 880 ! 881 !-- Allocate 1D arrays (dimension = number of turbines)882 883 884 885 886 887 888 889 890 891 ! 892 !-- Allocation of the 1D arrays for yaw control893 894 895 896 897 898 899 900 ! 901 !-- Allocate 1D arrays (dimension = number of rotor segments)902 903 904 905 906 907 908 909 910 911 912 913 ! 914 !-- Allocate 2D arrays (dimension = number of rotor rings and segments)915 916 917 918 919 920 921 922 ! 923 !-- Allocate additional 2D arrays924 925 926 927 928 ! 929 !-- Allocate 3D arrays (dimension = number of grid points)930 931 932 933 934 935 936 937 938 939 ! 940 !-- Allocate additional 3D arrays941 942 943 944 945 946 947 948 ! 949 !-- 950 951 952 953 954 955 956 957 958 959 960 IF ( TRIM( initializing_actions ) /= 'read_restart_data' )THEN961 962 963 964 965 pitch_angle_old(:)= 0.0_wp966 torque_gen(:)= 0.0_wp967 torque_gen_old(:)= 0.0_wp968 969 970 971 972 973 974 975 976 977 ! 978 !-- Allocate 1D arrays (dimension = number of rotor segments)979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 !------------------------------------------------------------------------------ !861 !!-- TODO: Folgendes im Header ausgeben! 862 ! IF ( myid == 0 ) THEN 863 ! PRINT*, 'nrings(1) = ', nrings(1) 864 ! PRINT*, '----------------------------------------------------------------------------------' 865 ! PRINT*, 'nsegs(:,1) = ', nsegs(:,1) 866 ! PRINT*, '----------------------------------------------------------------------------------' 867 ! PRINT*, 'nrings_max = ', nrings_max 868 ! PRINT*, 'nsegs_max = ', nsegs_max 869 ! PRINT*, 'nsegs_total(1) = ', nsegs_total(1) 870 ! ENDIF 871 872 873 ! 874 !-- Allocate 1D arrays (dimension = number of turbines): 875 ALLOCATE( i_hub(1:n_turbines) ) 876 ALLOCATE( i_smear(1:n_turbines) ) 877 ALLOCATE( j_hub(1:n_turbines) ) 878 ALLOCATE( j_smear(1:n_turbines) ) 879 ALLOCATE( k_hub(1:n_turbines) ) 880 ALLOCATE( k_smear(1:n_turbines) ) 881 ALLOCATE( torque_total(1:n_turbines) ) 882 ALLOCATE( thrust_rotor(1:n_turbines) ) 883 884 ! 885 !-- Allocation of the 1D arrays for yaw control: 886 ALLOCATE( yawdir(1:n_turbines) ) 887 ALLOCATE( u_inflow(1:n_turbines) ) 888 ALLOCATE( wdir(1:n_turbines) ) 889 ALLOCATE( u_inflow_l(1:n_turbines) ) 890 ALLOCATE( wdir_l(1:n_turbines) ) 891 ALLOCATE( yaw_angle_l(1:n_turbines) ) 892 893 ! 894 !-- Allocate 1D arrays (dimension = number of rotor segments): 895 ALLOCATE( alpha_attack(1:nsegs_max) ) 896 ALLOCATE( chord(1:nsegs_max) ) 897 ALLOCATE( phi_rel(1:nsegs_max) ) 898 ALLOCATE( thrust_seg(1:nsegs_max) ) 899 ALLOCATE( torque_seg_y(1:nsegs_max) ) 900 ALLOCATE( torque_seg_z(1:nsegs_max) ) 901 ALLOCATE( turb_cd(1:nsegs_max) ) 902 ALLOCATE( turb_cl(1:nsegs_max) ) 903 ALLOCATE( vrel(1:nsegs_max) ) 904 ALLOCATE( vtheta(1:nsegs_max) ) 905 906 ! 907 !-- Allocate 2D arrays (dimension = number of rotor rings and segments): 908 ALLOCATE( rbx(1:nrings_max,1:nsegs_max) ) 909 ALLOCATE( rby(1:nrings_max,1:nsegs_max) ) 910 ALLOCATE( rbz(1:nrings_max,1:nsegs_max) ) 911 ALLOCATE( thrust_ring(1:nrings_max,1:nsegs_max) ) 912 ALLOCATE( torque_ring_y(1:nrings_max,1:nsegs_max) ) 913 ALLOCATE( torque_ring_z(1:nrings_max,1:nsegs_max) ) 914 915 ! 916 !-- Allocate additional 2D arrays: 917 ALLOCATE( rotx(1:n_turbines,1:3) ) 918 ALLOCATE( roty(1:n_turbines,1:3) ) 919 ALLOCATE( rotz(1:n_turbines,1:3) ) 920 921 ! 922 !-- Allocate 3D arrays (dimension = number of grid points): 923 ALLOCATE( nac_cd_surf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 924 ALLOCATE( rot_tend_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 925 ALLOCATE( rot_tend_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 926 ALLOCATE( rot_tend_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 927 ALLOCATE( thrust(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 928 ALLOCATE( torque_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 929 ALLOCATE( torque_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 930 ALLOCATE( tow_cd_surf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 931 932 ! 933 !-- Allocate additional 3D arrays: 934 ALLOCATE( u_int(1:n_turbines,1:nrings_max,1:nsegs_max) ) 935 ALLOCATE( u_int_1_l(1:n_turbines,1:nrings_max,1:nsegs_max) ) 936 ALLOCATE( v_int(1:n_turbines,1:nrings_max,1:nsegs_max) ) 937 ALLOCATE( v_int_1_l(1:n_turbines,1:nrings_max,1:nsegs_max) ) 938 ALLOCATE( w_int(1:n_turbines,1:nrings_max,1:nsegs_max) ) 939 ALLOCATE( w_int_1_l(1:n_turbines,1:nrings_max,1:nsegs_max) ) 940 941 ! 942 !-- All of the arrays are initialized with a value of zero: 943 i_hub(:) = 0 944 i_smear(:) = 0 945 j_hub(:) = 0 946 j_smear(:) = 0 947 k_hub(:) = 0 948 k_smear(:) = 0 949 950 torque_total(:) = 0.0_wp 951 thrust_rotor(:) = 0.0_wp 952 953 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 954 generator_speed(:) = 0.0_wp 955 generator_speed_old(:) = 0.0_wp 956 generator_speed_f(:) = 0.0_wp 957 generator_speed_f_old(:) = 0.0_wp 958 pitch_angle_old(:) = 0.0_wp 959 torque_gen(:) = 0.0_wp 960 torque_gen_old(:) = 0.0_wp 961 ENDIF 962 963 yawdir(:) = 0.0_wp 964 wdir_l(:) = 0.0_wp 965 wdir(:) = 0.0_wp 966 u_inflow(:) = 0.0_wp 967 u_inflow_l(:) = 0.0_wp 968 yaw_angle_l(:) = 0.0_wp 969 970 ! 971 !-- Allocate 1D arrays (dimension = number of rotor segments): 972 alpha_attack(:) = 0.0_wp 973 chord(:) = 0.0_wp 974 phi_rel(:) = 0.0_wp 975 thrust_seg(:) = 0.0_wp 976 torque_seg_y(:) = 0.0_wp 977 torque_seg_z(:) = 0.0_wp 978 turb_cd(:) = 0.0_wp 979 turb_cl(:) = 0.0_wp 980 vrel(:) = 0.0_wp 981 vtheta(:) = 0.0_wp 982 983 rbx(:,:) = 0.0_wp 984 rby(:,:) = 0.0_wp 985 rbz(:,:) = 0.0_wp 986 thrust_ring(:,:) = 0.0_wp 987 torque_ring_y(:,:) = 0.0_wp 988 torque_ring_z(:,:) = 0.0_wp 989 990 rotx(:,:) = 0.0_wp 991 roty(:,:) = 0.0_wp 992 rotz(:,:) = 0.0_wp 993 994 nac_cd_surf(:,:,:) = 0.0_wp 995 rot_tend_x(:,:,:) = 0.0_wp 996 rot_tend_y(:,:,:) = 0.0_wp 997 rot_tend_z(:,:,:) = 0.0_wp 998 thrust(:,:,:) = 0.0_wp 999 torque_y(:,:,:) = 0.0_wp 1000 torque_z(:,:,:) = 0.0_wp 1001 tow_cd_surf(:,:,:) = 0.0_wp 1002 1003 u_int(:,:,:) = 0.0_wp 1004 u_int_1_l(:,:,:) = 0.0_wp 1005 v_int(:,:,:) = 0.0_wp 1006 v_int_1_l(:,:,:) = 0.0_wp 1007 w_int(:,:,:) = 0.0_wp 1008 w_int_1_l(:,:,:) = 0.0_wp 1009 1010 1011 END SUBROUTINE wtm_init_arrays 1012 1013 1014 !--------------------------------------------------------------------------------------------------! 1022 1015 ! Description: 1023 1016 ! ------------ 1024 1017 !> Initialization of the wind turbine model 1025 !------------------------------------------------------------------------------! 1026 SUBROUTINE wtm_init 1027 1028 1029 USE control_parameters, & 1030 ONLY: dz_stretch_level_start 1031 1032 USE exchange_horiz_mod, & 1033 ONLY: exchange_horiz 1034 1035 IMPLICIT NONE 1036 1037 1038 1039 INTEGER(iwp) :: i !< running index 1040 INTEGER(iwp) :: j !< running index 1041 INTEGER(iwp) :: k !< running index 1042 1043 1044 ! 1045 !-- Help variables for the smearing function 1046 REAL(wp) :: eps_kernel !< 1047 1048 ! 1049 !-- Help variables for calculation of the tower drag 1050 INTEGER(iwp) :: tower_n !< 1051 INTEGER(iwp) :: tower_s !< 1052 1053 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: circle_points !< 1054 1055 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacb !< 1056 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacl !< 1057 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacr !< 1058 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nact !< 1059 1060 IF ( debug_output ) CALL debug_message( 'wtm_init', 'start' ) 1061 1062 ALLOCATE( index_nacb(1:n_turbines) ) 1063 ALLOCATE( index_nacl(1:n_turbines) ) 1064 ALLOCATE( index_nacr(1:n_turbines) ) 1065 ALLOCATE( index_nact(1:n_turbines) ) 1066 1067 ! 1068 !------------------------------------------------------------------------------! 1069 !-- Calculation of parameters for the regularization kernel 1070 !-- (smearing of the forces) 1071 !------------------------------------------------------------------------------! 1072 ! 1073 !-- In the following, some of the required parameters for the smearing will 1074 !-- be calculated: 1075 1076 !-- The kernel is set equal to twice the grid spacing which has turned out to 1077 !-- be a reasonable value (see e.g. Troldborg et al. (2013), Wind Energy, 1078 !-- DOI: 10.1002/we.1608): 1079 eps_kernel = 2.0_wp * dx 1080 ! 1081 !-- The zero point (eps_min) of the polynomial function must be the following 1082 !-- if the integral of the polynomial function (for values < eps_min) shall 1083 !-- be equal to the integral of the Gaussian function used before: 1084 eps_min = ( 105.0_wp / 32.0_wp )**( 1.0_wp / 3.0_wp ) * & 1085 pi**( 1.0_wp / 6.0_wp ) * eps_kernel 1086 ! 1087 !-- Stretching (non-uniform grid spacing) is not considered in the wind 1088 !-- turbine model. Therefore, vertical stretching has to be applied above 1089 !-- the area where the wtm is active. ABS (...) is required because the 1090 !-- default value of dz_stretch_level_start is -9999999.9_wp (negative). 1091 IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL(hub_z(1:n_turbines)) + & 1092 MAXVAL(rotor_radius(1:n_turbines)) + & 1093 eps_min) THEN 1094 WRITE( message_string, * ) 'The lowest level where vertical ', & 1095 'stretching is applied &have to be ', & 1096 'greater than ',MAXVAL(hub_z(1:n_turbines)) & 1097 + MAXVAL(rotor_radius(1:n_turbines)) + eps_min 1098 CALL message( 'wtm_init', 'PA0484', 1, 2, 0, 6, 0 ) 1099 ENDIF 1100 ! 1101 !-- Square of eps_min: 1102 eps_min2 = eps_min**2 1103 ! 1104 !-- Parameters in the polynomial function: 1105 pol_a = 1.0_wp / eps_min**4 1106 pol_b = 2.0_wp / eps_min**2 1107 ! 1108 !-- Normalization factor which is the inverse of the integral of the smearing 1109 !-- function: 1110 eps_factor = 105.0_wp / ( 32.0_wp * pi * eps_min**3 ) 1111 1112 !-- Change tilt angle to rad: 1113 tilt_angle = tilt_angle * pi / 180.0_wp 1114 1115 ! 1116 !-- Change yaw angle to rad: 1117 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1118 yaw_angle(:) = yaw_angle(:) * pi / 180.0_wp 1018 !--------------------------------------------------------------------------------------------------! 1019 SUBROUTINE wtm_init 1020 1021 1022 USE control_parameters, & 1023 ONLY: dz_stretch_level_start 1024 1025 USE exchange_horiz_mod, & 1026 ONLY: exchange_horiz 1027 1028 IMPLICIT NONE 1029 1030 1031 1032 INTEGER(iwp) :: i !< running index 1033 INTEGER(iwp) :: j !< running index 1034 INTEGER(iwp) :: k !< running index 1035 1036 1037 ! 1038 !-- Help variables for the smearing function: 1039 REAL(wp) :: eps_kernel !< 1040 1041 ! 1042 !-- Help variables for calculation of the tower drag: 1043 INTEGER(iwp) :: tower_n !< 1044 INTEGER(iwp) :: tower_s !< 1045 1046 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacb !< 1047 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacl !< 1048 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nacr !< 1049 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: index_nact !< 1050 1051 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: circle_points !< 1052 1053 1054 IF ( debug_output ) CALL debug_message( 'wtm_init', 'start' ) 1055 1056 ALLOCATE( index_nacb(1:n_turbines) ) 1057 ALLOCATE( index_nacl(1:n_turbines) ) 1058 ALLOCATE( index_nacr(1:n_turbines) ) 1059 ALLOCATE( index_nact(1:n_turbines) ) 1060 1061 ! 1062 !--------------------------------------------------------------------------------------------------! 1063 !-- Calculation of parameters for the regularization kernel (smearing of the forces) 1064 !--------------------------------------------------------------------------------------------------! 1065 ! 1066 !-- In the following, some of the required parameters for the smearing will be calculated: 1067 1068 !-- The kernel is set equal to twice the grid spacing which has turned out to be a reasonable 1069 !-- value (see e.g. Troldborg et al. (2013), Wind Energy, DOI: 10.1002/we.1608): 1070 eps_kernel = 2.0_wp * dx 1071 ! 1072 !-- The zero point (eps_min) of the polynomial function must be the following if the integral of 1073 !-- the polynomial function (for values < eps_min) shall be equal to the integral of the Gaussian 1074 !-- function used before: 1075 eps_min = ( 105.0_wp / 32.0_wp )**( 1.0_wp / 3.0_wp ) * & 1076 pi**( 1.0_wp / 6.0_wp ) * eps_kernel 1077 ! 1078 !-- Stretching (non-uniform grid spacing) is not considered in the wind turbine model. 1079 !-- Therefore, vertical stretching has to be applied above the area where the wtm is active. 1080 !-- ABS (...) is required because the default value of dz_stretch_level_start is -9999999.9_wp 1081 !-- (negative): 1082 IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL( hub_z(1:n_turbines) ) + & 1083 MAXVAL( rotor_radius(1:n_turbines) ) + & 1084 eps_min) THEN 1085 WRITE( message_string, * ) 'The lowest level where vertical stretching is applied has ' // & 1086 'to be greater than ',MAXVAL( hub_z(1:n_turbines) ) & 1087 + MAXVAL( rotor_radius(1:n_turbines) ) + eps_min 1088 CALL message( 'wtm_init', 'PA0484', 1, 2, 0, 6, 0 ) 1089 ENDIF 1090 ! 1091 !-- Square of eps_min: 1092 eps_min2 = eps_min**2 1093 ! 1094 !-- Parameters in the polynomial function: 1095 pol_a = 1.0_wp / eps_min**4 1096 pol_b = 2.0_wp / eps_min**2 1097 ! 1098 !-- Normalization factor which is the inverse of the integral of the smearing function: 1099 eps_factor = 105.0_wp / ( 32.0_wp * pi * eps_min**3 ) 1100 1101 !-- Change tilt angle to rad: 1102 tilt_angle = tilt_angle * pi / 180.0_wp 1103 1104 ! 1105 !-- Change yaw angle to rad: 1106 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1107 yaw_angle(:) = yaw_angle(:) * pi / 180.0_wp 1108 ENDIF 1109 1110 1111 DO inot = 1, n_turbines 1112 ! 1113 !-- Rotate the rotor coordinates in case yaw and tilt are defined: 1114 CALL wtm_rotate_rotor( inot ) 1115 1116 ! 1117 !-- Determine the indices of the hub height: 1118 i_hub(inot) = INT( hub_x(inot) / dx ) 1119 j_hub(inot) = INT( ( hub_y(inot) + 0.5_wp * dy ) / dy ) 1120 k_hub(inot) = INT( ( hub_z(inot) + 0.5_wp * dz(1) ) / dz(1) ) 1121 1122 ! 1123 !-- Determining the area to which the smearing of the forces is applied. 1124 !-- As smearing now is effectively applied only for distances smaller than eps_min, the 1125 !-- smearing area can be further limited and regarded as a function of eps_min: 1126 i_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dx ) 1127 j_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dy ) 1128 k_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dz(1) ) 1129 1130 ENDDO 1131 1132 ! 1133 !-- Call the wtm_init_speed_control subroutine and calculate the local rotor_speed for the 1134 !-- respective processor: 1135 IF ( speed_control) THEN 1136 1137 CALL wtm_init_speed_control 1138 1139 IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN 1140 1141 DO inot = 1, n_turbines 1142 1143 IF ( nxl > i_hub(inot) ) THEN 1144 torque_gen(inot) = 0.0_wp 1145 generator_speed_f(inot) = 0.0_wp 1146 rotor_speed_l(inot) = 0.0_wp 1147 ENDIF 1148 1149 IF ( nxr < i_hub(inot) ) THEN 1150 torque_gen(inot) = 0.0_wp 1151 generator_speed_f(inot) = 0.0_wp 1152 rotor_speed_l(inot) = 0.0_wp 1153 ENDIF 1154 1155 IF ( nys > j_hub(inot) ) THEN 1156 torque_gen(inot) = 0.0_wp 1157 generator_speed_f(inot) = 0.0_wp 1158 rotor_speed_l(inot) = 0.0_wp 1159 ENDIF 1160 1161 IF ( nyn < j_hub(inot) ) THEN 1162 torque_gen(inot) = 0.0_wp 1163 generator_speed_f(inot) = 0.0_wp 1164 rotor_speed_l(inot) = 0.0_wp 1165 ENDIF 1166 1167 IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) THEN 1168 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN 1169 1170 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 1171 1172 ENDIF 1173 ENDIF 1174 1175 END DO 1176 1119 1177 ENDIF 1120 1178 1121 1122 DO inot = 1, n_turbines 1123 ! 1124 !-- Rotate the rotor coordinates in case yaw and tilt are defined 1125 CALL wtm_rotate_rotor( inot ) 1126 1127 ! 1128 !-- Determine the indices of the hub height 1129 i_hub(inot) = INT( hub_x(inot) / dx ) 1130 j_hub(inot) = INT( ( hub_y(inot) + 0.5_wp * dy ) / dy ) 1131 k_hub(inot) = INT( ( hub_z(inot) + 0.5_wp * dz(1) ) / dz(1) ) 1132 1133 ! 1134 !-- Determining the area to which the smearing of the forces is applied. 1135 !-- As smearing now is effectively applied only for distances smaller than 1136 !-- eps_min, the smearing area can be further limited and regarded as a 1137 !-- function of eps_min: 1138 i_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dx ) 1139 j_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dy ) 1140 k_smear(inot) = CEILING( ( rotor_radius(inot) + eps_min ) / dz(1) ) 1141 1142 ENDDO 1143 1144 ! 1145 !-- Call the wtm_init_speed_control subroutine and calculate the local 1146 !-- rotor_speed for the respective processor. 1147 IF ( speed_control) THEN 1148 1149 CALL wtm_init_speed_control 1150 1151 IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN 1152 1153 DO inot = 1, n_turbines 1154 1155 IF ( nxl > i_hub(inot) ) THEN 1156 torque_gen(inot) = 0.0_wp 1157 generator_speed_f(inot) = 0.0_wp 1158 rotor_speed_l(inot) = 0.0_wp 1159 ENDIF 1160 1161 IF ( nxr < i_hub(inot) ) THEN 1162 torque_gen(inot) = 0.0_wp 1163 generator_speed_f(inot) = 0.0_wp 1164 rotor_speed_l(inot) = 0.0_wp 1165 ENDIF 1166 1167 IF ( nys > j_hub(inot) ) THEN 1168 torque_gen(inot) = 0.0_wp 1169 generator_speed_f(inot) = 0.0_wp 1170 rotor_speed_l(inot) = 0.0_wp 1171 ENDIF 1172 1173 IF ( nyn < j_hub(inot) ) THEN 1174 torque_gen(inot) = 0.0_wp 1175 generator_speed_f(inot) = 0.0_wp 1176 rotor_speed_l(inot) = 0.0_wp 1177 ENDIF 1178 1179 IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) THEN 1180 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN 1181 1182 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 1183 1184 ENDIF 1185 ENDIF 1186 1187 END DO 1188 1189 ENDIF 1190 1191 ENDIF 1192 1193 ! 1194 !------------------------------------------------------------------------------! 1195 !-- Determine the area within each grid cell that overlaps with the area 1196 !-- of the nacelle and the tower (needed for calculation of the forces) 1197 !------------------------------------------------------------------------------! 1198 ! 1199 !-- Note: so far this is only a 2D version, in that the mean flow is 1200 !-- perpendicular to the rotor area. 1201 1202 ! 1203 !-- Allocation of the array containing information on the intersection points 1204 !-- between rotor disk and the numerical grid: 1205 upper_end = ( ny + 1 ) * 10000 1206 1207 ALLOCATE( circle_points(1:2,1:upper_end) ) 1208 1209 circle_points(:,:) = 0.0_wp 1210 1211 1212 DO inot = 1, n_turbines ! loop over number of turbines 1213 ! 1214 !-- Determine the grid index (u-grid) that corresponds to the location of 1215 !-- the rotor center (reduces the amount of calculations in the case that 1216 !-- the mean flow is perpendicular to the rotor area): 1217 i = i_hub(inot) 1218 1219 ! 1220 !-- Determine the left and the right edge of the nacelle (corresponding 1221 !-- grid point indices): 1222 index_nacl(inot) = INT( ( hub_y(inot) - nacelle_radius(inot) + 0.5_wp * dy ) / dy ) 1223 index_nacr(inot) = INT( ( hub_y(inot) + nacelle_radius(inot) + 0.5_wp * dy ) / dy ) 1224 ! 1225 !-- Determine the bottom and the top edge of the nacelle (corresponding 1226 !-- grid point indices).The grid point index has to be increased by 1, as 1227 !-- the first level for the u-component (index 0) is situated below the 1228 !-- surface. All points between z=0 and z=dz/s would already be contained 1229 !-- in grid box 1. 1230 index_nacb(inot) = INT( ( hub_z(inot) - nacelle_radius(inot) ) / dz(1) ) + 1 1231 index_nact(inot) = INT( ( hub_z(inot) + nacelle_radius(inot) ) / dz(1) ) + 1 1232 1233 ! 1234 !-- Determine the indices of the grid boxes containing the left and 1235 !-- the right boundaries of the tower: 1236 tower_n = ( hub_y(inot) + 0.5_wp * tower_diameter(inot) - 0.5_wp * dy ) / dy 1237 tower_s = ( hub_y(inot) - 0.5_wp * tower_diameter(inot) - 0.5_wp * dy ) / dy 1238 1239 ! 1240 !-- Determine the fraction of the grid box area overlapping with the tower 1241 !-- area and multiply it with the drag of the tower: 1242 IF ( ( nxlg <= i ) .AND. ( nxrg >= i ) ) THEN 1243 1244 DO j = nys, nyn 1245 ! 1246 !-- Loop from south to north boundary of tower 1247 IF ( ( j >= tower_s ) .AND. ( j <= tower_n ) ) THEN 1248 1249 DO k = nzb, nzt 1250 1251 IF ( k == k_hub(inot) ) THEN 1252 IF ( tower_n - tower_s >= 1 ) THEN 1253 ! 1254 !-- leftmost and rightmost grid box: 1255 IF ( j == tower_s ) THEN 1256 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1257 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1258 ( ( tower_s + 1.0_wp + 0.5_wp ) * dy - & 1259 ( hub_y(inot) - 0.5_wp * tower_diameter(inot) ) ) * & ! extension in y-direction 1260 tower_cd(inot) 1261 ELSEIF ( j == tower_n ) THEN 1262 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1263 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1264 ( ( hub_y(inot) + 0.5_wp * tower_diameter(inot) ) - & 1265 ( tower_n + 0.5_wp ) * dy ) * & ! extension in y-direction 1266 tower_cd(inot) 1267 ! 1268 !-- grid boxes inbetween 1269 !-- (where tow_cd_surf = grid box area): 1270 ELSE 1271 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1272 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& 1273 dy * tower_cd(inot) 1274 ENDIF 1275 ! 1276 !-- tower lies completely within one grid box: 1179 ENDIF 1180 1181 ! 1182 !--------------------------------------------------------------------------------------------------! 1183 !-- Determine the area within each grid cell that overlaps with the area of the nacelle and the 1184 !-- tower (needed for calculation of the forces) 1185 !--------------------------------------------------------------------------------------------------! 1186 ! 1187 !-- Note: so far this is only a 2D version, in that the mean flow is perpendicular to the rotor 1188 !-- area. 1189 1190 ! 1191 !-- Allocation of the array containing information on the intersection points between rotor disk 1192 !-- and the numerical grid: 1193 upper_end = ( ny + 1 ) * 10000 1194 1195 ALLOCATE( circle_points(1:2,1:upper_end) ) 1196 1197 circle_points(:,:) = 0.0_wp 1198 1199 1200 DO inot = 1, n_turbines ! loop over number of turbines 1201 ! 1202 !-- Determine the grid index (u-grid) that corresponds to the location of the rotor center 1203 !-- (reduces the amount of calculations in the case that the mean flow is perpendicular to the 1204 !-- rotor area): 1205 i = i_hub(inot) 1206 1207 ! 1208 !-- Determine the left and the right edge of the nacelle (corresponding grid point indices): 1209 index_nacl(inot) = INT( ( hub_y(inot) - nacelle_radius(inot) + 0.5_wp * dy ) / dy ) 1210 index_nacr(inot) = INT( ( hub_y(inot) + nacelle_radius(inot) + 0.5_wp * dy ) / dy ) 1211 ! 1212 !-- Determine the bottom and the top edge of the nacelle (corresponding grid point indices). 1213 !-- The grid point index has to be increased by 1, as the first level for the u-component 1214 !-- (index 0) is situated below the surface. All points between z=0 and z=dz/s would already 1215 !-- be contained in grid box 1: 1216 index_nacb(inot) = INT( ( hub_z(inot) - nacelle_radius(inot) ) / dz(1) ) + 1 1217 index_nact(inot) = INT( ( hub_z(inot) + nacelle_radius(inot) ) / dz(1) ) + 1 1218 1219 ! 1220 !-- Determine the indices of the grid boxes containing the left and the right boundaries of 1221 !-- the tower: 1222 tower_n = ( hub_y(inot) + 0.5_wp * tower_diameter(inot) - 0.5_wp * dy ) / dy 1223 tower_s = ( hub_y(inot) - 0.5_wp * tower_diameter(inot) - 0.5_wp * dy ) / dy 1224 1225 ! 1226 !-- Determine the fraction of the grid box area overlapping with the tower area and multiply 1227 !-- it with the drag of the tower: 1228 IF ( ( nxlg <= i ) .AND. ( nxrg >= i ) ) THEN 1229 1230 DO j = nys, nyn 1231 ! 1232 !-- Loop from south to north boundary of tower: 1233 IF ( ( j >= tower_s ) .AND. ( j <= tower_n ) ) THEN 1234 1235 DO k = nzb, nzt 1236 1237 IF ( k == k_hub(inot) ) THEN 1238 IF ( tower_n - tower_s >= 1 ) THEN 1239 ! 1240 !-- Leftmost and rightmost grid box: 1241 IF ( j == tower_s ) THEN 1242 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1243 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) ) * & ! extension in z-direction 1244 ( ( tower_s + 1.0_wp + 0.5_wp ) * dy - & 1245 ( hub_y(inot) - 0.5_wp * tower_diameter(inot) ) ) * & ! extension in y-direction 1246 tower_cd(inot) 1247 1248 ELSEIF ( j == tower_n ) THEN 1249 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1250 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) ) * & ! extension in z-direction 1251 ( ( hub_y(inot) + 0.5_wp * tower_diameter(inot) ) - & 1252 ( tower_n + 0.5_wp ) * dy ) * & ! extension in y-direction 1253 tower_cd(inot) 1254 ! 1255 !-- Grid boxes inbetween (where tow_cd_surf = grid box area): 1277 1256 ELSE 1278 tow_cd_surf(k,j,i) = ( hub_z(inot) - ( k_hub(inot) *&1279 dz(1) - 0.5_wp * dz(1) ) ) *&1280 tower_diameter(inot) *tower_cd(inot)1257 tow_cd_surf(k,j,i) = ( hub_z(inot) - & 1258 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) ) * dy * & 1259 tower_cd(inot) 1281 1260 ENDIF 1282 1261 ! 1283 !-- In case that k is smaller than k_hub the following actions 1284 !-- are carried out: 1285 ELSEIF ( k < k_hub(inot) ) THEN 1286 1287 IF ( ( tower_n - tower_s ) >= 1 ) THEN 1288 ! 1289 !-- leftmost and rightmost grid box: 1290 IF ( j == tower_s ) THEN 1291 tow_cd_surf(k,j,i) = dz(1) * ( & 1292 ( tower_s + 1 + 0.5_wp ) * dy - & 1293 ( hub_y(inot) - 0.5_wp * tower_diameter(inot) ) & 1294 ) * tower_cd(inot) 1295 ELSEIF ( j == tower_n ) THEN 1296 tow_cd_surf(k,j,i) = dz(1) * ( & 1297 ( hub_y(inot) + 0.5_wp * tower_diameter(inot) ) - & 1298 ( tower_n + 0.5_wp ) * dy & 1299 ) * tower_cd(inot) 1300 ! 1301 !-- grid boxes inbetween 1302 !-- (where tow_cd_surf = grid box area): 1303 ELSE 1304 tow_cd_surf(k,j,i) = dz(1) * dy * & 1305 tower_cd(inot) 1306 ENDIF 1307 ! 1308 !-- tower lies completely within one grid box: 1262 !-- Tower lies completely within one grid box: 1263 ELSE 1264 tow_cd_surf(k,j,i) = ( hub_z(inot) - ( k_hub(inot) * & 1265 dz(1) - 0.5_wp * dz(1) ) ) * & 1266 tower_diameter(inot) * tower_cd(inot) 1267 ENDIF 1268 ! 1269 !-- In case that k is smaller than k_hub the following actions are carried out: 1270 ELSEIF ( k < k_hub(inot) ) THEN 1271 1272 IF ( ( tower_n - tower_s ) >= 1 ) THEN 1273 ! 1274 !-- Leftmost and rightmost grid box: 1275 IF ( j == tower_s ) THEN 1276 tow_cd_surf(k,j,i) = dz(1) * ( ( tower_s + 1 + 0.5_wp ) * dy - & 1277 ( hub_y(inot) - 0.5_wp * tower_diameter(inot) ) ) & 1278 * tower_cd(inot) 1279 1280 ELSEIF ( j == tower_n ) THEN 1281 tow_cd_surf(k,j,i) = dz(1) * ( ( hub_y(inot) + 0.5_wp * & 1282 tower_diameter(inot) ) - ( tower_n + 0.5_wp ) & 1283 * dy ) * tower_cd(inot) 1284 ! 1285 !-- Grid boxes inbetween (where tow_cd_surf = grid box area): 1309 1286 ELSE 1310 tow_cd_surf(k,j,i) = dz(1) * tower_diameter(inot) * & 1311 tower_cd(inot) 1312 ENDIF ! end if larger than grid box 1313 1314 ENDIF ! end if k == k_hub 1315 1316 ENDDO ! end loop over k 1317 1318 ENDIF ! end if inside north and south boundary of tower 1319 1320 ENDDO ! end loop over j 1321 1322 ENDIF ! end if hub inside domain + ghostpoints 1323 1324 1325 CALL exchange_horiz( tow_cd_surf, nbgp ) 1326 1327 ! 1328 !-- Calculation of the nacelle area 1329 !-- CAUTION: Currently disabled due to segmentation faults on the FLOW HPC 1330 !-- cluster (Oldenburg) 1287 tow_cd_surf(k,j,i) = dz(1) * dy * tower_cd(inot) 1288 ENDIF 1289 ! 1290 !-- Tower lies completely within one grid box: 1291 ELSE 1292 tow_cd_surf(k,j,i) = dz(1) * tower_diameter(inot) * tower_cd(inot) 1293 ENDIF ! end if larger than grid box 1294 1295 ENDIF ! end if k == k_hub 1296 1297 ENDDO ! end loop over k 1298 1299 ENDIF ! end if inside north and south boundary of tower 1300 1301 ENDDO ! end loop over j 1302 1303 ENDIF ! end if hub inside domain + ghostpoints 1304 1305 1306 CALL exchange_horiz( tow_cd_surf, nbgp ) 1307 1308 ! 1309 !-- Calculation of the nacelle area 1310 !-- CAUTION: Currently disabled due to segmentation faults on the FLOW HPC cluster (Oldenburg) 1331 1311 !! 1332 !!-- Tabulate the points on the circle that are required in the following for 1333 !!-- the calculation of the Riemann integral (node points; they are called 1334 !!-- circle_points in the following): 1335 ! 1336 ! dy_int = dy / 10000.0_wp 1337 ! 1338 ! DO i_ip = 1, upper_end 1339 ! yvalue = dy_int * ( i_ip - 0.5_wp ) + 0.5_wp * dy !<--- segmentation fault 1340 ! sqrt_arg = nacelle_radius(inot)**2 - ( yvalue - hub_y(inot) )**2 !<--- segmentation fault 1341 ! IF ( sqrt_arg >= 0.0_wp ) THEN 1312 !!-- Tabulate the points on the circle that are required in the following for the calculation 1313 !!-- of the Riemann integral (node points; they are called circle_points in the following): 1314 ! 1315 ! dy_int = dy / 10000.0_wp 1316 ! 1317 ! DO i_ip = 1, upper_end 1318 ! yvalue = dy_int * ( i_ip - 0.5_wp ) + 0.5_wp * dy !<--- segmentation fault 1319 ! sqrt_arg = nacelle_radius(inot)**2 - ( yvalue - hub_y(inot) )**2 !<--- segmentation fault 1320 ! IF ( sqrt_arg >= 0.0_wp ) THEN 1342 1321 !! 1343 !!-- 1344 ! 1322 !!-- bottom intersection point 1323 ! circle_points(1,i_ip) = hub_z(inot) - SQRT( sqrt_arg ) 1345 1324 !! 1346 !!-- 1347 ! circle_points(2,i_ip) = hub_z(inot) + SQRT( sqrt_arg )!<--- segmentation fault1348 ! 1349 ! circle_points(:,i_ip) = -111111!<--- segmentation fault1350 ! 1351 ! 1352 ! 1353 ! 1354 ! 1325 !!-- top intersection point 1326 ! circle_points(2,i_ip) = hub_z(inot) + SQRT( sqrt_arg ) !<--- segmentation fault 1327 ! ELSE 1328 ! circle_points(:,i_ip) = -111111 !<--- segmentation fault 1329 ! ENDIF 1330 ! ENDDO 1331 ! 1332 ! 1333 ! DO j = nys, nyn 1355 1334 !! 1356 !!-- In case that the grid box is located completely outside the nacelle 1357 !!-- (y) it can automatically be stated that there is no overlap between 1358 !!-- the grid box and the nacelle and consequently we can set 1359 !!-- nac_cd_surf(:,j,i) = 0.0: 1360 ! IF ( ( j >= index_nacl(inot) ) .AND. ( j <= index_nacr(inot) ) ) THEN 1361 ! DO k = nzb+1, nzt 1335 !!-- In case that the grid box is located completely outside the nacelle (y) it can 1336 !!-- automatically be stated that there is no overlap between the grid box and the nacelle 1337 !!-- and consequently we can set nac_cd_surf(:,j,i) = 0.0: 1338 ! IF ( ( j >= index_nacl(inot) ) .AND. ( j <= index_nacr(inot) ) ) THEN 1339 ! DO k = nzb+1, nzt 1362 1340 !! 1363 !!-- In case that the grid box is located completely outside the 1364 !!-- nacelle (z) it can automatically be stated that there is no 1365 !!-- overlap between the grid box and the nacelle and consequently 1366 !!-- we can set nac_cd_surf(k,j,i) = 0.0: 1367 ! IF ( ( k >= index_nacb(inot) ) .OR. & 1368 ! ( k <= index_nact(inot) ) ) THEN 1341 !!-- In case that the grid box is located completely outside the nacelle (z) it can 1342 !!-- automatically be stated that there is no overlap between the grid box and the 1343 !!-- nacelle and consequently we can set nac_cd_surf(k,j,i) = 0.0: 1344 ! IF ( ( k >= index_nacb(inot) ) .OR. & 1345 ! ( k <= index_nact(inot) ) ) THEN 1369 1346 !! 1370 !!-- For all other cases Riemann integrals are calculated. 1371 !!-- Here, the points on the circle that have been determined 1372 !!-- above are used in order to calculate the overlap between the 1373 !!-- gridbox and the nacelle area (area approached by 10000 1374 !!-- rectangulars dz_int * dy_int): 1375 ! DO i_ipg = 1, 10000 1376 ! dz_int = dz 1377 ! i_ip = j * 10000 + i_ipg 1347 !!-- For all other cases Riemann integrals are calculated. Here, the points on the 1348 !!-- circle that have been determined above are used in order to calculate the 1349 !!-- overlap between the gridbox and the nacelle area (area approached by 10000 1350 !!-- rectangulars dz_int * dy_int): 1351 ! DO i_ipg = 1, 10000 1352 ! dz_int = dz 1353 ! i_ip = j * 10000 + i_ipg 1378 1354 !! 1379 !!-- Determine the vertical extension dz_int of the circle1380 !!-- within the currentgrid box:1381 ! 1382 ! 1383 ! 1384 ! 1385 ! 1386 ! 1387 ! 1388 ! 1389 ! 1390 ! 1391 ! 1392 ! 1393 ! 1394 ! 1395 ! dz_int = 0.0_wp1396 ! ENDIF1397 ! 1398 ! 1399 ! 1400 ! ENDIF1401 ! 1402 ! 1403 ! 1404 ! 1405 ! 1406 ! 1407 ! 1408 ! CALL exchange_horiz( nac_cd_surf, nbgp )!<--- segmentation fault1409 1410 ENDDO! end of loop over turbines1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1355 !!-- Determine the vertical extension dz_int of the circle within the current 1356 !!-- grid box: 1357 ! IF ( ( circle_points(2,i_ip) < zw(k) ) .AND. & !<--- segmentation fault 1358 ! ( circle_points(2,i_ip) >= zw(k-1) ) ) THEN 1359 ! dz_int = dz_int - & !<--- segmentation fault 1360 ! ( zw(k) - circle_points(2,i_ip) ) 1361 ! ENDIF 1362 ! IF ( ( circle_points(1,i_ip) <= zw(k) ) .AND. & !<--- segmentation fault 1363 ! ( circle_points(1,i_ip) > zw(k-1) ) ) THEN 1364 ! dz_int = dz_int - & 1365 ! ( circle_points(1,i_ip) - zw(k-1) ) 1366 ! ENDIF 1367 ! IF ( zw(k-1) > circle_points(2,i_ip) ) THEN 1368 ! dz_int = 0.0_wp 1369 ! ENDIF 1370 ! IF ( zw(k) < circle_points(1,i_ip) ) THEN 1371 ! dz_int = 0.0_wp 1372 ! ENDIF 1373 ! IF ( ( nxlg <= i ) .AND. ( nxrg >= i ) ) THEN 1374 ! nac_cd_surf(k,j,i) = nac_cd_surf(k,j,i) + & !<--- segmentation fault 1375 ! dy_int * dz_int * nacelle_cd(inot) 1376 ! ENDIF 1377 ! ENDDO 1378 ! ENDIF 1379 ! ENDDO 1380 ! ENDIF 1381 ! 1382 ! ENDDO 1383 ! 1384 ! CALL exchange_horiz( nac_cd_surf, nbgp ) !<--- segmentation fault 1385 1386 ENDDO ! end of loop over turbines 1387 1388 tow_cd_surf = tow_cd_surf / ( dx * dy * dz(1) ) ! Normalize tower drag 1389 nac_cd_surf = nac_cd_surf / ( dx * dy * dz(1) ) ! Normalize nacelle drag 1390 1391 CALL wtm_read_blade_tables 1392 1393 IF ( debug_output ) CALL debug_message( 'wtm_init', 'end' ) 1394 1395 END SUBROUTINE wtm_init 1396 1397 1398 1423 1399 SUBROUTINE wtm_init_output 1424 1425 1426 ! INTEGER(iwp) :: ntimesteps 1427 ! INTEGER(iwp) :: ntimesteps_max = 80000 1428 INTEGER(iwp) :: return_value 1429 1430 INTEGER(iwp) :: n !< running index 1431 1432 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ndim !< dummy to write dimension 1433 1434 1435 ! 1436 !-- Create NetCDF output file 1400 1401 1402 ! INTEGER(iwp) :: ntimesteps !< number of timesteps defined in NetCDF output file 1403 ! INTEGER(iwp) :: ntimesteps_max = 80000 !< number of maximum timesteps defined in NetCDF output file 1404 INTEGER(iwp) :: return_value !< returned status value of called function 1405 1406 INTEGER(iwp) :: n !< running index 1407 1408 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ndim !< dummy to write dimension 1409 1410 1411 ! 1412 !-- Create NetCDF output file: 1437 1413 #if defined( __netcdf4 ) 1438 1414 nc_filename = 'DATA_1D_TS_WTM_NETCDF' // TRIM( coupling_char ) 1439 1415 return_value = dom_def_file( nc_filename, 'netcdf4-serial' ) 1440 1416 #else 1441 message_string = 'Wind turbine model output requires NetCDF version 4. ' // &1417 message_string = 'Wind turbine model output requires NetCDF version 4. ' // & 1442 1418 'No output file will be created.' 1443 CALL message( 'wtm_init_output', 'PA0672', & 1444 0, 1, 0, 6, 0 ) 1419 CALL message( 'wtm_init_output', 'PA0672', 0, 1, 0, 6, 0 ) 1445 1420 #endif 1446 1421 IF ( myid == 0 ) THEN 1447 1422 ! 1448 !-- Define dimensions in output file 1423 !-- Define dimensions in output file: 1449 1424 ALLOCATE( ndim(1:n_turbines) ) 1450 1425 DO n = 1, n_turbines 1451 1426 ndim(n) = n 1452 1427 ENDDO 1453 return_value = dom_def_dim( nc_filename, &1454 dimension_name = 'turbine', &1455 output_type = 'int32', &1456 bounds = (/1_iwp, n_turbines/), &1428 return_value = dom_def_dim( nc_filename, & 1429 dimension_name = 'turbine', & 1430 output_type = 'int32', & 1431 bounds = (/1_iwp, n_turbines/), & 1457 1432 values_int32 = ndim ) 1458 1433 DEALLOCATE( ndim ) 1459 1434 1460 return_value = dom_def_dim( nc_filename, &1461 dimension_name = 'time', &1462 output_type = 'real32', &1463 bounds = (/1_iwp/), &1435 return_value = dom_def_dim( nc_filename, & 1436 dimension_name = 'time', & 1437 output_type = 'real32', & 1438 bounds = (/1_iwp/), & 1464 1439 values_realwp = (/0.0_wp/) ) 1465 1440 1466 1441 variable_name = 'x' 1467 return_value = dom_def_var( nc_filename, &1468 variable_name = variable_name, &1469 dimension_names = (/'turbine'/), &1442 return_value = dom_def_var( nc_filename, & 1443 variable_name = variable_name, & 1444 dimension_names = (/'turbine'/), & 1470 1445 output_type = 'real32' ) 1471 1446 1472 1447 variable_name = 'y' 1473 return_value = dom_def_var( nc_filename, &1474 variable_name = variable_name, &1475 dimension_names = (/'turbine'/), &1448 return_value = dom_def_var( nc_filename, & 1449 variable_name = variable_name, & 1450 dimension_names = (/'turbine'/), & 1476 1451 output_type = 'real32' ) 1477 1452 1478 1453 variable_name = 'z' 1479 return_value = dom_def_var( nc_filename, &1480 variable_name = variable_name, &1481 dimension_names = (/'turbine'/), &1454 return_value = dom_def_var( nc_filename, & 1455 variable_name = variable_name, & 1456 dimension_names = (/'turbine'/), & 1482 1457 output_type = 'real32' ) 1483 1458 1484 1459 1485 1460 variable_name = 'rotor_diameter' 1486 return_value = dom_def_var( nc_filename, &1487 variable_name = variable_name, &1488 dimension_names = (/'turbine'/), &1461 return_value = dom_def_var( nc_filename, & 1462 variable_name = variable_name, & 1463 dimension_names = (/'turbine'/), & 1489 1464 output_type = 'real32' ) 1490 1465 1491 1466 variable_name = 'tower_diameter' 1492 return_value = dom_def_var( nc_filename, &1493 variable_name = variable_name, &1494 dimension_names = (/'turbine'/), &1467 return_value = dom_def_var( nc_filename, & 1468 variable_name = variable_name, & 1469 dimension_names = (/'turbine'/), & 1495 1470 output_type = 'real32' ) 1496 1471 1497 return_value = dom_def_att( nc_filename, &1498 variable_name = 'time', &1499 attribute_name = 'units', &1472 return_value = dom_def_att( nc_filename, & 1473 variable_name = 'time', & 1474 attribute_name = 'units', & 1500 1475 value = 'seconds since ' // origin_date_time ) 1501 1476 1502 return_value = dom_def_att( nc_filename, &1503 variable_name = 'x', &1504 attribute_name = 'units', &1477 return_value = dom_def_att( nc_filename, & 1478 variable_name = 'x', & 1479 attribute_name = 'units', & 1505 1480 value = 'm' ) 1506 1481 1507 return_value = dom_def_att( nc_filename, & 1508 variable_name = 'y', & 1509 attribute_name = 'units', & 1510 value = 'm' ) 1511 1512 return_value = dom_def_att( nc_filename, & 1513 variable_name = 'z', & 1514 attribute_name = 'units', & 1482 return_value = dom_def_att( nc_filename, & 1483 variable_name = 'y', & 1484 attribute_name = 'units', & 1515 1485 value = 'm' ) 1516 1486 1517 return_value = dom_def_att( nc_filename, &1518 variable_name = ' rotor_diameter',&1519 attribute_name = 'units', &1487 return_value = dom_def_att( nc_filename, & 1488 variable_name = 'z', & 1489 attribute_name = 'units', & 1520 1490 value = 'm' ) 1521 1491 1522 return_value = dom_def_att( nc_filename, &1523 variable_name = ' tower_diameter',&1524 attribute_name = 'units', &1492 return_value = dom_def_att( nc_filename, & 1493 variable_name = 'rotor_diameter', & 1494 attribute_name = 'units', & 1525 1495 value = 'm' ) 1526 1496 1527 return_value = dom_def_att( nc_filename, & 1528 variable_name = 'x', & 1529 attribute_name = 'long_name', & 1497 return_value = dom_def_att( nc_filename, & 1498 variable_name = 'tower_diameter', & 1499 attribute_name = 'units', & 1500 value = 'm' ) 1501 1502 return_value = dom_def_att( nc_filename, & 1503 variable_name = 'x', & 1504 attribute_name = 'long_name', & 1530 1505 value = 'x location of rotor center' ) 1531 1506 1532 return_value = dom_def_att( nc_filename, &1533 variable_name = 'y', &1534 attribute_name = 'long_name', &1535 value = 'y location of rotor center' ) 1536 1537 return_value = dom_def_att( nc_filename, &1538 variable_name = 'z', &1539 attribute_name = 'long_name', &1540 value = 'z location of rotor center' ) 1541 1542 return_value = dom_def_att( nc_filename, &1543 variable_name = 'turbine_name', &1544 attribute_name = 'long_name', &1507 return_value = dom_def_att( nc_filename, & 1508 variable_name = 'y', & 1509 attribute_name = 'long_name', & 1510 value = 'y location of rotor center' ) 1511 1512 return_value = dom_def_att( nc_filename, & 1513 variable_name = 'z', & 1514 attribute_name = 'long_name', & 1515 value = 'z location of rotor center' ) 1516 1517 return_value = dom_def_att( nc_filename, & 1518 variable_name = 'turbine_name', & 1519 attribute_name = 'long_name', & 1545 1520 value = 'turbine name') 1546 1521 1547 return_value = dom_def_att( nc_filename, &1548 variable_name = 'rotor_diameter', &1549 attribute_name = 'long_name', &1522 return_value = dom_def_att( nc_filename, & 1523 variable_name = 'rotor_diameter', & 1524 attribute_name = 'long_name', & 1550 1525 value = 'rotor diameter') 1551 1526 1552 return_value = dom_def_att( nc_filename, &1553 variable_name = 'tower_diameter', &1554 attribute_name = 'long_name', &1527 return_value = dom_def_att( nc_filename, & 1528 variable_name = 'tower_diameter', & 1529 attribute_name = 'long_name', & 1555 1530 value = 'tower diameter') 1556 1531 1557 return_value = dom_def_att( nc_filename, &1558 variable_name = 'time', &1559 attribute_name = 'standard_name', &1532 return_value = dom_def_att( nc_filename, & 1533 variable_name = 'time', & 1534 attribute_name = 'standard_name', & 1560 1535 value = 'time') 1561 1536 1562 return_value = dom_def_att( nc_filename, &1563 variable_name = 'time', &1564 attribute_name = 'axis', &1537 return_value = dom_def_att( nc_filename, & 1538 variable_name = 'time', & 1539 attribute_name = 'axis', & 1565 1540 value = 'T') 1566 1541 1567 return_value = dom_def_att( nc_filename, &1568 variable_name = 'x', &1569 attribute_name = 'axis', &1542 return_value = dom_def_att( nc_filename, & 1543 variable_name = 'x', & 1544 attribute_name = 'axis', & 1570 1545 value = 'X' ) 1571 1546 1572 return_value = dom_def_att( nc_filename, &1573 variable_name = 'y', &1574 attribute_name = 'axis', &1547 return_value = dom_def_att( nc_filename, & 1548 variable_name = 'y', & 1549 attribute_name = 'axis', & 1575 1550 value = 'Y' ) 1576 1551 1577 1552 1578 1553 variable_name = 'generator_power' 1579 return_value = dom_def_var( nc_filename, &1580 variable_name = variable_name, &1581 dimension_names = (/ 'turbine', 'time ' /), &1554 return_value = dom_def_var( nc_filename, & 1555 variable_name = variable_name, & 1556 dimension_names = (/ 'turbine', 'time ' /), & 1582 1557 output_type = 'real32' ) 1583 1558 1584 return_value = dom_def_att( nc_filename, &1585 variable_name = variable_name, &1586 attribute_name = 'units', &1587 value = 'W' ) 1559 return_value = dom_def_att( nc_filename, & 1560 variable_name = variable_name, & 1561 attribute_name = 'units', & 1562 value = 'W' ) 1588 1563 1589 1564 variable_name = 'generator_speed' 1590 return_value = dom_def_var( nc_filename, &1591 variable_name = variable_name, &1592 dimension_names = (/ 'turbine', 'time ' /), &1565 return_value = dom_def_var( nc_filename, & 1566 variable_name = variable_name, & 1567 dimension_names = (/ 'turbine', 'time ' /), & 1593 1568 output_type = 'real32' ) 1594 1569 1595 return_value = dom_def_att( nc_filename, &1596 variable_name = variable_name, &1597 attribute_name = 'units', &1570 return_value = dom_def_att( nc_filename, & 1571 variable_name = variable_name, & 1572 attribute_name = 'units', & 1598 1573 value = 'rad/s' ) 1599 1574 1600 1575 variable_name = 'generator_torque' 1601 return_value = dom_def_var( nc_filename, &1602 variable_name = variable_name, &1603 dimension_names = (/ 'turbine', 'time ' /), &1576 return_value = dom_def_var( nc_filename, & 1577 variable_name = variable_name, & 1578 dimension_names = (/ 'turbine', 'time ' /), & 1604 1579 output_type = 'real32' ) 1605 1606 return_value = dom_def_att( nc_filename, &1607 variable_name = variable_name, &1608 attribute_name = 'units', &1580 1581 return_value = dom_def_att( nc_filename, & 1582 variable_name = variable_name, & 1583 attribute_name = 'units', & 1609 1584 value = 'Nm' ) 1610 1585 1611 1586 variable_name = 'pitch_angle' 1612 return_value = dom_def_var( nc_filename, &1613 variable_name = variable_name, &1614 dimension_names = (/ 'turbine', 'time ' /), &1587 return_value = dom_def_var( nc_filename, & 1588 variable_name = variable_name, & 1589 dimension_names = (/ 'turbine', 'time ' /), & 1615 1590 output_type = 'real32' ) 1616 1591 1617 return_value = dom_def_att( nc_filename, &1618 variable_name = variable_name, &1619 attribute_name = 'units', &1592 return_value = dom_def_att( nc_filename, & 1593 variable_name = variable_name, & 1594 attribute_name = 'units', & 1620 1595 value = 'degrees' ) 1621 1596 1622 1597 variable_name = 'rotor_power' 1623 return_value = dom_def_var( nc_filename, &1624 variable_name = variable_name, &1625 dimension_names = (/ 'turbine', 'time ' /), &1598 return_value = dom_def_var( nc_filename, & 1599 variable_name = variable_name, & 1600 dimension_names = (/ 'turbine', 'time ' /), & 1626 1601 output_type = 'real32' ) 1627 1628 return_value = dom_def_att( nc_filename, &1629 variable_name = variable_name, &1630 attribute_name = 'units', &1631 value = 'W' ) 1602 1603 return_value = dom_def_att( nc_filename, & 1604 variable_name = variable_name, & 1605 attribute_name = 'units', & 1606 value = 'W' ) 1632 1607 1633 1608 variable_name = 'rotor_speed' 1634 return_value = dom_def_var( nc_filename, &1635 variable_name = variable_name, &1636 dimension_names = (/ 'turbine', 'time ' /), &1609 return_value = dom_def_var( nc_filename, & 1610 variable_name = variable_name, & 1611 dimension_names = (/ 'turbine', 'time ' /), & 1637 1612 output_type = 'real32' ) 1638 1639 return_value = dom_def_att( nc_filename, &1640 variable_name = variable_name, &1641 attribute_name = 'units', &1613 1614 return_value = dom_def_att( nc_filename, & 1615 variable_name = variable_name, & 1616 attribute_name = 'units', & 1642 1617 value = 'rad/s' ) 1643 1618 1644 1619 variable_name = 'rotor_thrust' 1645 return_value = dom_def_var( nc_filename, &1646 variable_name = variable_name, &1647 dimension_names = (/ 'turbine', 'time ' /), &1620 return_value = dom_def_var( nc_filename, & 1621 variable_name = variable_name, & 1622 dimension_names = (/ 'turbine', 'time ' /), & 1648 1623 output_type = 'real32' ) 1649 1650 return_value = dom_def_att( nc_filename, &1651 variable_name = variable_name, &1652 attribute_name = 'units', &1653 value = 'N' ) 1624 1625 return_value = dom_def_att( nc_filename, & 1626 variable_name = variable_name, & 1627 attribute_name = 'units', & 1628 value = 'N' ) 1654 1629 1655 1630 variable_name = 'rotor_torque' 1656 return_value = dom_def_var( nc_filename, &1657 variable_name = variable_name, &1658 dimension_names = (/ 'turbine', 'time ' /), &1631 return_value = dom_def_var( nc_filename, & 1632 variable_name = variable_name, & 1633 dimension_names = (/ 'turbine', 'time ' /), & 1659 1634 output_type = 'real32' ) 1660 1661 return_value = dom_def_att( nc_filename, &1662 variable_name = variable_name, &1663 attribute_name = 'units', &1635 1636 return_value = dom_def_att( nc_filename, & 1637 variable_name = variable_name, & 1638 attribute_name = 'units', & 1664 1639 value = 'Nm' ) 1665 1640 1666 1641 variable_name = 'wind_direction' 1667 return_value = dom_def_var( nc_filename, &1668 variable_name = variable_name, &1669 dimension_names = (/ 'turbine', 'time ' /), &1642 return_value = dom_def_var( nc_filename, & 1643 variable_name = variable_name, & 1644 dimension_names = (/ 'turbine', 'time ' /), & 1670 1645 output_type = 'real32' ) 1671 1672 return_value = dom_def_att( nc_filename, &1673 variable_name = variable_name, &1674 attribute_name = 'units', &1675 value = 'degrees' ) 1646 1647 return_value = dom_def_att( nc_filename, & 1648 variable_name = variable_name, & 1649 attribute_name = 'units', & 1650 value = 'degrees' ) 1676 1651 1677 1652 variable_name = 'yaw_angle' 1678 return_value = dom_def_var( nc_filename, &1679 variable_name = variable_name, &1680 dimension_names = (/ 'turbine', 'time ' /), &1653 return_value = dom_def_var( nc_filename, & 1654 variable_name = variable_name, & 1655 dimension_names = (/ 'turbine', 'time ' /), & 1681 1656 output_type = 'real32' ) 1682 1683 return_value = dom_def_att( nc_filename, &1684 variable_name = variable_name, &1685 attribute_name = 'units', &1686 value = 'degrees' ) 1657 1658 return_value = dom_def_att( nc_filename, & 1659 variable_name = variable_name, & 1660 attribute_name = 'units', & 1661 value = 'degrees' ) 1687 1662 1688 1663 ENDIF 1689 1664 END SUBROUTINE 1690 1691 !------------------------------------------------------------------------------ !1665 1666 !--------------------------------------------------------------------------------------------------! 1692 1667 ! Description: 1693 1668 ! ------------ 1694 !> Read in layout of the rotor blade , the lift and drag tables 1695 !> and the distribution of lift and drag tables along the blade 1696 !------------------------------------------------------------------------------! 1697 ! 1698 SUBROUTINE wtm_read_blade_tables 1699 1700 1701 IMPLICIT NONE 1702 1703 INTEGER(iwp) :: ii !< running index 1704 INTEGER(iwp) :: jj !< running index 1705 1706 INTEGER(iwp) :: ierrn !< 1707 1708 CHARACTER(200) :: chmess !< Read in string 1709 1710 INTEGER(iwp) :: dlen !< no. rows of local table 1711 INTEGER(iwp) :: dlenbl !< no. rows of cd, cl table 1712 INTEGER(iwp) :: ialpha !< table position of current alpha value 1713 INTEGER(iwp) :: iialpha !< 1714 INTEGER(iwp) :: iir !< 1715 INTEGER(iwp) :: radres !< radial resolution 1716 INTEGER(iwp) :: t1 !< no. of airfoil 1717 INTEGER(iwp) :: t2 !< no. of airfoil 1718 INTEGER(iwp) :: trow !< 1719 INTEGER(iwp) :: dlenbl_int !< no. rows of interpolated cd, cl tables 1720 1721 REAL(wp) :: alpha_attack_i !< 1722 REAL(wp) :: weight_a !< 1723 REAL(wp) :: weight_b !< 1724 1725 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint1 !< 1726 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint2 !< 1727 1728 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd_sel1 !< 1729 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd_sel2 !< 1730 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl_sel1 !< 1731 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl_sel2 !< 1732 REAL(wp), DIMENSION(:), ALLOCATABLE :: read_cl_cd !< read in var array 1733 1734 REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_attack_tab !< 1735 REAL(wp), DIMENSION(:), ALLOCATABLE :: trad1 !< 1736 REAL(wp), DIMENSION(:), ALLOCATABLE :: trad2 !< 1737 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cd_table !< 1738 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_table !< 1739 1740 ALLOCATE ( read_cl_cd(1:2*n_airfoils+1) ) 1741 1742 ! 1743 !-- Read in the distribution of lift and drag tables along the blade, the 1744 !-- layout of the rotor blade and the lift and drag tables: 1745 1746 OPEN ( 90, FILE='WTM_DATA', STATUS='OLD', FORM='FORMATTED', IOSTAT=ierrn ) 1747 1748 IF ( ierrn /= 0 ) THEN 1749 message_string = 'file WTM_DATA does not exist' 1750 CALL message( 'wtm_init', 'PA0464', 1, 2, 0, 6, 0 ) 1669 !> Read in layout of the rotor blade , the lift and drag tables and the distribution of lift and 1670 !> drag tables along the blade 1671 !--------------------------------------------------------------------------------------------------! 1672 ! 1673 SUBROUTINE wtm_read_blade_tables 1674 1675 1676 IMPLICIT NONE 1677 1678 CHARACTER(200) :: chmess !< Read in string 1679 1680 INTEGER(iwp) :: ii !< running index 1681 INTEGER(iwp) :: jj !< running index 1682 1683 INTEGER(iwp) :: ierrn !< 1684 1685 INTEGER(iwp) :: dlen !< no. rows of local table 1686 INTEGER(iwp) :: dlenbl !< no. rows of cd, cl table 1687 INTEGER(iwp) :: dlenbl_int !< no. rows of interpolated cd, cl tables 1688 INTEGER(iwp) :: ialpha !< table position of current alpha value 1689 INTEGER(iwp) :: iialpha !< 1690 INTEGER(iwp) :: iir !< 1691 INTEGER(iwp) :: radres !< radial resolution 1692 INTEGER(iwp) :: t1 !< no. of airfoil 1693 INTEGER(iwp) :: t2 !< no. of airfoil 1694 INTEGER(iwp) :: trow !< 1695 1696 1697 REAL(wp) :: alpha_attack_i !< 1698 REAL(wp) :: weight_a !< 1699 REAL(wp) :: weight_b !< 1700 1701 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint1 !< 1702 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint2 !< 1703 1704 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd_sel1 !< 1705 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cd_sel2 !< 1706 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl_sel1 !< 1707 REAL(wp), DIMENSION(:), ALLOCATABLE :: turb_cl_sel2 !< 1708 REAL(wp), DIMENSION(:), ALLOCATABLE :: read_cl_cd !< read in var array 1709 1710 REAL(wp), DIMENSION(:), ALLOCATABLE :: alpha_attack_tab !< 1711 REAL(wp), DIMENSION(:), ALLOCATABLE :: trad1 !< 1712 REAL(wp), DIMENSION(:), ALLOCATABLE :: trad2 !< 1713 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cd_table !< 1714 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: turb_cl_table !< 1715 1716 ALLOCATE ( read_cl_cd(1:2 * n_airfoils + 1) ) 1717 1718 ! 1719 !-- Read in the distribution of lift and drag tables along the blade, the layout of the rotor 1720 !-- blade and the lift and drag tables: 1721 OPEN ( 90, FILE='WTM_DATA', STATUS='OLD', FORM='FORMATTED', IOSTAT=ierrn ) 1722 1723 IF ( ierrn /= 0 ) THEN 1724 message_string = 'file WTM_DATA does not exist' 1725 CALL message( 'wtm_init', 'PA0464', 1, 2, 0, 6, 0 ) 1726 ENDIF 1727 ! 1728 !-- Read distribution table: 1729 dlen = 0 1730 1731 READ ( 90, '(3/)' ) 1732 1733 rloop3: DO 1734 READ ( 90, *, IOSTAT=ierrn ) chmess 1735 IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop3 1736 dlen = dlen + 1 1737 ENDDO rloop3 1738 1739 ALLOCATE( trad1(1:dlen), trad2(1:dlen), ttoint1(1:dlen), ttoint2(1:dlen) ) 1740 1741 DO jj = 1, dlen+1 1742 BACKSPACE ( 90, IOSTAT=ierrn ) 1743 ENDDO 1744 1745 DO jj = 1, dlen 1746 READ ( 90, * ) trad1(jj), trad2(jj), ttoint1(jj), ttoint2(jj) 1747 ENDDO 1748 1749 ! 1750 !-- Read layout table: 1751 dlen = 0 1752 1753 READ ( 90, '(3/)') 1754 1755 rloop1: DO 1756 READ ( 90, *, IOSTAT=ierrn ) chmess 1757 IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop1 1758 dlen = dlen + 1 1759 ENDDO rloop1 1760 1761 ALLOCATE( lrd(1:dlen), ard(1:dlen), crd(1:dlen) ) 1762 DO jj = 1, dlen + 1 1763 BACKSPACE ( 90, IOSTAT=ierrn ) 1764 ENDDO 1765 DO jj = 1, dlen 1766 READ ( 90, * ) lrd(jj), ard(jj), crd(jj) 1767 ENDDO 1768 1769 ! 1770 !-- Read tables (turb_cl(alpha),turb_cd(alpha) for the different profiles: 1771 dlen = 0 1772 1773 READ ( 90, '(3/)' ) 1774 1775 rloop2: DO 1776 READ ( 90, *, IOSTAT=ierrn ) chmess 1777 IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop2 1778 dlen = dlen + 1 1779 ENDDO rloop2 1780 1781 ALLOCATE( alpha_attack_tab(1:dlen), turb_cl_table(1:dlen,1:n_airfoils), & 1782 turb_cd_table(1:dlen,1:n_airfoils) ) 1783 1784 DO jj = 1, dlen + 1 1785 BACKSPACE ( 90, IOSTAT=ierrn ) 1786 ENDDO 1787 1788 DO jj = 1, dlen 1789 READ ( 90, * ) read_cl_cd 1790 alpha_attack_tab(jj) = read_cl_cd(1) 1791 DO ii = 1, n_airfoils 1792 turb_cl_table(jj,ii) = read_cl_cd(ii * 2) 1793 turb_cd_table(jj,ii) = read_cl_cd(ii * 2 + 1) 1794 ENDDO 1795 1796 ENDDO 1797 1798 dlenbl = dlen 1799 1800 CLOSE ( 90 ) 1801 1802 ! 1803 !-- For each possible radial position (resolution: 0.1 m --> 631 values if rotor_radius(1)=63m) 1804 !-- and each possible angle of attack (resolution: 0.1 degrees --> 3601 values!) determine the 1805 !-- lift and drag coefficient by interpolating between the tabulated values of each table 1806 !-- (interpolate to current angle of attack) and between the tables (interpolate to current 1807 !-- radial position): 1808 ALLOCATE( turb_cl_sel1(1:dlenbl) ) 1809 ALLOCATE( turb_cl_sel2(1:dlenbl) ) 1810 ALLOCATE( turb_cd_sel1(1:dlenbl) ) 1811 ALLOCATE( turb_cd_sel2(1:dlenbl) ) 1812 1813 radres = INT( rotor_radius(1) * 10.0_wp ) + 1_iwp 1814 dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp 1815 1816 ALLOCATE( turb_cl_tab(1:dlenbl_int,1:radres) ) 1817 ALLOCATE( turb_cd_tab(1:dlenbl_int,1:radres) ) 1818 1819 DO iir = 1, radres ! loop over radius 1820 1821 cur_r = ( iir - 1_iwp ) * 0.1_wp 1822 ! 1823 !-- Find position in table 1: 1824 lct = MINLOC( ABS( trad1 - cur_r ) ) 1825 ! lct(1) = lct(1) 1826 1827 IF ( ( trad1(lct(1)) - cur_r ) > 0.0 ) THEN 1828 lct(1) = lct(1) - 1 1751 1829 ENDIF 1752 ! 1753 !-- Read distribution table: 1754 1755 dlen = 0 1756 1757 READ ( 90, '(3/)' ) 1758 1759 rloop3: DO 1760 READ ( 90, *, IOSTAT=ierrn ) chmess 1761 IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop3 1762 dlen = dlen + 1 1763 ENDDO rloop3 1764 1765 ALLOCATE( trad1(1:dlen), trad2(1:dlen), ttoint1(1:dlen), ttoint2(1:dlen)) 1766 1767 DO jj = 1,dlen+1 1768 BACKSPACE ( 90, IOSTAT=ierrn ) 1769 ENDDO 1770 1771 DO jj = 1,dlen 1772 READ ( 90, * ) trad1(jj), trad2(jj), ttoint1(jj), ttoint2(jj) 1773 ENDDO 1774 1775 ! 1776 !-- Read layout table: 1777 1778 dlen = 0 1779 1780 READ ( 90, '(3/)') 1781 1782 rloop1: DO 1783 READ ( 90, *, IOSTAT=ierrn ) chmess 1784 IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop1 1785 dlen = dlen + 1 1786 ENDDO rloop1 1787 1788 ALLOCATE( lrd(1:dlen), ard(1:dlen), crd(1:dlen) ) 1789 DO jj = 1, dlen+1 1790 BACKSPACE ( 90, IOSTAT=ierrn ) 1791 ENDDO 1792 DO jj = 1, dlen 1793 READ ( 90, * ) lrd(jj), ard(jj), crd(jj) 1794 ENDDO 1795 1796 ! 1797 !-- Read tables (turb_cl(alpha),turb_cd(alpha) for the different profiles: 1798 1799 dlen = 0 1800 1801 READ ( 90, '(3/)' ) 1802 1803 rloop2: DO 1804 READ ( 90, *, IOSTAT=ierrn ) chmess 1805 IF ( ierrn < 0 .OR. chmess == '#' .OR. chmess == '') EXIT rloop2 1806 dlen = dlen + 1 1807 ENDDO rloop2 1808 1809 ALLOCATE( alpha_attack_tab(1:dlen), turb_cl_table(1:dlen,1:n_airfoils), & 1810 turb_cd_table(1:dlen,1:n_airfoils) ) 1811 1812 DO jj = 1,dlen+1 1813 BACKSPACE ( 90, IOSTAT=ierrn ) 1814 ENDDO 1815 1816 DO jj = 1,dlen 1817 READ ( 90, * ) read_cl_cd 1818 alpha_attack_tab(jj) = read_cl_cd(1) 1819 DO ii= 1, n_airfoils 1820 turb_cl_table(jj,ii) = read_cl_cd(ii*2) 1821 turb_cd_table(jj,ii) = read_cl_cd(ii*2+1) 1830 1831 trow = lct(1) 1832 ! 1833 !-- Calculate weights for radius interpolation: 1834 weight_a = ( trad2(trow) - cur_r ) / ( trad2(trow) - trad1(trow) ) 1835 weight_b = ( cur_r - trad1(trow) ) / ( trad2(trow) - trad1(trow) ) 1836 t1 = ttoint1(trow) 1837 t2 = ttoint2(trow) 1838 1839 IF ( t1 == t2 ) THEN ! if both are the same, the weights are NaN 1840 weight_a = 0.5_wp ! then do interpolate in between same twice 1841 weight_b = 0.5_wp ! using 0.5 as weight 1842 ENDIF 1843 1844 IF ( t1 == 0 .AND. t2 == 0 ) THEN 1845 turb_cd_sel1 = 0.0_wp 1846 turb_cd_sel2 = 0.0_wp 1847 turb_cl_sel1 = 0.0_wp 1848 turb_cl_sel2 = 0.0_wp 1849 1850 turb_cd_tab(1,iir) = 0.0_wp ! For -180 degrees (iialpha=1) the values 1851 turb_cl_tab(1,iir) = 0.0_wp ! for each radius has to be set 1852 ! explicitly 1853 ELSE 1854 turb_cd_sel1 = turb_cd_table(:,t1) 1855 turb_cd_sel2 = turb_cd_table(:,t2) 1856 turb_cl_sel1 = turb_cl_table(:,t1) 1857 turb_cl_sel2 = turb_cl_table(:,t2) 1858 ! 1859 !-- For -180 degrees (iialpha=1) the values for each radius has to be set explicitly: 1860 turb_cd_tab(1,iir) = ( weight_a * turb_cd_table(1,t1) + weight_b * turb_cd_table(1,t2) ) 1861 turb_cl_tab(1,iir) = ( weight_a * turb_cl_table(1,t1) + weight_b * turb_cl_table(1,t2) ) 1862 ENDIF 1863 1864 DO iialpha = 2, dlenbl_int ! loop over angles 1865 1866 alpha_attack_i = -180.0_wp + REAL( iialpha-1 ) * accu_cl_cd_tab 1867 ialpha = 1 1868 1869 DO WHILE ( ( alpha_attack_i > alpha_attack_tab(ialpha) ) .AND. ( ialpha < dlen ) ) 1870 ialpha = ialpha + 1 1822 1871 ENDDO 1823 1872 1824 ENDDO 1825 1826 dlenbl = dlen 1827 1828 CLOSE ( 90 ) 1829 1830 ! 1831 !-- For each possible radial position (resolution: 0.1 m --> 631 values if rotor_radius(1)=63m) 1832 !-- and each possible angle of attack (resolution: 0.1 degrees --> 3601 values!) 1833 !-- determine the lift and drag coefficient by interpolating between the 1834 !-- tabulated values of each table (interpolate to current angle of attack) 1835 !-- and between the tables (interpolate to current radial position): 1836 1837 ALLOCATE( turb_cl_sel1(1:dlenbl) ) 1838 ALLOCATE( turb_cl_sel2(1:dlenbl) ) 1839 ALLOCATE( turb_cd_sel1(1:dlenbl) ) 1840 ALLOCATE( turb_cd_sel2(1:dlenbl) ) 1841 1842 radres = INT( rotor_radius(1) * 10.0_wp ) + 1_iwp 1843 dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp 1844 1845 ALLOCATE( turb_cl_tab(1:dlenbl_int,1:radres) ) 1846 ALLOCATE( turb_cd_tab(1:dlenbl_int,1:radres) ) 1847 1848 DO iir = 1, radres ! loop over radius 1849 1850 cur_r = ( iir - 1_iwp ) * 0.1_wp 1851 ! 1852 !-- Find position in table 1 1853 lct = MINLOC( ABS( trad1 - cur_r ) ) 1854 ! lct(1) = lct(1) 1855 1856 IF ( ( trad1(lct(1)) - cur_r ) > 0.0 ) THEN 1857 lct(1) = lct(1) - 1 1858 ENDIF 1859 1860 trow = lct(1) 1861 ! 1862 !-- Calculate weights for radius interpolation 1863 weight_a = ( trad2(trow) - cur_r ) / ( trad2(trow) - trad1(trow) ) 1864 weight_b = ( cur_r - trad1(trow) ) / ( trad2(trow) - trad1(trow) ) 1865 t1 = ttoint1(trow) 1866 t2 = ttoint2(trow) 1867 1868 IF ( t1 == t2 ) THEN ! if both are the same, the weights are NaN 1869 weight_a = 0.5_wp ! then do interpolate in between same twice 1870 weight_b = 0.5_wp ! using 0.5 as weight 1871 ENDIF 1872 1873 IF ( t1 == 0 .AND. t2 == 0 ) THEN 1874 turb_cd_sel1 = 0.0_wp 1875 turb_cd_sel2 = 0.0_wp 1876 turb_cl_sel1 = 0.0_wp 1877 turb_cl_sel2 = 0.0_wp 1878 1879 turb_cd_tab(1,iir) = 0.0_wp ! For -180 degrees (iialpha=1) the values 1880 turb_cl_tab(1,iir) = 0.0_wp ! for each radius has to be set 1881 ! explicitly 1882 ELSE 1883 turb_cd_sel1 = turb_cd_table(:,t1) 1884 turb_cd_sel2 = turb_cd_table(:,t2) 1885 turb_cl_sel1 = turb_cl_table(:,t1) 1886 turb_cl_sel2 = turb_cl_table(:,t2) 1887 ! 1888 !-- For -180 degrees (iialpha=1) the values for each radius has to be set 1889 !-- explicitly 1890 turb_cd_tab(1,iir) = ( weight_a * turb_cd_table(1,t1) + weight_b & 1891 * turb_cd_table(1,t2) ) 1892 turb_cl_tab(1,iir) = ( weight_a * turb_cl_table(1,t1) + weight_b & 1893 * turb_cl_table(1,t2) ) 1894 ENDIF 1895 1896 DO iialpha = 2, dlenbl_int ! loop over angles 1897 1898 alpha_attack_i = -180.0_wp + REAL( iialpha-1 ) * accu_cl_cd_tab 1899 ialpha = 1 1900 1901 DO WHILE ( ( alpha_attack_i > alpha_attack_tab(ialpha) ) .AND. (ialpha < dlen ) ) 1902 ialpha = ialpha + 1 1903 ENDDO 1904 1905 ! 1906 !-- Interpolation of lift and drag coefficiencts on fine grid of radius 1907 !-- segments and angles of attack 1908 1909 turb_cl_tab(iialpha,iir) = ( alpha_attack_tab(ialpha) - & 1910 alpha_attack_i ) / & 1911 ( alpha_attack_tab(ialpha) - & 1912 alpha_attack_tab(ialpha-1) ) * & 1913 ( weight_a * turb_cl_sel1(ialpha-1) + & 1914 weight_b * turb_cl_sel2(ialpha-1) ) +& 1915 ( alpha_attack_i - & 1916 alpha_attack_tab(ialpha-1) ) / & 1917 ( alpha_attack_tab(ialpha) - & 1918 alpha_attack_tab(ialpha-1) ) * & 1919 ( weight_a * turb_cl_sel1(ialpha) + & 1920 weight_b * turb_cl_sel2(ialpha) ) 1921 turb_cd_tab(iialpha,iir) = ( alpha_attack_tab(ialpha) - & 1922 alpha_attack_i ) / & 1923 ( alpha_attack_tab(ialpha) - & 1924 alpha_attack_tab(ialpha-1) ) * & 1925 ( weight_a * turb_cd_sel1(ialpha-1) + & 1926 weight_b * turb_cd_sel2(ialpha-1) ) +& 1927 ( alpha_attack_i - & 1928 alpha_attack_tab(ialpha-1) ) / & 1929 ( alpha_attack_tab(ialpha) - & 1930 alpha_attack_tab(ialpha-1) ) * & 1931 ( weight_a * turb_cd_sel1(ialpha) + & 1932 weight_b * turb_cd_sel2(ialpha) ) 1933 1934 ENDDO ! end loop over angles of attack 1935 1936 ENDDO ! end loop over radius 1937 1938 1939 END SUBROUTINE wtm_read_blade_tables 1940 1941 1942 !------------------------------------------------------------------------------! 1873 ! 1874 !-- Interpolation of lift and drag coefficiencts on fine grid of radius segments and angles 1875 !-- of attack: 1876 turb_cl_tab(iialpha,iir) = ( alpha_attack_tab(ialpha) - alpha_attack_i ) / & 1877 ( alpha_attack_tab(ialpha) - alpha_attack_tab(ialpha-1) ) * & 1878 ( weight_a * turb_cl_sel1(ialpha-1) + & 1879 weight_b * turb_cl_sel2(ialpha-1) ) + & 1880 ( alpha_attack_i - alpha_attack_tab(ialpha-1) ) / & 1881 ( alpha_attack_tab(ialpha) - alpha_attack_tab(ialpha-1) ) * & 1882 ( weight_a * turb_cl_sel1(ialpha) + & 1883 weight_b * turb_cl_sel2(ialpha) ) 1884 turb_cd_tab(iialpha,iir) = ( alpha_attack_tab(ialpha) - alpha_attack_i ) / & 1885 ( alpha_attack_tab(ialpha) - alpha_attack_tab(ialpha-1) ) * & 1886 ( weight_a * turb_cd_sel1(ialpha-1) + & 1887 weight_b * turb_cd_sel2(ialpha-1) ) + & 1888 ( alpha_attack_i - alpha_attack_tab(ialpha-1) ) / & 1889 ( alpha_attack_tab(ialpha) - alpha_attack_tab(ialpha-1) ) * & 1890 ( weight_a * turb_cd_sel1(ialpha) + & 1891 weight_b * turb_cd_sel2(ialpha) ) 1892 1893 ENDDO ! end loop over angles of attack 1894 1895 ENDDO ! end loop over radius 1896 1897 1898 END SUBROUTINE wtm_read_blade_tables 1899 1900 1901 !--------------------------------------------------------------------------------------------------! 1943 1902 ! Description: 1944 1903 ! ------------ 1945 !> The projection matrix for the coordinate system of therotor disc in respect 1946 !> to the yaw and tilt angle of the rotor is calculated 1947 !------------------------------------------------------------------------------! 1948 SUBROUTINE wtm_rotate_rotor( inot ) 1949 1950 1951 IMPLICIT NONE 1952 1953 INTEGER(iwp) :: inot 1954 ! 1955 !-- Calculation of the rotation matrix for the application of the tilt to 1956 !-- the rotors 1957 rot_eigen_rad(1) = SIN( yaw_angle(inot) ) ! x-component of the radial eigenvector 1958 rot_eigen_rad(2) = COS( yaw_angle(inot) ) ! y-component of the radial eigenvector 1959 rot_eigen_rad(3) = 0.0_wp ! z-component of the radial eigenvector 1960 1961 rot_eigen_azi(1) = 0.0_wp ! x-component of the azimuth eigenvector 1962 rot_eigen_azi(2) = 0.0_wp ! y-component of the azimuth eigenvector 1963 rot_eigen_azi(3) = 1.0_wp ! z-component of the azimuth eigenvector 1964 1965 rot_eigen_nor(1) = COS( yaw_angle(inot) ) ! x-component of the normal eigenvector 1966 rot_eigen_nor(2) = -SIN( yaw_angle(inot) ) ! y-component of the normal eigenvector 1967 rot_eigen_nor(3) = 0.0_wp ! z-component of the normal eigenvector 1968 1969 ! 1970 !-- Calculation of the coordinate transformation matrix to apply a tilt to 1971 !-- the rotor. If tilt = 0, rot_coord_trans is a unit matrix. 1972 1973 rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2 * & 1974 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1975 rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2) * & 1976 ( 1.0_wp - COS( tilt_angle ) ) - & 1977 rot_eigen_rad(3) * SIN( tilt_angle ) 1978 rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3) * & 1979 ( 1.0_wp - COS( tilt_angle ) ) + & 1980 rot_eigen_rad(2) * SIN( tilt_angle ) 1981 rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1) * & 1982 ( 1.0_wp - COS( tilt_angle ) ) + & 1983 rot_eigen_rad(3) * SIN( tilt_angle ) 1984 rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2 * & 1985 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1986 rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3) * & 1987 ( 1.0_wp - COS( tilt_angle ) ) - & 1988 rot_eigen_rad(1) * SIN( tilt_angle ) 1989 rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1) * & 1990 ( 1.0_wp - COS( tilt_angle ) ) - & 1991 rot_eigen_rad(2) * SIN( tilt_angle ) 1992 rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2) * & 1993 ( 1.0_wp - COS( tilt_angle ) ) + & 1994 rot_eigen_rad(1) * SIN( tilt_angle ) 1995 rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2 * & 1996 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1997 1998 ! 1999 !-- Vectors for the Transformation of forces from the rotor's spheric 2000 !-- coordinate system to the cartesian coordinate system 2001 rotx(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_nor ) 2002 roty(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_rad ) 2003 rotz(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_azi ) 2004 2005 END SUBROUTINE wtm_rotate_rotor 2006 2007 2008 !------------------------------------------------------------------------------! 1904 !> The projection matrix for the coordinate system of therotor disc in respect to the yaw and tilt 1905 !> angle of the rotor is calculated 1906 !--------------------------------------------------------------------------------------------------! 1907 SUBROUTINE wtm_rotate_rotor( inot ) 1908 1909 1910 IMPLICIT NONE 1911 1912 INTEGER(iwp) :: inot 1913 ! 1914 !-- Calculation of the rotation matrix for the application of the tilt to the rotors: 1915 rot_eigen_rad(1) = SIN( yaw_angle(inot) ) ! x-component of the radial eigenvector 1916 rot_eigen_rad(2) = COS( yaw_angle(inot) ) ! y-component of the radial eigenvector 1917 rot_eigen_rad(3) = 0.0_wp ! z-component of the radial eigenvector 1918 1919 rot_eigen_azi(1) = 0.0_wp ! x-component of the azimuth eigenvector 1920 rot_eigen_azi(2) = 0.0_wp ! y-component of the azimuth eigenvector 1921 rot_eigen_azi(3) = 1.0_wp ! z-component of the azimuth eigenvector 1922 1923 rot_eigen_nor(1) = COS( yaw_angle(inot) ) ! x-component of the normal eigenvector 1924 rot_eigen_nor(2) = -SIN( yaw_angle(inot) ) ! y-component of the normal eigenvector 1925 rot_eigen_nor(3) = 0.0_wp ! z-component of the normal eigenvector 1926 1927 ! 1928 !-- Calculation of the coordinate transformation matrix to apply a tilt to the rotor. 1929 !-- If tilt = 0, rot_coord_trans is a unit matrix. 1930 rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2 * & 1931 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1932 rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2) * & 1933 ( 1.0_wp - COS( tilt_angle ) ) - & 1934 rot_eigen_rad(3) * SIN( tilt_angle ) 1935 rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3) * & 1936 ( 1.0_wp - COS( tilt_angle ) ) + & 1937 rot_eigen_rad(2) * SIN( tilt_angle ) 1938 rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1) * & 1939 ( 1.0_wp - COS( tilt_angle ) ) + & 1940 rot_eigen_rad(3) * SIN( tilt_angle ) 1941 rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2 * & 1942 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1943 rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3) * & 1944 ( 1.0_wp - COS( tilt_angle ) ) - & 1945 rot_eigen_rad(1) * SIN( tilt_angle ) 1946 rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1) * & 1947 ( 1.0_wp - COS( tilt_angle ) ) - & 1948 rot_eigen_rad(2) * SIN( tilt_angle ) 1949 rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2) * & 1950 ( 1.0_wp - COS( tilt_angle ) ) + & 1951 rot_eigen_rad(1) * SIN( tilt_angle ) 1952 rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2 * & 1953 ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 1954 1955 ! 1956 !-- Vectors for the Transformation of forces from the rotor's spheric coordinate system to the 1957 !-- cartesian coordinate system: 1958 rotx(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_nor ) 1959 roty(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_rad ) 1960 rotz(inot,:) = MATMUL( rot_coord_trans(inot,:,:), rot_eigen_azi ) 1961 1962 END SUBROUTINE wtm_rotate_rotor 1963 1964 1965 !--------------------------------------------------------------------------------------------------! 2009 1966 ! Description: 2010 1967 ! ------------ 2011 1968 !> Calculation of the forces generated by the wind turbine 2012 !------------------------------------------------------------------------------! 2013 SUBROUTINE wtm_forces 2014 2015 2016 IMPLICIT NONE 2017 2018 2019 INTEGER(iwp) :: i, j, k !< loop indices 2020 INTEGER(iwp) :: inot !< turbine loop index (turbine id) 2021 INTEGER(iwp) :: iialpha, iir !< 2022 INTEGER(iwp) :: rseg !< 2023 INTEGER(iwp) :: ring !< 2024 INTEGER(iwp) :: ii, jj, kk !< 2025 2026 REAL(wp) :: flag !< flag to mask topography grid points 2027 REAL(wp) :: sin_rot, cos_rot !< 2028 REAL(wp) :: sin_yaw, cos_yaw !< 2029 2030 REAL(wp) :: aa, bb, cc, dd !< interpolation distances 2031 REAL(wp) :: gg !< interpolation volume var 2032 2033 REAL(wp) :: dist_u_3d, dist_v_3d, dist_w_3d !< smearing distances 2034 2035 2036 ! 2037 ! Variables for pitch control 2038 LOGICAL :: pitch_sw = .FALSE. 2039 2040 INTEGER(iwp), DIMENSION(1) :: lct = 0 2041 REAL(wp), DIMENSION(1) :: rad_d = 0.0_wp 2042 2043 REAL(wp) :: tl_factor !< factor for tip loss correction 2044 2045 2046 CALL cpu_log( log_point_s(61), 'wtm_forces', 'start' ) 2047 2048 2049 IF ( time_since_reference_point >= time_turbine_on ) THEN 2050 2051 ! 2052 !-- Set forces to zero for each new time step: 2053 thrust(:,:,:) = 0.0_wp 2054 torque_y(:,:,:) = 0.0_wp 2055 torque_z(:,:,:) = 0.0_wp 2056 torque_total(:) = 0.0_wp 2057 rot_tend_x(:,:,:) = 0.0_wp 2058 rot_tend_y(:,:,:) = 0.0_wp 2059 rot_tend_z(:,:,:) = 0.0_wp 2060 thrust_rotor(:) = 0.0_wp 2061 ! 2062 !-- Loop over number of turbines: 2063 DO inot = 1, n_turbines 2064 2065 cos_yaw = COS(yaw_angle(inot)) 2066 sin_yaw = SIN(yaw_angle(inot)) 2067 ! 2068 !-- Loop over rings of each turbine: 2069 2070 !$OMP PARALLEL PRIVATE (ring, rseg, thrust_seg, torque_seg_y, torque_seg_z, sin_rot, & 2071 !$OMP& cos_rot, re, rbx, rby, rbz, ii, jj, kk, aa, bb, cc, dd, gg) 2072 !$OMP DO 2073 DO ring = 1, nrings(inot) 2074 2075 thrust_seg(:) = 0.0_wp 2076 torque_seg_y(:) = 0.0_wp 2077 torque_seg_z(:) = 0.0_wp 2078 ! 2079 !-- Determine distance between each ring (center) and the hub: 2080 cur_r = (ring - 0.5_wp) * delta_r(inot) 2081 2082 ! 2083 !-- Loop over segments of each ring of each turbine: 2084 DO rseg = 1, nsegs(ring,inot) 2085 ! 2086 !-- !-----------------------------------------------------------! 2087 !-- !-- Determine coordinates of the ring segments --! 2088 !-- !-----------------------------------------------------------! 2089 ! 2090 !-- Determine angle of ring segment towards zero degree angle of 2091 !-- rotor system (at zero degree rotor direction vectors aligned 2092 !-- with y-axis): 2093 phi_rotor = rseg * 2.0_wp * pi / nsegs(ring,inot) 2094 cos_rot = COS( phi_rotor ) 2095 sin_rot = SIN( phi_rotor ) 2096 2097 !-- Now the direction vectors can be determined with respect to 2098 !-- the yaw and tilt angle: 2099 re(1) = cos_rot * sin_yaw 2100 re(2) = cos_rot * cos_yaw 2101 re(3) = sin_rot 2102 2103 rote = MATMUL( rot_coord_trans(inot,:,:), re ) 2104 ! 2105 !-- Coordinates of the single segments (center points): 2106 rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1) 2107 rby(ring,rseg) = hub_y(inot) + cur_r * rote(2) 2108 rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3) 2109 2110 !-- !-----------------------------------------------------------! 2111 !-- !-- Interpolation of the velocity components from the --! 2112 !-- !-- cartesian grid point to the coordinates of each ring --! 2113 !-- !-- segment (follows a method used in the particle model) --! 2114 !-- !-----------------------------------------------------------! 2115 2116 u_int(inot,ring,rseg) = 0.0_wp 2117 u_int_1_l(inot,ring,rseg) = 0.0_wp 2118 2119 v_int(inot,ring,rseg) = 0.0_wp 2120 v_int_1_l(inot,ring,rseg) = 0.0_wp 2121 2122 w_int(inot,ring,rseg) = 0.0_wp 2123 w_int_1_l(inot,ring,rseg) = 0.0_wp 2124 2125 ! 2126 !-- Interpolation of the u-component: 2127 2128 ii = rbx(ring,rseg) * ddx 2129 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 2130 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 2131 ! 2132 !-- Interpolate only if all required information is available on 2133 !-- the current PE: 2134 IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN 2135 IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN 2136 2137 aa = ( ( ii + 1 ) * dx - rbx(ring,rseg) ) * & 2138 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2139 bb = ( rbx(ring,rseg) - ii * dx ) * & 2140 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2141 cc = ( ( ii+1 ) * dx - rbx(ring,rseg) ) * & 2142 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2143 dd = ( rbx(ring,rseg) - ii * dx ) * & 2144 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2145 gg = dx * dy 2146 2147 u_int_l = ( aa * u(kk,jj,ii) + & 2148 bb * u(kk,jj,ii+1) + & 2149 cc * u(kk,jj+1,ii) + & 2150 dd * u(kk,jj+1,ii+1) & 2151 ) / gg 2152 2153 u_int_u = ( aa * u(kk+1,jj,ii) + & 2154 bb * u(kk+1,jj,ii+1) + & 2155 cc * u(kk+1,jj+1,ii) + & 2156 dd * u(kk+1,jj+1,ii+1) & 2157 ) / gg 2158 2159 u_int_1_l(inot,ring,rseg) = u_int_l + & 2160 ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 2161 ( u_int_u - u_int_l ) 2162 2163 ELSE 2164 u_int_1_l(inot,ring,rseg) = 0.0_wp 2165 ENDIF 1969 !--------------------------------------------------------------------------------------------------! 1970 SUBROUTINE wtm_forces 1971 1972 1973 IMPLICIT NONE 1974 1975 1976 INTEGER(iwp) :: i, j, k !< loop indices 1977 INTEGER(iwp) :: ii, jj, kk !< 1978 INTEGER(iwp) :: inot !< turbine loop index (turbine id) 1979 INTEGER(iwp) :: iialpha, iir !< 1980 INTEGER(iwp) :: rseg !< 1981 INTEGER(iwp) :: ring !< 1982 1983 REAL(wp) :: flag !< flag to mask topography grid points 1984 REAL(wp) :: sin_rot, cos_rot !< 1985 REAL(wp) :: sin_yaw, cos_yaw !< 1986 1987 REAL(wp) :: aa, bb, cc, dd !< interpolation distances 1988 REAL(wp) :: gg !< interpolation volume var 1989 1990 REAL(wp) :: dist_u_3d, dist_v_3d, dist_w_3d !< smearing distances 1991 1992 1993 ! 1994 !-- Variables for pitch control: 1995 INTEGER(iwp), DIMENSION(1) :: lct = 0 1996 1997 LOGICAL :: pitch_sw = .FALSE. 1998 1999 REAL(wp), DIMENSION(1) :: rad_d = 0.0_wp 2000 2001 REAL(wp) :: tl_factor !< factor for tip loss correction 2002 2003 2004 CALL cpu_log( log_point_s(61), 'wtm_forces', 'start' ) 2005 2006 2007 IF ( time_since_reference_point >= time_turbine_on ) THEN 2008 2009 ! 2010 !-- Set forces to zero for each new time step: 2011 thrust(:,:,:) = 0.0_wp 2012 torque_y(:,:,:) = 0.0_wp 2013 torque_z(:,:,:) = 0.0_wp 2014 torque_total(:) = 0.0_wp 2015 rot_tend_x(:,:,:) = 0.0_wp 2016 rot_tend_y(:,:,:) = 0.0_wp 2017 rot_tend_z(:,:,:) = 0.0_wp 2018 thrust_rotor(:) = 0.0_wp 2019 ! 2020 !-- Loop over number of turbines: 2021 DO inot = 1, n_turbines 2022 2023 cos_yaw = COS(yaw_angle(inot)) 2024 sin_yaw = SIN(yaw_angle(inot)) 2025 ! 2026 !-- Loop over rings of each turbine: 2027 !$OMP PARALLEL PRIVATE (ring, rseg, thrust_seg, torque_seg_y, torque_seg_z, sin_rot, & 2028 !$OMP& cos_rot, re, rbx, rby, rbz, ii, jj, kk, aa, bb, cc, dd, gg) 2029 !$OMP DO 2030 DO ring = 1, nrings(inot) 2031 2032 thrust_seg(:) = 0.0_wp 2033 torque_seg_y(:) = 0.0_wp 2034 torque_seg_z(:) = 0.0_wp 2035 ! 2036 !-- Determine distance between each ring (center) and the hub: 2037 cur_r = (ring - 0.5_wp) * delta_r(inot) 2038 2039 ! 2040 !-- Loop over segments of each ring of each turbine: 2041 DO rseg = 1, nsegs(ring,inot) 2042 ! 2043 !-- !----------------------------------------------------------------------------------! 2044 !-- !-- Determine coordinates of the ring segments --! 2045 !-- !----------------------------------------------------------------------------------! 2046 ! 2047 !-- Determine angle of ring segment towards zero degree angle of rotor system 2048 !-- (at zero degree rotor direction vectors aligned with y-axis): 2049 phi_rotor = rseg * 2.0_wp * pi / nsegs(ring,inot) 2050 cos_rot = COS( phi_rotor ) 2051 sin_rot = SIN( phi_rotor ) 2052 2053 !-- Now the direction vectors can be determined with respect to the yaw and tilt angle: 2054 re(1) = cos_rot * sin_yaw 2055 re(2) = cos_rot * cos_yaw 2056 re(3) = sin_rot 2057 2058 rote = MATMUL( rot_coord_trans(inot,:,:), re ) 2059 ! 2060 !-- Coordinates of the single segments (center points): 2061 rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1) 2062 rby(ring,rseg) = hub_y(inot) + cur_r * rote(2) 2063 rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3) 2064 2065 !-- !----------------------------------------------------------------------------------! 2066 !-- !-- Interpolation of the velocity components from the cartesian grid point to --! 2067 !-- !-- the coordinates of each ring segment (follows a method used in the --! 2068 !-- !-- particle model) --! 2069 !-- !----------------------------------------------------------------------------------! 2070 2071 u_int(inot,ring,rseg) = 0.0_wp 2072 u_int_1_l(inot,ring,rseg) = 0.0_wp 2073 2074 v_int(inot,ring,rseg) = 0.0_wp 2075 v_int_1_l(inot,ring,rseg) = 0.0_wp 2076 2077 w_int(inot,ring,rseg) = 0.0_wp 2078 w_int_1_l(inot,ring,rseg) = 0.0_wp 2079 2080 ! 2081 !-- Interpolation of the u-component: 2082 ii = rbx(ring,rseg) * ddx 2083 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 2084 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 2085 ! 2086 !-- Interpolate only if all required information is available on the current PE: 2087 IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN 2088 IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN 2089 2090 aa = ( ( ii + 1 ) * dx - rbx(ring,rseg) ) * & 2091 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2092 bb = ( rbx(ring,rseg) - ii * dx ) * & 2093 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2094 cc = ( ( ii+1 ) * dx - rbx(ring,rseg) ) * & 2095 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2096 dd = ( rbx(ring,rseg) - ii * dx ) * & 2097 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2098 gg = dx * dy 2099 2100 u_int_l = ( aa * u(kk,jj,ii) + bb * u(kk,jj,ii+1) + cc * u(kk,jj+1,ii) & 2101 + dd * u(kk,jj+1,ii+1) ) / gg 2102 2103 u_int_u = ( aa * u(kk+1,jj,ii) + bb * u(kk+1,jj,ii+1) + cc * u(kk+1,jj+1,ii) & 2104 + dd * u(kk+1,jj+1,ii+1) ) / gg 2105 2106 u_int_1_l(inot,ring,rseg) = u_int_l + ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 2107 ( u_int_u - u_int_l ) 2108 2166 2109 ELSE 2167 2110 u_int_1_l(inot,ring,rseg) = 0.0_wp 2168 2111 ENDIF 2169 2112 2170 2171 ! 2172 !-- Interpolation of the v-component: 2173 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 2174 jj = rby(ring,rseg) * ddy 2175 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 2176 ! 2177 !-- Interpolate only if all required information is available on 2178 !-- the current PE: 2179 IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN 2180 IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN 2181 2182 aa = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2183 ( ( jj + 1 ) * dy - rby(ring,rseg) ) 2184 bb = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2185 ( ( jj + 1 ) * dy - rby(ring,rseg) ) 2186 cc = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2187 ( rby(ring,rseg) - jj * dy ) 2188 dd = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2189 ( rby(ring,rseg) - jj * dy ) 2190 gg = dx * dy 2191 2192 v_int_l = ( aa * v(kk,jj,ii) + & 2193 bb * v(kk,jj,ii+1) + & 2194 cc * v(kk,jj+1,ii) + & 2195 dd * v(kk,jj+1,ii+1) & 2196 ) / gg 2197 2198 v_int_u = ( aa * v(kk+1,jj,ii) + & 2199 bb * v(kk+1,jj,ii+1) + & 2200 cc * v(kk+1,jj+1,ii) + & 2201 dd * v(kk+1,jj+1,ii+1) & 2202 ) / gg 2203 2204 v_int_1_l(inot,ring,rseg) = v_int_l + & 2205 ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 2206 ( v_int_u - v_int_l ) 2207 2208 ELSE 2209 v_int_1_l(inot,ring,rseg) = 0.0_wp 2210 ENDIF 2113 ELSE 2114 u_int_1_l(inot,ring,rseg) = 0.0_wp 2115 ENDIF 2116 2117 2118 ! 2119 !-- Interpolation of the v-component: 2120 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 2121 jj = rby(ring,rseg) * ddy 2122 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 2123 ! 2124 !-- Interpolate only if all required information is available on the current PE: 2125 IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN 2126 IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN 2127 2128 aa = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2129 ( ( jj + 1 ) * dy - rby(ring,rseg) ) 2130 bb = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2131 ( ( jj + 1 ) * dy - rby(ring,rseg) ) 2132 cc = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2133 ( rby(ring,rseg) - jj * dy ) 2134 dd = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2135 ( rby(ring,rseg) - jj * dy ) 2136 gg = dx * dy 2137 2138 v_int_l = ( aa * v(kk,jj,ii) + bb * v(kk,jj,ii+1) + cc * v(kk,jj+1,ii) & 2139 + dd * v(kk,jj+1,ii+1) ) / gg 2140 2141 v_int_u = ( aa * v(kk+1,jj,ii) + bb * v(kk+1,jj,ii+1) + cc * v(kk+1,jj+1,ii) & 2142 + dd * v(kk+1,jj+1,ii+1) ) / gg 2143 2144 v_int_1_l(inot,ring,rseg) = v_int_l + ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 2145 ( v_int_u - v_int_l ) 2146 2211 2147 ELSE 2212 2148 v_int_1_l(inot,ring,rseg) = 0.0_wp 2213 2149 ENDIF 2214 2150 2215 2216 ! 2217 !-- Interpolation of the w-component: 2218 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 2219 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 2220 kk = rbz(ring,rseg) / dz(1) 2221 ! 2222 !-- Interpolate only if all required information is available on 2223 !-- the current PE: 2224 IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN 2225 IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN 2226 2227 aa = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2228 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2229 bb = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2230 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2231 cc = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2232 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2233 dd = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2234 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2235 gg = dx * dy 2236 2237 w_int_l = ( aa * w(kk,jj,ii) + & 2238 bb * w(kk,jj,ii+1) + & 2239 cc * w(kk,jj+1,ii) + & 2240 dd * w(kk,jj+1,ii+1) & 2241 ) / gg 2242 2243 w_int_u = ( aa * w(kk+1,jj,ii) + & 2244 bb * w(kk+1,jj,ii+1) + & 2245 cc * w(kk+1,jj+1,ii) + & 2246 dd * w(kk+1,jj+1,ii+1) & 2247 ) / gg 2248 2249 w_int_1_l(inot,ring,rseg) = w_int_l + & 2250 ( rbz(ring,rseg) - zw(kk) ) / dz(1) * & 2251 ( w_int_u - w_int_l ) 2252 ELSE 2253 w_int_1_l(inot,ring,rseg) = 0.0_wp 2254 ENDIF 2151 ELSE 2152 v_int_1_l(inot,ring,rseg) = 0.0_wp 2153 ENDIF 2154 2155 2156 ! 2157 !-- Interpolation of the w-component: 2158 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 2159 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 2160 kk = rbz(ring,rseg) / dz(1) 2161 ! 2162 !-- Interpolate only if all required information is available on the current PE: 2163 IF ( ( ii >= nxl ) .AND. ( ii <= nxr ) ) THEN 2164 IF ( ( jj >= nys ) .AND. ( jj <= nyn ) ) THEN 2165 2166 aa = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2167 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2168 bb = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2169 ( ( jj + 1 + 0.5_wp ) * dy - rby(ring,rseg) ) 2170 cc = ( ( ii + 1 + 0.5_wp ) * dx - rbx(ring,rseg) ) * & 2171 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2172 dd = ( rbx(ring,rseg) - ( ii + 0.5_wp ) * dx ) * & 2173 ( rby(ring,rseg) - ( jj + 0.5_wp ) * dy ) 2174 gg = dx * dy 2175 2176 w_int_l = ( aa * w(kk,jj,ii) + bb * w(kk,jj,ii+1) + cc * w(kk,jj+1,ii) & 2177 + dd * w(kk,jj+1,ii+1) ) / gg 2178 2179 w_int_u = ( aa * w(kk+1,jj,ii) + bb * w(kk+1,jj,ii+1) + cc * w(kk+1,jj+1,ii) & 2180 + dd * w(kk+1,jj+1,ii+1) ) / gg 2181 2182 w_int_1_l(inot,ring,rseg) = w_int_l + ( rbz(ring,rseg) - zw(kk) ) / dz(1) * & 2183 ( w_int_u - w_int_l ) 2255 2184 ELSE 2256 2185 w_int_1_l(inot,ring,rseg) = 0.0_wp 2257 2186 ENDIF 2258 2187 2259 ENDDO 2188 ELSE 2189 w_int_1_l(inot,ring,rseg) = 0.0_wp 2190 ENDIF 2191 2260 2192 ENDDO 2261 !$OMP END PARALLEL2262 2193 2263 2194 ENDDO 2264 2265 ! 2266 !-- Exchange between PEs (information required on each PE): 2195 !$OMP END PARALLEL 2196 2197 ENDDO 2198 2199 ! 2200 !-- Exchange between PEs (information required on each PE): 2267 2201 #if defined( __parallel ) 2268 CALL MPI_ALLREDUCE( u_int_1_l, u_int, n_turbines * MAXVAL(nrings) *&2269 2270 CALL MPI_ALLREDUCE( v_int_1_l, v_int, n_turbines * MAXVAL(nrings) *&2271 2272 CALL MPI_ALLREDUCE( w_int_1_l, w_int, n_turbines * MAXVAL(nrings) *&2273 2202 CALL MPI_ALLREDUCE( u_int_1_l, u_int, n_turbines * MAXVAL(nrings) * & 2203 MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) 2204 CALL MPI_ALLREDUCE( v_int_1_l, v_int, n_turbines * MAXVAL(nrings) * & 2205 MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) 2206 CALL MPI_ALLREDUCE( w_int_1_l, w_int, n_turbines * MAXVAL(nrings) * & 2207 MAXVAL(nsegs), MPI_REAL, MPI_SUM, comm2d, ierr ) 2274 2208 #else 2275 2276 2277 2209 u_int = u_int_1_l 2210 v_int = v_int_1_l 2211 w_int = w_int_1_l 2278 2212 #endif 2279 2213 2280 2214 2281 2215 ! 2282 !-- Loop over number of turbines: 2283 2284 DO inot = 1, n_turbines 2216 !-- Loop over number of turbines: 2217 DO inot = 1, n_turbines 2285 2218 pit_loop: DO 2286 2219 2287 IF ( pitch_sw ) THEN 2288 torque_total(inot) = 0.0_wp 2289 thrust_rotor(inot) = 0.0_wp 2290 pitch_angle(inot) = pitch_angle(inot) + 0.25_wp 2291 ! IF ( myid == 0 ) PRINT*, 'Pitch', inot, pitch_angle(inot) 2292 ELSE 2293 cos_yaw = COS(yaw_angle(inot)) 2294 sin_yaw = SIN(yaw_angle(inot)) 2295 IF ( pitch_control ) THEN 2296 pitch_angle(inot) = MAX(pitch_angle_old(inot) - pitch_rate * & 2297 dt_3d , 0.0_wp ) 2220 IF ( pitch_sw ) THEN 2221 torque_total(inot) = 0.0_wp 2222 thrust_rotor(inot) = 0.0_wp 2223 pitch_angle(inot) = pitch_angle(inot) + 0.25_wp 2224 ! IF ( myid == 0 ) PRINT*, 'Pitch', inot, pitch_angle(inot) 2225 ELSE 2226 cos_yaw = COS( yaw_angle(inot) ) 2227 sin_yaw = SIN( yaw_angle(inot) ) 2228 IF ( pitch_control ) THEN 2229 pitch_angle(inot) = MAX( pitch_angle_old(inot) - pitch_rate * dt_3d , 0.0_wp ) 2230 ENDIF 2231 ENDIF 2232 2233 ! 2234 !-- Loop over rings of each turbine: 2235 !$OMP PARALLEL PRIVATE (ring, rseg, sin_rot, cos_rot, re, rea, ren, rote, rota, rotn, & 2236 !$OMP& vtheta, phi_rel, lct, rad_d, alpha_attack, vrel, & 2237 !$OMP& chord, iialpha, iir, turb_cl, tl_factor, thrust_seg, & 2238 !$OMP& torque_seg_y, turb_cd, torque_seg_z, thrust_ring, & 2239 !$OMP& torque_ring_y, torque_ring_z) 2240 !$OMP DO 2241 DO ring = 1, nrings(inot) 2242 ! 2243 !-- Determine distance between each ring (center) and the hub: 2244 cur_r = ( ring - 0.5_wp ) * delta_r(inot) 2245 ! 2246 !-- Loop over segments of each ring of each turbine: 2247 DO rseg = 1, nsegs(ring,inot) 2248 ! 2249 !-- Determine angle of ring segment towards zero degree angle of rotor system 2250 !-- (at zero degree rotor direction vectors aligned with y-axis): 2251 phi_rotor = rseg * 2.0_wp * pi / nsegs(ring,inot) 2252 cos_rot = COS( phi_rotor ) 2253 sin_rot = SIN( phi_rotor ) 2254 ! 2255 !-- Now the direction vectors can be determined with respect to the yaw and tilt angle: 2256 re(1) = cos_rot * sin_yaw 2257 re(2) = cos_rot * cos_yaw 2258 re(3) = sin_rot 2259 2260 ! The current unit vector in azimuthal direction: 2261 rea(1) = - sin_rot * sin_yaw 2262 rea(2) = - sin_rot * cos_yaw 2263 rea(3) = cos_rot 2264 2265 ! 2266 !-- To respect the yawing angle for the calculations of velocities and forces the 2267 !-- unit vectors perpendicular to the rotor area in direction of the positive yaw 2268 !-- angle are defined: 2269 ren(1) = cos_yaw 2270 ren(2) = - sin_yaw 2271 ren(3) = 0.0_wp 2272 ! 2273 !-- Multiplication with the coordinate transformation matrix gives the final unit 2274 !-- vector with consideration of the rotor tilt: 2275 rote = MATMUL( rot_coord_trans(inot,:,:), re ) 2276 rota = MATMUL( rot_coord_trans(inot,:,:), rea ) 2277 rotn = MATMUL( rot_coord_trans(inot,:,:), ren ) 2278 ! 2279 !-- Coordinates of the single segments (center points): 2280 rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1) 2281 2282 rby(ring,rseg) = hub_y(inot) + cur_r * rote(2) 2283 2284 rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3) 2285 2286 ! 2287 !-- !----------------------------------------------------------------------------------! 2288 !-- !-- Calculation of various angles and relative velocities --! 2289 !-- !----------------------------------------------------------------------------------! 2290 ! 2291 !-- In the following the 3D-velocity field is projected its components perpendicular 2292 !-- and parallel to the rotor area. 2293 !-- The calculation of forces will be done in the rotor-coordinates y' and z. 2294 !-- The yaw angle will be reintroduced when the force is applied on the hydrodynamic 2295 !-- equations. 2296 ! 2297 !-- Projection of the xy-velocities relative to the rotor area: 2298 ! 2299 !-- Velocity perpendicular to the rotor area: 2300 u_rot = u_int(inot,ring,rseg) * rotn(1) + v_int(inot,ring,rseg)*rotn(2) + & 2301 w_int(inot,ring,rseg)*rotn(3) 2302 ! 2303 2304 !-- Projection of the 3D-velocity vector in the azimuthal direction: 2305 vtheta(rseg) = rota(1) * u_int(inot,ring,rseg) + rota(2) * v_int(inot,ring,rseg) + & 2306 rota(3) * w_int(inot,ring,rseg) 2307 ! 2308 2309 !-- Determination of the angle phi_rel between the rotor plane and the direction of the 2310 !-- flow relative to the rotor: 2311 phi_rel(rseg) = ATAN2( u_rot , ( rotor_speed(inot) * cur_r - vtheta(rseg) ) ) 2312 2313 ! 2314 !-- Interpolation of the local pitch angle from tabulated values to the current 2315 !-- radial position: 2316 lct = minloc( ABS( cur_r-lrd ) ) 2317 rad_d = cur_r-lrd(lct) 2318 2319 IF ( cur_r == 0.0_wp ) THEN 2320 alpha_attack(rseg) = 0.0_wp 2321 2322 ELSE IF ( cur_r >= lrd(size(ard)) ) THEN 2323 alpha_attack(rseg) = ( ard(size(ard) ) + ard(size(ard) -1 ) ) / 2.0_wp 2324 2325 ELSE 2326 alpha_attack(rseg) = ( ard( lct(1) ) * ( ( lrd( lct(1) + 1 ) - cur_r ) / & 2327 ( lrd( lct(1) + 1 ) - lrd( lct(1) ) ) ) ) + & 2328 ( ard( lct(1) + 1 ) * ( ( cur_r - lrd( lct(1) ) ) / & 2329 ( lrd( lct(1) + 1 ) - lrd( lct(1) ) ) ) ) 2330 ENDIF 2331 2332 ! 2333 !-- In Fortran radian instead of degree is used as unit for all angles. 2334 !-- Therefore, a transformation from angles given in degree to angles given in radian 2335 !-- is necessary here: 2336 alpha_attack(rseg) = alpha_attack(rseg) * ( ( 2.0_wp * pi ) / 360.0_wp ) 2337 ! 2338 !-- Substraction of the local pitch angle to obtain the local angle of attack: 2339 alpha_attack(rseg) = phi_rel(rseg) - alpha_attack(rseg) 2340 ! 2341 !-- Preliminary transformation back from angles given in radian to angles given in 2342 !-- degree: 2343 alpha_attack(rseg) = alpha_attack(rseg) * ( 360.0_wp / ( 2.0_wp * pi ) ) 2344 ! 2345 !-- Correct with collective pitch angle: 2346 alpha_attack(rseg) = alpha_attack(rseg) - pitch_angle(inot) 2347 2348 ! 2349 !-- Determination of the magnitude of the flow velocity relative to the rotor: 2350 vrel(rseg) = SQRT( u_rot**2 + ( rotor_speed(inot) * cur_r - vtheta(rseg) )**2 ) 2351 2352 ! 2353 !-- !----------------------------------------------------------------------------------! 2354 !-- !-- Interpolation of chord as well as lift and drag --! 2355 !-- !-- coefficients from tabulated values --! 2356 !-- !----------------------------------------------------------------------------------! 2357 2358 ! 2359 !-- Interpolation of the chord_length from tabulated values to the current radial 2360 !-- position: 2361 IF ( cur_r == 0.0_wp ) THEN 2362 chord(rseg) = 0.0_wp 2363 2364 ELSE IF ( cur_r >= lrd( size(crd) ) ) THEN 2365 chord(rseg) = ( crd( size(crd) ) + ard( size(crd) - 1 ) ) / 2.0_wp 2366 2367 ELSE 2368 chord(rseg) = ( crd( lct(1) ) * ( ( lrd( lct(1) + 1 ) - cur_r ) / & 2369 ( lrd( lct(1) + 1 ) - lrd( lct(1) ) ) ) ) + ( crd( lct(1) + 1 ) & 2370 * ( ( cur_r - lrd( lct(1) ) ) / ( lrd( lct(1) + 1 ) - & 2371 lrd( lct(1) ) ) ) ) 2372 ENDIF 2373 2374 ! 2375 !-- Determine index of current angle of attack, needed for finding the appropriate 2376 !-- interpolated values of the lift and drag coefficients 2377 !-- (-180.0 degrees = 1, +180.0 degrees = 3601, so one index every 0.1 degrees): 2378 iialpha = CEILING( ( alpha_attack(rseg) + 180.0_wp ) & 2379 * ( 1.0_wp / accu_cl_cd_tab ) ) + 1.0_wp 2380 ! 2381 !-- Determine index of current radial position, needed for finding the appropriate 2382 !-- interpolated values of the lift and drag coefficients (one index every 0.1 m): 2383 iir = CEILING( cur_r * 10.0_wp ) 2384 ! 2385 !-- Read in interpolated values of the lift and drag coefficients for the current 2386 !-- radial position and angle of attack: 2387 turb_cl(rseg) = turb_cl_tab(iialpha,iir) 2388 turb_cd(rseg) = turb_cd_tab(iialpha,iir) 2389 2390 ! 2391 !-- Final transformation back from angles given in degree to angles given in radian: 2392 alpha_attack(rseg) = alpha_attack(rseg) * ( ( 2.0_wp * pi ) / 360.0_wp ) 2393 2394 IF ( tip_loss_correction ) THEN 2395 ! 2396 !-- Tip loss correction following Schito. 2397 !-- Schito applies the tip loss correction only to the lift force. 2398 !-- Therefore, the tip loss correction is only applied to the lift coefficient and 2399 !-- not to the drag coefficient in our case. 2400 IF ( phi_rel(rseg) == 0.0_wp ) THEN 2401 tl_factor = 1.0_wp 2402 2403 ELSE 2404 tl_factor = ( 2.0 / pi ) * & 2405 ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) / & 2406 ( 2.0 * cur_r * ABS( SIN( phi_rel(rseg) ) ) ) ) ) ) 2407 ENDIF 2408 2409 turb_cl(rseg) = tl_factor * turb_cl(rseg) 2410 2411 ENDIF 2412 ! 2413 !-- !----------------------------------------------------------------------------------! 2414 !-- !-- Calculation of the forces --! 2415 !-- !----------------------------------------------------------------------------------! 2416 2417 ! 2418 !-- Calculate the pre_factor for the thrust and torque forces: 2419 pre_factor = 0.5_wp * ( vrel(rseg)**2 ) * 3.0_wp * chord(rseg) * delta_r(inot) & 2420 / nsegs(ring,inot) 2421 2422 ! 2423 !-- Calculate the thrust force (x-component of the total force) for each ring segment: 2424 thrust_seg(rseg) = pre_factor * ( turb_cl(rseg) * COS (phi_rel(rseg) ) + & 2425 turb_cd(rseg) * SIN( phi_rel(rseg) ) ) 2426 2427 ! 2428 !-- Determination of the second of the additional forces acting on the flow in the 2429 !-- azimuthal direction: force vector as basis for torque (torque itself would be the 2430 !-- vector product of the radius vector and the force vector): 2431 torque_seg = pre_factor * ( turb_cl(rseg) * SIN (phi_rel(rseg) ) - & 2432 turb_cd(rseg) * COS( phi_rel(rseg) ) ) 2433 ! 2434 !-- Decomposition of the force vector into two parts: One acting along the 2435 !-- y-direction and one acting along the z-direction of the rotor coordinate system: 2436 torque_seg_y(rseg) = -torque_seg * sin_rot 2437 torque_seg_z(rseg) = torque_seg * cos_rot 2438 2439 ! 2440 !-- Add the segment thrust to the thrust of the whole rotor: 2441 !$OMP CRITICAL 2442 thrust_rotor(inot) = thrust_rotor(inot) + thrust_seg(rseg) 2443 2444 2445 torque_total(inot) = torque_total(inot) + (torque_seg * cur_r) 2446 !$OMP END CRITICAL 2447 2448 ENDDO !-- end of loop over ring segments 2449 2450 ! 2451 !-- Restore the forces into arrays containing all the segments of each ring: 2452 thrust_ring(ring,:) = thrust_seg(:) 2453 torque_ring_y(ring,:) = torque_seg_y(:) 2454 torque_ring_z(ring,:) = torque_seg_z(:) 2455 2456 2457 ENDDO !-- end of loop over rings 2458 !$OMP END PARALLEL 2459 2460 2461 CALL cpu_log( log_point_s(62), 'wtm_controller', 'start' ) 2462 2463 2464 IF ( speed_control ) THEN 2465 ! 2466 !-- Calculation of the current generator speed for rotor speed control: 2467 2468 ! 2469 !-- The acceleration of the rotor speed is calculated from the force balance of the 2470 !-- accelerating torque and the torque of the rotating rotor and generator: 2471 om_rate = ( torque_total(inot) * air_density * gear_efficiency - & 2472 gear_ratio * torque_gen_old(inot) ) / ( rotor_inertia + & 2473 gear_ratio * gear_ratio * generator_inertia ) * dt_3d 2474 2475 ! 2476 !-- The generator speed is given by the product of gear gear_ratio and rotor speed: 2477 generator_speed(inot) = gear_ratio * ( rotor_speed(inot) + om_rate ) 2478 2479 ENDIF 2480 2481 IF ( pitch_control ) THEN 2482 2483 ! 2484 !-- If the current generator speed is above rated, the pitch is not saturated and the 2485 !-- change from the last time step is within the maximum pitch rate, then the pitch loop 2486 !-- is repeated with a pitch gain: 2487 IF ( ( generator_speed(inot) > generator_speed_rated ) .AND. & 2488 ( pitch_angle(inot) < 25.0_wp ) .AND. & 2489 ( pitch_angle(inot) < pitch_angle_old(inot) + pitch_rate * dt_3d ) ) THEN 2490 pitch_sw = .TRUE. 2491 ! 2492 !-- Go back to beginning of pit_loop: 2493 CYCLE pit_loop 2494 ENDIF 2495 2496 ! 2497 !-- The current pitch is saved for the next time step: 2498 pitch_angle_old(inot) = pitch_angle(inot) 2499 pitch_sw = .FALSE. 2500 ENDIF 2501 EXIT pit_loop 2502 ENDDO pit_loop ! Recursive pitch control loop 2503 2504 2505 ! 2506 !-- Call the rotor speed controller: 2507 IF ( speed_control ) THEN 2508 ! 2509 !-- Find processor at i_hub, j_hub: 2510 IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) THEN 2511 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN 2512 CALL wtm_speed_control( inot ) 2298 2513 ENDIF 2299 2514 ENDIF 2300 2515 2301 ! 2302 !-- Loop over rings of each turbine: 2303 !$OMP PARALLEL PRIVATE (ring, rseg, sin_rot, cos_rot, re, rea, ren, rote, rota, rotn, & 2304 !$OMP& vtheta, phi_rel, lct, rad_d, alpha_attack, vrel, & 2305 !$OMP& chord, iialpha, iir, turb_cl, tl_factor, thrust_seg, & 2306 !$OMP& torque_seg_y, turb_cd, torque_seg_z, thrust_ring, & 2307 !$OMP& torque_ring_y, torque_ring_z) 2308 !$OMP DO 2309 DO ring = 1, nrings(inot) 2310 ! 2311 !-- Determine distance between each ring (center) and the hub: 2312 cur_r = (ring - 0.5_wp) * delta_r(inot) 2313 ! 2314 !-- Loop over segments of each ring of each turbine: 2315 DO rseg = 1, nsegs(ring,inot) 2316 ! 2317 !-- Determine angle of ring segment towards zero degree angle of 2318 !-- rotor system (at zero degree rotor direction vectors aligned 2319 !-- with y-axis): 2320 phi_rotor = rseg * 2.0_wp * pi / nsegs(ring,inot) 2321 cos_rot = COS(phi_rotor) 2322 sin_rot = SIN(phi_rotor) 2323 ! 2324 !-- Now the direction vectors can be determined with respect to 2325 !-- the yaw and tilt angle: 2326 re(1) = cos_rot * sin_yaw 2327 re(2) = cos_rot * cos_yaw 2328 re(3) = sin_rot 2329 2330 ! The current unit vector in azimuthal direction: 2331 rea(1) = - sin_rot * sin_yaw 2332 rea(2) = - sin_rot * cos_yaw 2333 rea(3) = cos_rot 2334 2335 ! 2336 !-- To respect the yawing angle for the calculations of 2337 !-- velocities and forces the unit vectors perpendicular to the 2338 !-- rotor area in direction of the positive yaw angle are defined: 2339 ren(1) = cos_yaw 2340 ren(2) = - sin_yaw 2341 ren(3) = 0.0_wp 2342 ! 2343 !-- Multiplication with the coordinate transformation matrix 2344 !-- gives the final unit vector with consideration of the rotor 2345 !-- tilt: 2346 rote = MATMUL( rot_coord_trans(inot,:,:), re ) 2347 rota = MATMUL( rot_coord_trans(inot,:,:), rea ) 2348 rotn = MATMUL( rot_coord_trans(inot,:,:), ren ) 2349 ! 2350 !-- Coordinates of the single segments (center points): 2351 rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1) 2352 2353 rby(ring,rseg) = hub_y(inot) + cur_r * rote(2) 2354 2355 rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3) 2356 2357 ! 2358 !-- !-----------------------------------------------------------! 2359 !-- !-- Calculation of various angles and relative velocities --! 2360 !-- !-----------------------------------------------------------! 2361 ! 2362 !-- In the following the 3D-velocity field is projected its 2363 !-- components perpendicular and parallel to the rotor area 2364 !-- The calculation of forces will be done in the rotor- 2365 !-- coordinates y' and z. 2366 !-- The yaw angle will be reintroduced when the force is applied 2367 !-- on the hydrodynamic equations 2368 ! 2369 !-- Projection of the xy-velocities relative to the rotor area 2370 ! 2371 !-- Velocity perpendicular to the rotor area: 2372 u_rot = u_int(inot,ring,rseg)*rotn(1) + & 2373 v_int(inot,ring,rseg)*rotn(2) + & 2374 w_int(inot,ring,rseg)*rotn(3) 2375 ! 2376 !-- Projection of the 3D-velocity vector in the azimuthal 2377 !-- direction: 2378 vtheta(rseg) = rota(1) * u_int(inot,ring,rseg) + & 2379 rota(2) * v_int(inot,ring,rseg) + & 2380 rota(3) * w_int(inot,ring,rseg) 2381 ! 2382 !-- Determination of the angle phi_rel between the rotor plane 2383 !-- and the direction of the flow relative to the rotor: 2384 2385 phi_rel(rseg) = ATAN2( u_rot , & 2386 ( rotor_speed(inot) * cur_r - & 2387 vtheta(rseg) ) ) 2388 2389 ! 2390 !-- Interpolation of the local pitch angle from tabulated values 2391 !-- to the current radial position: 2392 2393 lct=minloc(ABS(cur_r-lrd)) 2394 rad_d=cur_r-lrd(lct) 2395 2396 IF (cur_r == 0.0_wp) THEN 2397 alpha_attack(rseg) = 0.0_wp 2398 ELSE IF (cur_r >= lrd(size(ard))) THEN 2399 alpha_attack(rseg) = ( ard(size(ard)) + & 2400 ard(size(ard)-1) ) / 2.0_wp 2401 ELSE 2402 alpha_attack(rseg) = ( ard(lct(1)) * & 2403 ( ( lrd(lct(1)+1) - cur_r ) / & 2404 ( lrd(lct(1)+1) - lrd(lct(1)) ) & 2405 ) ) + ( ard(lct(1)+1) * & 2406 ( ( cur_r - lrd(lct(1)) ) / & 2407 ( lrd(lct(1)+1) - lrd(lct(1)) ) ) ) 2408 ENDIF 2409 2410 ! 2411 !-- In Fortran radian instead of degree is used as unit for all 2412 !-- angles. Therefore, a transformation from angles given in 2413 !-- degree to angles given in radian is necessary here: 2414 alpha_attack(rseg) = alpha_attack(rseg) * & 2415 ( (2.0_wp*pi) / 360.0_wp ) 2416 ! 2417 !-- Substraction of the local pitch angle to obtain the local 2418 !-- angle of attack: 2419 alpha_attack(rseg) = phi_rel(rseg) - alpha_attack(rseg) 2420 ! 2421 !-- Preliminary transformation back from angles given in radian 2422 !-- to angles given in degree: 2423 alpha_attack(rseg) = alpha_attack(rseg) * & 2424 ( 360.0_wp / (2.0_wp*pi) ) 2425 ! 2426 !-- Correct with collective pitch angle: 2427 alpha_attack(rseg) = alpha_attack(rseg) - pitch_angle(inot) 2428 2429 ! 2430 !-- Determination of the magnitude of the flow velocity relative 2431 !-- to the rotor: 2432 vrel(rseg) = SQRT( u_rot**2 + & 2433 ( rotor_speed(inot) * cur_r - & 2434 vtheta(rseg) )**2 ) 2435 2436 ! 2437 !-- !-----------------------------------------------------------! 2438 !-- !-- Interpolation of chord as well as lift and drag --! 2439 !-- !-- coefficients from tabulated values --! 2440 !-- !-----------------------------------------------------------! 2441 2442 ! 2443 !-- Interpolation of the chord_length from tabulated values to 2444 !-- the current radial position: 2445 2446 IF (cur_r == 0.0_wp) THEN 2447 chord(rseg) = 0.0_wp 2448 ELSE IF (cur_r >= lrd(size(crd))) THEN 2449 chord(rseg) = (crd(size(crd)) + ard(size(crd)-1)) / 2.0_wp 2450 ELSE 2451 chord(rseg) = ( crd(lct(1)) * & 2452 ( ( lrd(lct(1)+1) - cur_r ) / & 2453 ( lrd(lct(1)+1) - lrd(lct(1)) ) ) ) + & 2454 ( crd(lct(1)+1) * & 2455 ( ( cur_r-lrd(lct(1)) ) / & 2456 ( lrd(lct(1)+1) - lrd(lct(1)) ) ) ) 2457 ENDIF 2458 2459 ! 2460 !-- Determine index of current angle of attack, needed for 2461 !-- finding the appropriate interpolated values of the lift and 2462 !-- drag coefficients (-180.0 degrees = 1, +180.0 degrees = 3601, 2463 !-- so one index every 0.1 degrees): 2464 iialpha = CEILING( ( alpha_attack(rseg) + 180.0_wp ) & 2465 * ( 1.0_wp / accu_cl_cd_tab ) ) + 1.0_wp 2466 ! 2467 !-- Determine index of current radial position, needed for 2468 !-- finding the appropriate interpolated values of the lift and 2469 !-- drag coefficients (one index every 0.1 m): 2470 iir = CEILING( cur_r * 10.0_wp ) 2471 ! 2472 !-- Read in interpolated values of the lift and drag coefficients 2473 !-- for the current radial position and angle of attack: 2474 turb_cl(rseg) = turb_cl_tab(iialpha,iir) 2475 turb_cd(rseg) = turb_cd_tab(iialpha,iir) 2476 2477 ! 2478 !-- Final transformation back from angles given in degree to 2479 !-- angles given in radian: 2480 alpha_attack(rseg) = alpha_attack(rseg) * & 2481 ( (2.0_wp*pi) / 360.0_wp ) 2482 2483 IF ( tip_loss_correction ) THEN 2484 ! 2485 !-- Tip loss correction following Schito 2486 !-- Schito applies the tip loss correction only to the lift force 2487 !-- Therefore, the tip loss correction is only applied to the lift 2488 !-- coefficient and not to the drag coefficient in our case 2516 ENDIF 2517 2518 2519 CALL cpu_log( log_point_s(62), 'wtm_controller', 'stop' ) 2520 2521 CALL cpu_log( log_point_s(63), 'wtm_smearing', 'start' ) 2522 2523 2524 !-- !----------------------------------------------------------------------------------------! 2525 !-- !-- Regularization kernel --! 2526 !-- !-- Smearing of the forces and interpolation to cartesian grid --! 2527 !-- !----------------------------------------------------------------------------------------! 2528 ! 2529 !-- The aerodynamic blade forces need to be distributed smoothly on several mesh points in 2530 !-- order to avoid singular behaviour. 2531 ! 2532 !-- Summation over sum of weighted forces. The weighting factor (calculated in user_init) 2533 !-- includes information on the distance between the center of the grid cell and the rotor 2534 !-- segment under consideration. 2535 ! 2536 !-- To save computing time, apply smearing only for the relevant part of the model domain: 2537 ! 2489 2538 !-- 2490 IF ( phi_rel(rseg) == 0.0_wp ) THEN 2491 tl_factor = 1.0_wp 2492 ELSE 2493 tl_factor = ( 2.0 / pi ) * & 2494 ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) / & 2495 ( 2.0 * cur_r * abs( sin( phi_rel(rseg) ) ) ) ) ) ) 2496 ENDIF 2497 2498 turb_cl(rseg) = tl_factor * turb_cl(rseg) 2499 2500 ENDIF 2501 ! 2502 !-- !-----------------------------------------------------! 2503 !-- !-- Calculation of the forces --! 2504 !-- !-----------------------------------------------------! 2505 2506 ! 2507 !-- Calculate the pre_factor for the thrust and torque forces: 2508 2509 pre_factor = 0.5_wp * (vrel(rseg)**2) * 3.0_wp * & 2510 chord(rseg) * delta_r(inot) / nsegs(ring,inot) 2511 2512 ! 2513 !-- Calculate the thrust force (x-component of the total force) 2514 !-- for each ring segment: 2515 thrust_seg(rseg) = pre_factor * & 2516 ( turb_cl(rseg) * COS(phi_rel(rseg)) + & 2517 turb_cd(rseg) * SIN(phi_rel(rseg)) ) 2518 2519 ! 2520 !-- Determination of the second of the additional forces acting 2521 !-- on the flow in the azimuthal direction: force vector as basis 2522 !-- for torque (torque itself would be the vector product of the 2523 !-- radius vector and the force vector): 2524 torque_seg = pre_factor * & 2525 ( turb_cl(rseg) * SIN(phi_rel(rseg)) - & 2526 turb_cd(rseg) * COS(phi_rel(rseg)) ) 2527 ! 2528 !-- Decomposition of the force vector into two parts: 2529 !-- One acting along the y-direction and one acting along the 2530 !-- z-direction of the rotor coordinate system: 2531 2532 torque_seg_y(rseg) = -torque_seg * sin_rot 2533 torque_seg_z(rseg) = torque_seg * cos_rot 2534 2535 ! 2536 !-- Add the segment thrust to the thrust of the whole rotor 2537 !$OMP CRITICAL 2538 thrust_rotor(inot) = thrust_rotor(inot) + & 2539 thrust_seg(rseg) 2540 2541 2542 torque_total(inot) = torque_total(inot) + (torque_seg * cur_r) 2543 !$OMP END CRITICAL 2544 2545 ENDDO !-- end of loop over ring segments 2546 2547 ! 2548 !-- Restore the forces into arrays containing all the segments of 2549 !-- each ring: 2550 thrust_ring(ring,:) = thrust_seg(:) 2551 torque_ring_y(ring,:) = torque_seg_y(:) 2552 torque_ring_z(ring,:) = torque_seg_z(:) 2553 2554 2555 ENDDO !-- end of loop over rings 2556 !$OMP END PARALLEL 2557 2558 2559 CALL cpu_log( log_point_s(62), 'wtm_controller', 'start' ) 2560 2561 2562 IF ( speed_control ) THEN 2563 ! 2564 !-- Calculation of the current generator speed for rotor speed control 2565 2566 ! 2567 !-- The acceleration of the rotor speed is calculated from 2568 !-- the force balance of the accelerating torque 2569 !-- and the torque of the rotating rotor and generator 2570 om_rate = ( torque_total(inot) * air_density * gear_efficiency - & 2571 gear_ratio * torque_gen_old(inot) ) / & 2572 ( rotor_inertia + & 2573 gear_ratio * gear_ratio * generator_inertia ) * dt_3d 2574 2575 ! 2576 !-- The generator speed is given by the product of gear gear_ratio 2577 !-- and rotor speed 2578 generator_speed(inot) = gear_ratio * ( rotor_speed(inot) + om_rate ) 2579 2539 !-- Calculation of the boundaries: 2540 i_smear(inot) = CEILING( ( rotor_radius(inot) * ABS( roty(inot,1) ) + eps_min ) / dx ) 2541 j_smear(inot) = CEILING( ( rotor_radius(inot) * ABS( roty(inot,2) ) + eps_min ) / dy ) 2542 2543 !$OMP PARALLEL PRIVATE (i, j, k, ring, rseg, flag, dist_u_3d, dist_v_3d, dist_w_3d) 2544 !$OMP DO 2545 DO i = MAX( nxl, i_hub(inot) - i_smear(inot) ), MIN( nxr, i_hub(inot) + i_smear(inot) ) 2546 DO j = MAX( nys, j_hub(inot) - j_smear(inot) ), MIN( nyn, j_hub(inot) + j_smear(inot) ) 2547 ! DO k = MAX( nzb_u_inner(j,i)+1, k_hub(inot) - k_smear(inot) ), & 2548 ! k_hub(inot) + k_smear(inot) 2549 DO k = nzb + 1, k_hub(inot) + k_smear(inot) 2550 DO ring = 1, nrings(inot) 2551 DO rseg = 1, nsegs(ring,inot) 2552 ! 2553 !-- Predetermine flag to mask topography: 2554 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 2555 2556 ! 2557 !-- Determine the square of the distance between the current grid point and 2558 !-- each rotor area segment: 2559 dist_u_3d = ( i * dx - rbx(ring,rseg) )**2 + & 2560 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2561 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2562 2563 dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2564 ( j * dy - rby(ring,rseg) )**2 + & 2565 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2566 2567 dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2568 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2569 ( k * dz(1) - rbz(ring,rseg) )**2 2570 2571 ! 2572 !-- 3D-smearing of the forces with a polynomial function (much faster than 2573 !-- the old Gaussian function), using some parameters that have been 2574 !-- calculated in user_init. The function is only similar to Gaussian 2575 !-- function for squared distances <= eps_min2: 2576 IF ( dist_u_3d <= eps_min2 ) THEN 2577 thrust(k,j,i) = thrust(k,j,i) + thrust_ring(ring,rseg) * & 2578 ( ( pol_a * dist_u_3d - pol_b ) * & 2579 dist_u_3d + 1.0_wp ) * eps_factor * flag 2580 ENDIF 2581 2582 IF ( dist_v_3d <= eps_min2 ) THEN 2583 torque_y(k,j,i) = torque_y(k,j,i) + torque_ring_y(ring,rseg) * & 2584 ( ( pol_a * dist_v_3d - pol_b ) * & 2585 dist_v_3d + 1.0_wp ) * eps_factor * flag 2586 ENDIF 2587 2588 IF ( dist_w_3d <= eps_min2 ) THEN 2589 torque_z(k,j,i) = torque_z(k,j,i) + torque_ring_z(ring,rseg) * & 2590 ( ( pol_a * dist_w_3d - pol_b ) * & 2591 dist_w_3d + 1.0_wp ) * eps_factor * flag 2592 ENDIF 2593 2594 ENDDO ! End of loop over rseg 2595 ENDDO ! End of loop over ring 2596 2597 ! 2598 !-- Rotation of force components: 2599 rot_tend_x(k,j,i) = rot_tend_x(k,j,i) + ( thrust(k,j,i) * rotx(inot,1) + & 2600 torque_y(k,j,i) * roty(inot,1) + torque_z(k,j,i) * & 2601 rotz(inot,1) ) * flag 2602 2603 rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + ( thrust(k,j,i) * rotx(inot,2) + & 2604 torque_y(k,j,i) * roty(inot,2) + torque_z(k,j,i) * & 2605 rotz(inot,2) ) * flag 2606 2607 rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + ( thrust(k,j,i) * rotx(inot,3) + & 2608 torque_y(k,j,i) * roty(inot,3) + torque_z(k,j,i) * & 2609 rotz(inot,3) ) * flag 2610 2611 ENDDO ! End of loop over k 2612 ENDDO ! End of loop over j 2613 ENDDO ! End of loop over i 2614 !$OMP END PARALLEL 2615 2616 CALL cpu_log( log_point_s(63), 'wtm_smearing', 'stop' ) 2617 2618 ENDDO !-- end of loop over turbines 2619 2620 2621 IF ( yaw_control ) THEN 2622 ! 2623 !-- Allocate arrays for yaw control at first call. Can't be allocated before dt_3d is set. 2624 IF ( start_up ) THEN 2625 wdlon= MAX( 1 , NINT( 30.0_wp / dt_3d ) ) ! 30s running mean array 2626 ALLOCATE( wd30(1:n_turbines,1:WDLON) ) 2627 wd30 = 999.0_wp ! Set to dummy value 2628 ALLOCATE( wd30_l(1:WDLON) ) 2629 2630 wdsho = MAX( 1 , NINT( 2.0_wp / dt_3d ) ) ! 2s running mean array 2631 ALLOCATE( wd2(1:n_turbines,1:wdsho) ) 2632 wd2 = 999.0_wp ! Set to dummy value 2633 ALLOCATE( wd2_l(1:wdsho) ) 2634 start_up = .FALSE. 2635 ENDIF 2636 2637 ! 2638 !-- Calculate the inflow wind speed: 2639 !-- 2640 !-- Loop over number of turbines: 2641 DO inot = 1, n_turbines 2642 ! 2643 !-- Find processor at i_hub, j_hub: 2644 IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) THEN 2645 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN 2646 u_inflow_l(inot) = u(k_hub(inot),j_hub(inot),i_hub(inot)) 2647 wdir_l(inot) = -1.0_wp * ATAN2( 0.5_wp * & 2648 ( v(k_hub(inot), j_hub(inot), i_hub(inot) + 1) + & 2649 v(k_hub(inot), j_hub(inot), i_hub(inot)) ) , 0.5_wp * & 2650 ( u(k_hub(inot), j_hub(inot) + 1, i_hub(inot)) + & 2651 u(k_hub(inot), j_hub(inot), i_hub(inot)) ) ) 2652 2653 CALL wtm_yawcontrol( inot ) 2654 2655 yaw_angle_l(inot) = yaw_angle(inot) 2656 2657 ENDIF 2580 2658 ENDIF 2581 2659 2582 IF ( pitch_control ) THEN 2583 2584 ! 2585 !-- If the current generator speed is above rated, the pitch is not 2586 !-- saturated and the change from the last time step is within the 2587 !-- maximum pitch rate, then the pitch loop is repeated with a pitch 2588 !-- gain 2589 IF ( ( generator_speed(inot) > generator_speed_rated ) .AND. & 2590 ( pitch_angle(inot) < 25.0_wp ) .AND. & 2591 ( pitch_angle(inot) < pitch_angle_old(inot) + & 2592 pitch_rate * dt_3d ) ) THEN 2593 pitch_sw = .TRUE. 2594 ! 2595 !-- Go back to beginning of pit_loop 2596 CYCLE pit_loop 2597 ENDIF 2598 2599 ! 2600 !-- The current pitch is saved for the next time step 2601 pitch_angle_old(inot) = pitch_angle(inot) 2602 pitch_sw = .FALSE. 2603 ENDIF 2604 EXIT pit_loop 2605 ENDDO pit_loop ! Recursive pitch control loop 2606 2607 2608 ! 2609 !-- Call the rotor speed controller 2610 2611 IF ( speed_control ) THEN 2612 ! 2613 !-- Find processor at i_hub, j_hub 2614 IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) & 2615 THEN 2616 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) )& 2617 THEN 2618 CALL wtm_speed_control( inot ) 2619 ENDIF 2620 ENDIF 2621 2622 ENDIF 2623 2624 2625 CALL cpu_log( log_point_s(62), 'wtm_controller', 'stop' ) 2626 2627 CALL cpu_log( log_point_s(63), 'wtm_smearing', 'start' ) 2628 2629 2630 !-- !-----------------------------------------------------------------! 2631 !-- !-- Regularization kernel --! 2632 !-- !-- Smearing of the forces and interpolation to cartesian grid --! 2633 !-- !-----------------------------------------------------------------! 2634 ! 2635 !-- The aerodynamic blade forces need to be distributed smoothly on 2636 !-- several mesh points in order to avoid singular behaviour 2637 ! 2638 !-- Summation over sum of weighted forces. The weighting factor 2639 !-- (calculated in user_init) includes information on the distance 2640 !-- between the center of the grid cell and the rotor segment under 2641 !-- consideration 2642 ! 2643 !-- To save computing time, apply smearing only for the relevant part 2644 !-- of the model domain: 2645 ! 2646 !-- 2647 !-- Calculation of the boundaries: 2648 i_smear(inot) = CEILING( ( rotor_radius(inot) * ABS( roty(inot,1) ) + & 2649 eps_min ) / dx ) 2650 j_smear(inot) = CEILING( ( rotor_radius(inot) * ABS( roty(inot,2) ) + & 2651 eps_min ) / dy ) 2652 2653 !$OMP PARALLEL PRIVATE (i, j, k, ring, rseg, flag, dist_u_3d, dist_v_3d, dist_w_3d) 2654 !$OMP DO 2655 DO i = MAX( nxl, i_hub(inot) - i_smear(inot) ), & 2656 MIN( nxr, i_hub(inot) + i_smear(inot) ) 2657 DO j = MAX( nys, j_hub(inot) - j_smear(inot) ), & 2658 MIN( nyn, j_hub(inot) + j_smear(inot) ) 2659 ! DO k = MAX( nzb_u_inner(j,i)+1, k_hub(inot) - k_smear(inot) ), & 2660 ! k_hub(inot) + k_smear(inot) 2661 DO k = nzb+1, k_hub(inot) + k_smear(inot) 2662 DO ring = 1, nrings(inot) 2663 DO rseg = 1, nsegs(ring,inot) 2664 ! 2665 !-- Predetermine flag to mask topography 2666 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 2667 2668 ! 2669 !-- Determine the square of the distance between the 2670 !-- current grid point and each rotor area segment: 2671 dist_u_3d = ( i * dx - rbx(ring,rseg) )**2 + & 2672 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2673 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2674 dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2675 ( j * dy - rby(ring,rseg) )**2 + & 2676 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2677 dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2678 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2679 ( k * dz(1) - rbz(ring,rseg) )**2 2680 2681 ! 2682 !-- 3D-smearing of the forces with a polynomial function 2683 !-- (much faster than the old Gaussian function), using 2684 !-- some parameters that have been calculated in user_init. 2685 !-- The function is only similar to Gaussian function for 2686 !-- squared distances <= eps_min2: 2687 IF ( dist_u_3d <= eps_min2 ) THEN 2688 thrust(k,j,i) = thrust(k,j,i) + & 2689 thrust_ring(ring,rseg) * & 2690 ( ( pol_a * dist_u_3d - pol_b ) * & 2691 dist_u_3d + 1.0_wp ) * eps_factor *& 2692 flag 2693 ENDIF 2694 IF ( dist_v_3d <= eps_min2 ) THEN 2695 torque_y(k,j,i) = torque_y(k,j,i) + & 2696 torque_ring_y(ring,rseg) * & 2697 ( ( pol_a * dist_v_3d - pol_b ) * & 2698 dist_v_3d + 1.0_wp ) * eps_factor * & 2699 flag 2700 ENDIF 2701 IF ( dist_w_3d <= eps_min2 ) THEN 2702 torque_z(k,j,i) = torque_z(k,j,i) + & 2703 torque_ring_z(ring,rseg) * & 2704 ( ( pol_a * dist_w_3d - pol_b ) * & 2705 dist_w_3d + 1.0_wp ) * eps_factor * & 2706 flag 2707 ENDIF 2708 2709 ENDDO ! End of loop over rseg 2710 ENDDO ! End of loop over ring 2711 2712 ! 2713 !-- Rotation of force components: 2714 rot_tend_x(k,j,i) = rot_tend_x(k,j,i) + ( & 2715 thrust(k,j,i)*rotx(inot,1) + & 2716 torque_y(k,j,i)*roty(inot,1) + & 2717 torque_z(k,j,i)*rotz(inot,1) & 2718 ) * flag 2719 2720 rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + ( & 2721 thrust(k,j,i)*rotx(inot,2) + & 2722 torque_y(k,j,i)*roty(inot,2) + & 2723 torque_z(k,j,i)*rotz(inot,2) & 2724 ) * flag 2725 2726 rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + ( & 2727 thrust(k,j,i)*rotx(inot,3) + & 2728 torque_y(k,j,i)*roty(inot,3) + & 2729 torque_z(k,j,i)*rotz(inot,3) & 2730 ) * flag 2731 2732 ENDDO ! End of loop over k 2733 ENDDO ! End of loop over j 2734 ENDDO ! End of loop over i 2735 !$OMP END PARALLEL 2736 2737 CALL cpu_log( log_point_s(63), 'wtm_smearing', 'stop' ) 2738 2739 ENDDO !-- end of loop over turbines 2740 2741 2742 IF ( yaw_control ) THEN 2743 ! 2744 !-- Allocate arrays for yaw control at first call 2745 !-- Can't be allocated before dt_3d is set 2746 IF ( start_up ) THEN 2747 WDLON = MAX( 1 , NINT( 30.0_wp / dt_3d ) ) ! 30s running mean array 2748 ALLOCATE( wd30(1:n_turbines,1:WDLON) ) 2749 wd30 = 999.0_wp ! Set to dummy value 2750 ALLOCATE( wd30_l(1:WDLON) ) 2751 2752 WDSHO = MAX( 1 , NINT( 2.0_wp / dt_3d ) ) ! 2s running mean array 2753 ALLOCATE( wd2(1:n_turbines,1:WDSHO) ) 2754 wd2 = 999.0_wp ! Set to dummy value 2755 ALLOCATE( wd2_l(1:WDSHO) ) 2756 start_up = .FALSE. 2757 ENDIF 2758 2759 ! 2760 !-- Calculate the inflow wind speed 2761 !-- 2762 !-- Loop over number of turbines: 2763 DO inot = 1, n_turbines 2764 ! 2765 !-- Find processor at i_hub, j_hub 2766 IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) & 2767 THEN 2768 IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) )& 2769 THEN 2770 2771 u_inflow_l(inot) = u(k_hub(inot),j_hub(inot),i_hub(inot)) 2772 2773 wdir_l(inot) = -1.0_wp * ATAN2( & 2774 0.5_wp * ( v(k_hub(inot),j_hub(inot),i_hub(inot)+1) + & 2775 v(k_hub(inot),j_hub(inot),i_hub(inot)) ) , & 2776 0.5_wp * ( u(k_hub(inot),j_hub(inot)+1,i_hub(inot)) + & 2777 u(k_hub(inot),j_hub(inot),i_hub(inot)) ) ) 2778 2779 CALL wtm_yawcontrol( inot ) 2780 2781 yaw_angle_l(inot) = yaw_angle(inot) 2782 2783 ENDIF 2784 ENDIF 2785 2786 ENDDO !-- end of loop over turbines 2787 2788 ! 2789 !-- Transfer of information to the other cpus 2790 #if defined( __parallel ) 2791 CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, n_turbines, MPI_REAL, & 2792 MPI_SUM, comm2d, ierr ) 2793 CALL MPI_ALLREDUCE( wdir_l, wdir, n_turbines, MPI_REAL, MPI_SUM, & 2794 comm2d, ierr ) 2795 CALL MPI_ALLREDUCE( yaw_angle_l, yaw_angle, n_turbines, MPI_REAL, & 2796 MPI_SUM, comm2d, ierr ) 2660 ENDDO ! end of loop over turbines 2661 2662 ! 2663 !-- Transfer of information to the other cpus: 2664 #if defined( __parallel ) 2665 CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, n_turbines, MPI_REAL, & 2666 MPI_SUM, comm2d, ierr ) 2667 CALL MPI_ALLREDUCE( wdir_l, wdir, n_turbines, MPI_REAL, MPI_SUM, & 2668 comm2d, ierr ) 2669 CALL MPI_ALLREDUCE( yaw_angle_l, yaw_angle, n_turbines, MPI_REAL, & 2670 MPI_SUM, comm2d, ierr ) 2797 2671 #else 2798 u_inflow= u_inflow_l2799 wdir= wdir_l2800 2672 u_inflow = u_inflow_l 2673 wdir = wdir_l 2674 yaw_angle = yaw_angle_l 2801 2675 2802 2676 #endif 2803 DOinot = 1, n_turbines2804 ! 2805 !-- Update rotor orientation2806 2807 2808 2809 2810 ENDIF ! end of yaw control2811 2812 2813 ! 2814 !-- Transfer of information to the other cpus2815 ! CALL MPI_ALLREDUCE( generator_speed, generator_speed_old, n_turbines,&2816 ! 2817 #if defined( __parallel ) 2818 CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, n_turbines,&2819 2820 CALL MPI_ALLREDUCE( rotor_speed_l, rotor_speed, n_turbines,&2821 2822 CALL MPI_ALLREDUCE( generator_speed_f, generator_speed_f_old, n_turbines,&2823 2677 DO inot = 1, n_turbines 2678 ! 2679 !-- Update rotor orientation: 2680 CALL wtm_rotate_rotor( inot ) 2681 2682 ENDDO ! End of loop over turbines 2683 2684 ENDIF ! end of yaw control 2685 2686 IF ( speed_control ) THEN 2687 ! 2688 !-- Transfer of information to the other cpus: 2689 ! CALL MPI_ALLREDUCE( generator_speed, generator_speed_old, n_turbines, & 2690 ! MPI_REAL,MPI_SUM, comm2d, ierr ) 2691 #if defined( __parallel ) 2692 CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, n_turbines, & 2693 MPI_REAL, MPI_SUM, comm2d, ierr ) 2694 CALL MPI_ALLREDUCE( rotor_speed_l, rotor_speed, n_turbines, & 2695 MPI_REAL, MPI_SUM, comm2d, ierr ) 2696 CALL MPI_ALLREDUCE( generator_speed_f, generator_speed_f_old, n_turbines, & 2697 MPI_REAL, MPI_SUM, comm2d, ierr ) 2824 2698 #else 2825 torque_gen_old= torque_gen2826 rotor_speed= rotor_speed_l2827 2699 torque_gen_old = torque_gen 2700 rotor_speed = rotor_speed_l 2701 generator_speed_f_old = generator_speed_f 2828 2702 #endif 2829 2703 2830 ENDIF2831 2832 2704 ENDIF 2833 2705 2834 CALL cpu_log( log_point_s(61), 'wtm_forces', 'stop' ) 2835 2836 2837 END SUBROUTINE wtm_forces 2838 2839 2840 !------------------------------------------------------------------------------! 2706 ENDIF 2707 2708 CALL cpu_log( log_point_s(61), 'wtm_forces', 'stop' ) 2709 2710 2711 END SUBROUTINE wtm_forces 2712 2713 2714 !--------------------------------------------------------------------------------------------------! 2841 2715 ! Description: 2842 2716 ! ------------ 2843 2717 !> Yaw controller for the wind turbine model 2844 !------------------------------------------------------------------------------! 2845 SUBROUTINE wtm_yawcontrol( inot ) 2846 2847 USE kinds 2848 2849 IMPLICIT NONE 2850 2851 INTEGER(iwp) :: inot 2852 INTEGER(iwp) :: i_wd_30 2853 REAL(wp) :: missal 2854 2855 i_wd_30 = 0_iwp 2856 2857 ! 2858 !-- The yaw controller computes a 30s running mean of the wind direction. 2859 !-- If the difference between turbine alignment and wind direction exceeds 2860 !-- 5 degrees, the turbine is yawed. The mechanism stops as soon as the 2s-running 2861 !-- mean of the missalignment is smaller than 0.5 degrees. 2862 !-- Attention: If the timestep during the simulation changes significantly 2863 !-- the lengths of the running means change and it does not correspond to 2864 !-- 30s/2s anymore. 2865 !-- ! Needs to be modified for these situations ! 2866 !-- For wind from the east, the averaging of the wind direction could cause 2867 !-- problems and the yaw controller is probably flawed. -> Routine for 2868 !-- averaging needs to be improved! 2869 ! 2870 !-- Check if turbine is not yawing 2871 IF ( .NOT. doyaw(inot) ) THEN 2872 ! 2873 !-- Write current wind direction into array 2874 wd30_l = wd30(inot,:) 2875 wd30_l = CSHIFT( wd30_l, SHIFT=-1 ) 2876 wd30_l(1) = wdir(inot) 2877 ! 2878 !-- Check if array is full ( no more dummies ) 2879 IF ( .NOT. ANY( wd30_l == 999.) ) THEN 2880 2881 missal = SUM( wd30_l ) / SIZE( wd30_l ) - yaw_angle(inot) 2882 ! 2883 !-- Check if turbine is missaligned by more than yaw_misalignment_max 2884 IF ( ABS( missal ) > yaw_misalignment_max ) THEN 2885 ! 2886 !-- Check in which direction to yaw 2887 yawdir(inot) = SIGN( 1.0_wp, missal ) 2888 ! 2889 !-- Start yawing of turbine 2890 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2891 doyaw(inot) = .TRUE. 2892 wd30_l = 999. ! fill with dummies again 2893 ENDIF 2718 !--------------------------------------------------------------------------------------------------! 2719 SUBROUTINE wtm_yawcontrol( inot ) 2720 2721 USE kinds 2722 2723 IMPLICIT NONE 2724 2725 INTEGER(iwp) :: inot 2726 INTEGER(iwp) :: i_wd_30 2727 2728 REAL(wp) :: missal 2729 2730 i_wd_30 = 0_iwp 2731 2732 ! 2733 !-- The yaw controller computes a 30s running mean of the wind direction. If the difference 2734 !-- between turbine alignment and wind direction exceeds 5 degrees, the turbine is yawed. The 2735 !-- mechanism stops as soon as the 2s-running mean of the missalignment is smaller than 0.5 2736 !-- degrees. Attention: If the timestep during the simulation changes significantly the lengths 2737 !-- of the running means change and it does not correspond to 30s/2s anymore. 2738 !-- ! Needs to be modified for these situations ! 2739 !-- For wind from the east, the averaging of the wind direction could cause problems and the yaw 2740 !-- controller is probably flawed. -> Routine for averaging needs to be improved! 2741 ! 2742 !-- Check if turbine is not yawing: 2743 IF ( .NOT. doyaw(inot) ) THEN 2744 ! 2745 !-- Write current wind direction into array: 2746 wd30_l = wd30(inot,:) 2747 wd30_l = CSHIFT( wd30_l, SHIFT = -1 ) 2748 wd30_l(1) = wdir(inot) 2749 ! 2750 !-- Check if array is full ( no more dummies ): 2751 IF ( .NOT. ANY( wd30_l == 999.) ) THEN 2752 2753 missal = SUM( wd30_l ) / SIZE( wd30_l ) - yaw_angle(inot) 2754 ! 2755 !-- Check if turbine is missaligned by more than yaw_misalignment_max: 2756 IF ( ABS( missal ) > yaw_misalignment_max ) THEN 2757 ! 2758 !-- Check in which direction to yaw: 2759 yawdir(inot) = SIGN( 1.0_wp, missal ) 2760 ! 2761 !-- Start yawing of turbine: 2762 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2763 doyaw(inot) = .TRUE. 2764 wd30_l = 999. ! fill with dummies again 2894 2765 ENDIF 2895 2896 wd30(inot,:) = wd30_l 2897 2898 ! 2899 !-- If turbine is already yawing: 2900 !-- Initialize 2 s running mean and yaw until the missalignment is smaller 2901 !-- than yaw_misalignment_min 2902 2766 ENDIF 2767 2768 wd30(inot,:) = wd30_l 2769 2770 ! 2771 !-- If turbine is already yawing: 2772 !-- Initialize 2 s running mean and yaw until the missalignment is smaller than 2773 !-- yaw_misalignment_min 2774 2775 ELSE 2776 ! 2777 !-- Initialize 2 s running mean: 2778 2779 wd2_l = wd2(inot,:) 2780 wd2_l = CSHIFT( wd2_l, SHIFT = -1 ) 2781 wd2_l(1) = wdir(inot) 2782 ! 2783 !-- Check if array is full ( no more dummies ): 2784 IF ( .NOT. ANY( wd2_l == 999.0_wp ) ) THEN 2785 ! 2786 !-- Calculate missalignment of turbine: 2787 missal = SUM( wd2_l - yaw_angle(inot) ) / SIZE( wd2_l ) 2788 ! 2789 !-- Check if missalignment is still larger than 0.5 degree and if the yaw direction is 2790 !-- still right: 2791 IF ( ( ABS( missal ) > yaw_misalignment_min ) .AND. & 2792 ( yawdir(inot) == SIGN( 1.0_wp, missal ) ) ) THEN 2793 ! 2794 !-- Continue yawing: 2795 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2796 2797 ELSE 2798 ! 2799 !-- Stop yawing: 2800 doyaw(inot) = .FALSE. 2801 wd2_l = 999.0_wp ! fill with dummies again 2802 ENDIF 2903 2803 ELSE 2904 2804 ! 2905 !-- Initialize 2 s running mean 2906 wd2_l = wd2(inot,:) 2907 wd2_l = CSHIFT( wd2_l, SHIFT = -1 ) 2908 wd2_l(1) = wdir(inot) 2909 ! 2910 !-- Check if array is full ( no more dummies ) 2911 IF ( .NOT. ANY( wd2_l == 999.0_wp ) ) THEN 2912 ! 2913 !-- Calculate missalignment of turbine 2914 missal = SUM( wd2_l - yaw_angle(inot) ) / SIZE( wd2_l ) 2915 ! 2916 !-- Check if missalignment is still larger than 0.5 degree and if the 2917 !-- yaw direction is still right 2918 IF ( ( ABS( missal ) > yaw_misalignment_min ) .AND. & 2919 ( yawdir(inot) == SIGN( 1.0_wp, missal ) ) ) THEN 2920 ! 2921 !-- Continue yawing 2922 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2923 ELSE 2924 ! 2925 !-- Stop yawing 2926 doyaw(inot) = .FALSE. 2927 wd2_l = 999.0_wp ! fill with dummies again 2928 ENDIF 2929 ELSE 2930 ! 2931 !-- Continue yawing 2932 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2933 ENDIF 2934 2935 wd2(inot,:) = wd2_l 2936 2805 !-- Continue yawing: 2806 yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d 2937 2807 ENDIF 2938 2808 2939 END SUBROUTINE wtm_yawcontrol 2940 2941 2942 !------------------------------------------------------------------------------! 2809 wd2(inot,:) = wd2_l 2810 2811 ENDIF 2812 2813 END SUBROUTINE wtm_yawcontrol 2814 2815 2816 !--------------------------------------------------------------------------------------------------! 2943 2817 ! Description: 2944 2818 ! ------------ 2945 2819 !> Initialization of the speed control 2946 !------------------------------------------------------------------------------! 2947 SUBROUTINE wtm_init_speed_control 2948 2949 2950 IMPLICIT NONE 2951 2952 ! 2953 !-- If speed control is set, remaining variables and control_parameters for 2954 !-- the control algorithm are calculated 2955 ! 2956 !-- Calculate slope constant for region 15 2957 slope15 = ( region_2_slope * region_2_min * region_2_min ) / & 2958 ( region_2_min - region_15_min ) 2959 ! 2960 !-- Calculate upper limit of slipage region 2961 vs_sysp = generator_speed_rated / 1.1_wp 2962 ! 2963 !-- Calculate slope of slipage region 2964 region_25_slope = ( generator_power_rated / generator_speed_rated ) / & 2965 ( generator_speed_rated - vs_sysp ) 2966 ! 2967 !-- Calculate lower limit of slipage region 2968 region_25_min = ( region_25_slope - SQRT( region_25_slope * & 2969 ( region_25_slope - 4.0_wp * region_2_slope * vs_sysp ) ) ) / & 2970 ( 2.0_wp * region_2_slope ) 2971 ! 2972 !-- Frequency for the simple low pass filter 2973 Fcorner = 0.25_wp 2974 ! 2975 !-- At the first timestep the torque is set to its maximum to prevent 2976 !-- an overspeeding of the rotor 2977 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 2978 torque_gen_old(:) = generator_torque_max 2979 ENDIF 2980 2981 END SUBROUTINE wtm_init_speed_control 2982 2983 2984 !------------------------------------------------------------------------------! 2820 !--------------------------------------------------------------------------------------------------! 2821 SUBROUTINE wtm_init_speed_control 2822 2823 2824 IMPLICIT NONE 2825 2826 ! 2827 !-- If speed control is set, remaining variables and control_parameters for the control algorithm 2828 !-- are calculated: 2829 ! 2830 !-- Calculate slope constant for region 15: 2831 slope15 = ( region_2_slope * region_2_min * region_2_min ) / & 2832 ( region_2_min - region_15_min ) 2833 ! 2834 !-- Calculate upper limit of slipage region: 2835 vs_sysp = generator_speed_rated / 1.1_wp 2836 ! 2837 !-- Calculate slope of slipage region: 2838 region_25_slope = ( generator_power_rated / generator_speed_rated ) / & 2839 ( generator_speed_rated - vs_sysp ) 2840 ! 2841 !-- Calculate lower limit of slipage region: 2842 region_25_min = ( region_25_slope - SQRT( region_25_slope * & 2843 ( region_25_slope - 4.0_wp * region_2_slope * vs_sysp ) ) ) / & 2844 ( 2.0_wp * region_2_slope ) 2845 ! 2846 !-- Frequency for the simple low pass filter: 2847 fcorner = 0.25_wp 2848 ! 2849 !-- At the first timestep the torque is set to its maximum to prevent an overspeeding of the rotor: 2850 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 2851 torque_gen_old(:) = generator_torque_max 2852 ENDIF 2853 2854 END SUBROUTINE wtm_init_speed_control 2855 2856 2857 !--------------------------------------------------------------------------------------------------! 2985 2858 ! Description: 2986 2859 ! ------------ 2987 2860 !> Simple controller for the regulation of the rotor speed 2988 !------------------------------------------------------------------------------! 2989 SUBROUTINE wtm_speed_control( inot ) 2990 2991 2992 IMPLICIT NONE 2993 2994 INTEGER(iwp) :: inot 2995 2996 2997 ! 2998 !-- The controller is based on the fortran script from Jonkman 2999 !-- et al. 2009 "Definition of a 5 MW Reference Wind Turbine for 3000 !-- offshore system developement" 3001 3002 ! 3003 !-- The generator speed is filtered by a low pass filter 3004 !-- for the control of the generator torque 3005 lp_coeff = EXP( -2.0_wp * 3.14_wp * dt_3d * Fcorner ) 3006 generator_speed_f(inot) = ( 1.0_wp - lp_coeff ) * generator_speed(inot) + lp_coeff *& 3007 generator_speed_f_old(inot) 3008 3009 IF ( generator_speed_f(inot) <= region_15_min ) THEN 3010 ! 3011 !-- Region 1: Generator torque is set to zero to accelerate the rotor: 3012 torque_gen(inot) = 0 3013 3014 ELSEIF ( generator_speed_f(inot) <= region_2_min ) THEN 3015 ! 3016 !-- Region 1.5: Generator torque is increasing linearly with rotor speed: 3017 torque_gen(inot) = slope15 * ( generator_speed_f(inot) - region_15_min ) 3018 3019 ELSEIF ( generator_speed_f(inot) <= region_25_min ) THEN 3020 ! 3021 !-- Region 2: Generator torque is increased by the square of the generator 3022 !-- speed to keep the TSR optimal: 3023 torque_gen(inot) = region_2_slope * generator_speed_f(inot) * generator_speed_f(inot) 3024 3025 ELSEIF ( generator_speed_f(inot) < generator_speed_rated ) THEN 3026 ! 3027 !-- Region 2.5: Slipage region between 2 and 3: 3028 torque_gen(inot) = region_25_slope * ( generator_speed_f(inot) - vs_sysp ) 3029 3030 ELSE 3031 ! 3032 !-- Region 3: Generator torque is antiproportional to the rotor speed to 3033 !-- keep the power constant: 3034 torque_gen(inot) = generator_power_rated / generator_speed_f(inot) 3035 3036 ENDIF 3037 ! 3038 !-- Calculate torque rate and confine with a max 3039 trq_rate = ( torque_gen(inot) - torque_gen_old(inot) ) / dt_3d 3040 trq_rate = MIN( MAX( trq_rate, -1.0_wp * generator_torque_rate_max ), & 3041 generator_torque_rate_max ) 3042 ! 3043 !-- Calculate new gen torque and confine with max torque 3044 torque_gen(inot) = torque_gen_old(inot) + trq_rate * dt_3d 3045 torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max ) 3046 ! 3047 !-- Overwrite values for next timestep 3048 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 3049 3050 3051 END SUBROUTINE wtm_speed_control 3052 3053 3054 !------------------------------------------------------------------------------! 2861 !--------------------------------------------------------------------------------------------------! 2862 SUBROUTINE wtm_speed_control( inot ) 2863 2864 2865 IMPLICIT NONE 2866 2867 INTEGER(iwp):: inot 2868 2869 2870 ! 2871 !-- The controller is based on the fortran script from Jonkman et al. 2009 "Definition of a 5 MW 2872 !-- Reference Wind Turbine for offshore system developement" 2873 2874 ! 2875 !-- The generator speed is filtered by a low pass filter for the control of the generator torque: 2876 lp_coeff = EXP( -2.0_wp * 3.14_wp * dt_3d * fcorner) 2877 generator_speed_f(inot) = ( 1.0_wp - lp_coeff ) * generator_speed(inot) + lp_coeff * & 2878 generator_speed_f_old(inot) 2879 2880 IF ( generator_speed_f(inot) <= region_15_min ) THEN 2881 ! 2882 !-- Region 1: Generator torque is set to zero to accelerate the rotor: 2883 torque_gen(inot) = 0 2884 2885 ELSEIF ( generator_speed_f(inot) <= region_2_min ) THEN 2886 ! 2887 !-- Region 1.5: Generator torque is increasing linearly with rotor speed: 2888 torque_gen(inot) = slope15 * ( generator_speed_f(inot) - region_15_min ) 2889 2890 ELSEIF ( generator_speed_f(inot) <= region_25_min ) THEN 2891 ! 2892 !-- Region 2: Generator torque is increased by the square of the generator speed to keep the 2893 !-- TSR optimal: 2894 torque_gen(inot) = region_2_slope * generator_speed_f(inot) * generator_speed_f(inot) 2895 2896 ELSEIF ( generator_speed_f(inot) < generator_speed_rated ) THEN 2897 ! 2898 !-- Region 2.5: Slipage region between 2 and 3: 2899 torque_gen(inot) = region_25_slope * ( generator_speed_f(inot) - vs_sysp ) 2900 2901 ELSE 2902 ! 2903 !-- Region 3: Generator torque is antiproportional to the rotor speed to keep the power 2904 !-- constant: 2905 torque_gen(inot) = generator_power_rated / generator_speed_f(inot) 2906 2907 ENDIF 2908 ! 2909 !-- Calculate torque rate and confine with a max: 2910 trq_rate = ( torque_gen(inot) - torque_gen_old(inot) ) / dt_3d 2911 trq_rate = MIN( MAX( trq_rate, -1.0_wp * generator_torque_rate_max ), generator_torque_rate_max ) 2912 ! 2913 !-- Calculate new gen torque and confine with max torque: 2914 torque_gen(inot) = torque_gen_old(inot) + trq_rate * dt_3d 2915 torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max ) 2916 ! 2917 !-- Overwrite values for next timestep: 2918 rotor_speed_l(inot) = generator_speed(inot) / gear_ratio 2919 2920 2921 END SUBROUTINE wtm_speed_control 2922 2923 2924 !--------------------------------------------------------------------------------------------------! 3055 2925 ! Description: 3056 2926 ! ------------ 3057 !> Application of the additional forces generated by the wind turbine on the 3058 !> flow components(tendency terms)2927 !> Application of the additional forces generated by the wind turbine on the flow components 2928 !> (tendency terms) 3059 2929 !> Call for all grid points 3060 !------------------------------------------------------------------------------! 3061 SUBROUTINE wtm_actions( location ) 3062 3063 3064 CHARACTER (LEN=*) :: location !< 3065 3066 INTEGER(iwp) :: i !< running index 3067 INTEGER(iwp) :: j !< running index 3068 INTEGER(iwp) :: k !< running index 3069 3070 3071 SELECT CASE ( location ) 3072 3073 CASE ( 'before_timestep' ) 3074 3075 CALL wtm_forces 3076 3077 CASE ( 'u-tendency' ) 3078 ! 3079 !-- Apply the x-component of the force to the u-component of the flow: 3080 IF ( time_since_reference_point >= time_turbine_on ) THEN 3081 DO i = nxlg, nxrg 3082 DO j = nysg, nyng 3083 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3084 ! 3085 !-- Calculate the thrust generated by the nacelle and the tower 3086 tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * & 3087 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3088 tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * & 3089 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3090 3091 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & 3092 - tend_nac_x - tend_tow_x ) & 3093 * MERGE( 1.0_wp, 0.0_wp, & 3094 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 3095 ENDDO 2930 !--------------------------------------------------------------------------------------------------! 2931 SUBROUTINE wtm_actions( location ) 2932 2933 2934 CHARACTER(LEN=*) :: location !< 2935 2936 INTEGER(iwp) :: i !< running index 2937 INTEGER(iwp) :: j !< running index 2938 INTEGER(iwp) :: k !< running index 2939 2940 2941 SELECT CASE ( location ) 2942 2943 CASE ( 'before_timestep' ) 2944 2945 CALL wtm_forces 2946 2947 CASE ( 'u-tendency' ) 2948 ! 2949 2950 !-- Apply the x-component of the force to the u-component of the flow: 2951 IF ( time_since_reference_point >= time_turbine_on ) THEN 2952 DO i = nxlg, nxrg 2953 DO j = nysg, nyng 2954 DO k = nzb + 1, MAXVAL( k_hub ) + MAXVAL( k_smear ) 2955 ! 2956 !-- Calculate the thrust generated by the nacelle and the tower: 2957 tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * SIGN( u(k,j,i)**2 , u(k,j,i) ) 2958 tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * SIGN( u(k,j,i)**2 , u(k,j,i) ) 2959 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) - tend_nac_x - tend_tow_x ) * & 2960 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 3096 2961 ENDDO 3097 2962 ENDDO 3098 ENDIF 3099 3100 CASE ( 'v-tendency' ) 3101 ! 3102 !-- Apply the y-component of the force to the v-component of the flow: 3103 IF ( time_since_reference_point >= time_turbine_on ) THEN 3104 DO i = nxlg, nxrg 3105 DO j = nysg, nyng 3106 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3107 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & 3108 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3109 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & 3110 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3111 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & 3112 - tend_nac_y - tend_tow_y ) & 3113 * MERGE( 1.0_wp, 0.0_wp, & 3114 BTEST( wall_flags_total_0(k,j,i), 2 ) ) 3115 ENDDO 2963 ENDDO 2964 ENDIF 2965 2966 CASE ( 'v-tendency' ) 2967 ! 2968 !-- Apply the y-component of the force to the v-component of the flow: 2969 IF ( time_since_reference_point >= time_turbine_on ) THEN 2970 DO i = nxlg, nxrg 2971 DO j = nysg, nyng 2972 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 2973 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * SIGN( v(k,j,i)**2 , v(k,j,i) ) 2974 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * SIGN( v(k,j,i)**2 , v(k,j,i) ) 2975 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) - tend_nac_y - tend_tow_y ) * & 2976 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 3116 2977 ENDDO 3117 2978 ENDDO 3118 ENDIF 3119 3120 CASE ( 'w-tendency' ) 3121 ! 3122 !-- Apply the z-component of the force to the w-component of the flow: 3123 IF ( time_since_reference_point >= time_turbine_on ) THEN 3124 DO i = nxlg, nxrg 3125 DO j = nysg, nyng 3126 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3127 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) & 3128 * MERGE( 1.0_wp, 0.0_wp, & 3129 BTEST( wall_flags_total_0(k,j,i), 3 ) ) 3130 ENDDO 2979 ENDDO 2980 ENDIF 2981 2982 CASE ( 'w-tendency' ) 2983 ! 2984 !-- Apply the z-component of the force to the w-component of the flow: 2985 IF ( time_since_reference_point >= time_turbine_on ) THEN 2986 DO i = nxlg, nxrg 2987 DO j = nysg, nyng 2988 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 2989 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) * & 2990 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 3131 2991 ENDDO 3132 2992 ENDDO 3133 ENDIF 3134 3135 3136 CASE DEFAULT 3137 CONTINUE 3138 3139 END SELECT 3140 3141 3142 END SUBROUTINE wtm_actions 3143 3144 3145 !------------------------------------------------------------------------------! 2993 ENDDO 2994 ENDIF 2995 2996 2997 CASE DEFAULT 2998 CONTINUE 2999 3000 END SELECT 3001 3002 3003 END SUBROUTINE wtm_actions 3004 3005 3006 !--------------------------------------------------------------------------------------------------! 3146 3007 ! Description: 3147 3008 ! ------------ 3148 !> Application of the additional forces generated by the wind turbine on the 3149 !> flow components(tendency terms)3009 !> Application of the additional forces generated by the wind turbine on the flow components 3010 !> (tendency terms) 3150 3011 !> Call for grid point i,j 3151 !------------------------------------------------------------------------------! 3152 SUBROUTINE wtm_actions_ij( i, j, location ) 3153 3154 3155 CHARACTER (LEN=*) :: location !< 3156 INTEGER(iwp) :: i !< running index 3157 INTEGER(iwp) :: j !< running index 3158 INTEGER(iwp) :: k !< running index 3159 3160 SELECT CASE ( location ) 3161 3162 CASE ( 'before_timestep' ) 3163 3164 CALL wtm_forces 3165 3166 CASE ( 'u-tendency' ) 3167 ! 3168 !-- Apply the x-component of the force to the u-component of the flow: 3169 IF ( time_since_reference_point >= time_turbine_on ) THEN 3170 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3171 ! 3172 !-- Calculate the thrust generated by the nacelle and the tower 3173 tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * & 3174 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3175 tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * & 3176 SIGN( u(k,j,i)**2 , u(k,j,i) ) 3177 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) & 3178 - tend_nac_x - tend_tow_x ) & 3179 * MERGE( 1.0_wp, 0.0_wp, & 3180 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 3012 !--------------------------------------------------------------------------------------------------! 3013 SUBROUTINE wtm_actions_ij( i, j, location ) 3014 3015 3016 CHARACTER (LEN=*) :: location !< 3017 3018 INTEGER(iwp) :: i !< running index 3019 INTEGER(iwp) :: j !< running index 3020 INTEGER(iwp) :: k !< running index 3021 3022 SELECT CASE ( location ) 3023 3024 CASE ( 'before_timestep' ) 3025 3026 CALL wtm_forces 3027 3028 CASE ( 'u-tendency' ) 3029 ! 3030 !-- Apply the x-component of the force to the u-component of the flow: 3031 IF ( time_since_reference_point >= time_turbine_on ) THEN 3032 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3033 ! 3034 !-- Calculate the thrust generated by the nacelle and the tower: 3035 tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) * SIGN( u(k,j,i)**2 , u(k,j,i) ) 3036 tend_tow_x = 0.5_wp * tow_cd_surf(k,j,i) * SIGN( u(k,j,i)**2 , u(k,j,i) ) 3037 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i) - tend_nac_x - tend_tow_x ) * & 3038 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 3039 ENDDO 3040 ENDIF 3041 3042 CASE ( 'v-tendency' ) 3043 ! 3044 !-- Apply the y-component of the force to the v-component of the flow: 3045 IF ( time_since_reference_point >= time_turbine_on ) THEN 3046 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3047 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * SIGN( v(k,j,i)**2 , v(k,j,i) ) 3048 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * SIGN( v(k,j,i)**2 , v(k,j,i) ) 3049 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) - tend_nac_y - tend_tow_y ) * & 3050 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 3181 3051 ENDDO 3182 ENDIF 3183 3184 CASE ( 'v-tendency' ) 3185 ! 3186 !-- Apply the y-component of the force to the v-component of the flow: 3187 IF ( time_since_reference_point >= time_turbine_on ) THEN 3188 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3189 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) * & 3190 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3191 tend_tow_y = 0.5_wp * tow_cd_surf(k,j,i) * & 3192 SIGN( v(k,j,i)**2 , v(k,j,i) ) 3193 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i) & 3194 - tend_nac_y - tend_tow_y ) & 3195 * MERGE( 1.0_wp, 0.0_wp, & 3196 BTEST( wall_flags_total_0(k,j,i), 2 ) ) 3197 ENDDO 3198 ENDIF 3199 3200 CASE ( 'w-tendency' ) 3201 ! 3202 !-- Apply the z-component of the force to the w-component of the flow: 3203 IF ( time_since_reference_point >= time_turbine_on ) THEN 3204 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3205 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) & 3206 * MERGE( 1.0_wp, 0.0_wp, & 3207 BTEST( wall_flags_total_0(k,j,i), 3 ) ) 3208 ENDDO 3209 ENDIF 3210 3211 3212 CASE DEFAULT 3213 CONTINUE 3214 3215 END SELECT 3216 3217 3218 END SUBROUTINE wtm_actions_ij 3219 3220 3221 SUBROUTINE wtm_data_output 3222 3223 3224 INTEGER(iwp) :: t_ind = 0 !< time index 3225 3226 INTEGER(iwp) :: return_value !< returned status value of called function 3227 3228 IF ( myid == 0 ) THEN 3229 3230 ! 3231 !-- At the first call of this routine write the spatial coordinates. 3232 IF ( .NOT. initial_write_coordinates ) THEN 3233 ALLOCATE ( output_values_1d_target(1:n_turbines) ) 3234 output_values_1d_target = hub_x(1:n_turbines) 3235 output_values_1d_pointer => output_values_1d_target 3236 return_value = dom_write_var( nc_filename, & 3237 'x', & 3238 values_realwp_1d = output_values_1d_pointer, & 3239 bounds_start = (/1/), & 3240 bounds_end = (/n_turbines/) ) 3241 3242 output_values_1d_target = hub_y(1:n_turbines) 3243 output_values_1d_pointer => output_values_1d_target 3244 return_value = dom_write_var( nc_filename, & 3245 'y', & 3246 values_realwp_1d = output_values_1d_pointer, & 3247 bounds_start = (/1/), & 3248 bounds_end = (/n_turbines/) ) 3249 3250 output_values_1d_target = hub_z(1:n_turbines) 3251 output_values_1d_pointer => output_values_1d_target 3252 return_value = dom_write_var( nc_filename, & 3253 'z', & 3254 values_realwp_1d = output_values_1d_pointer, & 3255 bounds_start = (/1/), & 3256 bounds_end = (/n_turbines/) ) 3257 3258 output_values_1d_target = rotor_radius(1:n_turbines) * 2.0_wp 3259 output_values_1d_pointer => output_values_1d_target 3260 return_value = dom_write_var( nc_filename, & 3261 'rotor_diameter', & 3262 values_realwp_1d = output_values_1d_pointer, & 3263 bounds_start = (/1/), & 3264 bounds_end = (/n_turbines/) ) 3265 3266 output_values_1d_target = tower_diameter(1:n_turbines) 3267 output_values_1d_pointer => output_values_1d_target 3268 return_value = dom_write_var( nc_filename, & 3269 'tower_diameter', & 3270 values_realwp_1d = output_values_1d_pointer, & 3271 bounds_start = (/1/), & 3272 bounds_end = (/n_turbines/) ) 3273 3274 initial_write_coordinates = .TRUE. 3275 DEALLOCATE ( output_values_1d_target ) 3276 ENDIF 3277 3278 t_ind = t_ind + 1 3279 3052 ENDIF 3053 3054 CASE ( 'w-tendency' ) 3055 ! 3056 !-- Apply the z-component of the force to the w-component of the flow: 3057 IF ( time_since_reference_point >= time_turbine_on ) THEN 3058 DO k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear) 3059 tend(k,j,i) = tend(k,j,i) - rot_tend_z(k,j,i) * & 3060 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 3061 ENDDO 3062 ENDIF 3063 3064 3065 CASE DEFAULT 3066 CONTINUE 3067 3068 END SELECT 3069 3070 3071 END SUBROUTINE wtm_actions_ij 3072 3073 3074 SUBROUTINE wtm_data_output 3075 3076 3077 INTEGER(iwp) :: return_value !< returned status value of called function 3078 INTEGER(iwp) :: t_ind = 0 !< time index 3079 3080 IF ( myid == 0 ) THEN 3081 3082 ! 3083 !-- At the first call of this routine write the spatial coordinates: 3084 IF ( .NOT. initial_write_coordinates ) THEN 3280 3085 ALLOCATE ( output_values_1d_target(1:n_turbines) ) 3281 output_values_1d_target = rotor_speed(:)3086 output_values_1d_target = hub_x(1:n_turbines) 3282 3087 output_values_1d_pointer => output_values_1d_target 3283 3284 return_value = dom_write_var( nc_filename, & 3285 'rotor_speed', & 3286 values_realwp_1d = output_values_1d_pointer, & 3287 bounds_start = (/1, t_ind/), & 3288 bounds_end = (/n_turbines, t_ind /) ) 3289 3290 output_values_1d_target = generator_speed(:) 3291 output_values_1d_pointer => output_values_1d_target 3292 return_value = dom_write_var( nc_filename, & 3293 'generator_speed', & 3294 values_realwp_1d = output_values_1d_pointer, & 3295 bounds_start = (/1, t_ind/), & 3296 bounds_end = (/n_turbines, t_ind /) ) 3297 3298 output_values_1d_target = torque_gen_old(:) 3299 output_values_1d_pointer => output_values_1d_target 3300 3301 return_value = dom_write_var( nc_filename, & 3302 'generator_torque', & 3303 values_realwp_1d = output_values_1d_pointer, & 3304 bounds_start = (/1, t_ind/), & 3305 bounds_end = (/n_turbines, t_ind /) ) 3306 3307 output_values_1d_target = torque_total(:) 3308 output_values_1d_pointer => output_values_1d_target 3309 3310 return_value = dom_write_var( nc_filename, & 3311 'rotor_torque', & 3312 values_realwp_1d = output_values_1d_pointer, & 3313 bounds_start = (/1, t_ind/), & 3314 bounds_end = (/n_turbines, t_ind /) ) 3315 3316 output_values_1d_target = pitch_angle(:) 3317 output_values_1d_pointer => output_values_1d_target 3318 3319 return_value = dom_write_var( nc_filename, & 3320 'pitch_angle', & 3321 values_realwp_1d = output_values_1d_pointer, & 3322 bounds_start = (/1, t_ind/), & 3323 bounds_end = (/n_turbines, t_ind /) ) 3324 3325 output_values_1d_target = torque_gen_old(:)*generator_speed(:)*generator_efficiency 3326 output_values_1d_pointer => output_values_1d_target 3327 3328 return_value = dom_write_var( nc_filename, & 3329 'generator_power', & 3330 values_realwp_1d = output_values_1d_pointer, & 3331 bounds_start = (/1, t_ind/), & 3332 bounds_end = (/n_turbines, t_ind /) ) 3333 3334 DO inot = 1, n_turbines 3335 output_values_1d_target(inot) = torque_total(inot)*rotor_speed(inot)*air_density 3336 ENDDO 3337 output_values_1d_pointer => output_values_1d_target 3338 3339 return_value = dom_write_var( nc_filename, & 3340 'rotor_power', & 3341 values_realwp_1d = output_values_1d_pointer, & 3342 bounds_start = (/1, t_ind/), & 3343 bounds_end = (/n_turbines, t_ind /) ) 3344 3345 output_values_1d_target = thrust_rotor(:) 3346 output_values_1d_pointer => output_values_1d_target 3347 3348 return_value = dom_write_var( nc_filename, & 3349 'rotor_thrust', & 3350 values_realwp_1d = output_values_1d_pointer, & 3351 bounds_start = (/1, t_ind/), & 3352 bounds_end = (/n_turbines, t_ind /) ) 3353 3354 output_values_1d_target = wdir(:)*180.0_wp/pi 3355 output_values_1d_pointer => output_values_1d_target 3356 3357 return_value = dom_write_var( nc_filename, & 3358 'wind_direction', & 3359 values_realwp_1d = output_values_1d_pointer, & 3360 bounds_start = (/1, t_ind/), & 3361 bounds_end = (/n_turbines, t_ind /) ) 3362 3363 output_values_1d_target = (yaw_angle(:))*180.0_wp/pi 3364 output_values_1d_pointer => output_values_1d_target 3365 3366 return_value = dom_write_var( nc_filename, & 3367 'yaw_angle', & 3368 values_realwp_1d = output_values_1d_pointer, & 3369 bounds_start = (/1, t_ind/), & 3370 bounds_end = (/n_turbines, t_ind /) ) 3371 3372 output_values_0d_target = time_since_reference_point 3373 output_values_0d_pointer => output_values_0d_target 3374 3375 return_value = dom_write_var( nc_filename, & 3376 'time', & 3377 values_realwp_0d = output_values_0d_pointer, & 3378 bounds_start = (/t_ind/), & 3379 bounds_end = (/t_ind/) ) 3380 3088 return_value = dom_write_var( nc_filename, & 3089 'x', & 3090 values_realwp_1d = output_values_1d_pointer, & 3091 bounds_start = (/1/), & 3092 bounds_end = (/n_turbines/) ) 3093 3094 output_values_1d_target = hub_y(1:n_turbines) 3095 output_values_1d_pointer => output_values_1d_target 3096 return_value = dom_write_var( nc_filename, & 3097 'y', & 3098 values_realwp_1d = output_values_1d_pointer, & 3099 bounds_start = (/1/), & 3100 bounds_end = (/n_turbines/) ) 3101 3102 output_values_1d_target = hub_z(1:n_turbines) 3103 output_values_1d_pointer => output_values_1d_target 3104 return_value = dom_write_var( nc_filename, & 3105 'z', & 3106 values_realwp_1d = output_values_1d_pointer, & 3107 bounds_start = (/1/), & 3108 bounds_end = (/n_turbines/) ) 3109 3110 output_values_1d_target = rotor_radius(1:n_turbines) * 2.0_wp 3111 output_values_1d_pointer => output_values_1d_target 3112 return_value = dom_write_var( nc_filename, & 3113 'rotor_diameter', & 3114 values_realwp_1d = output_values_1d_pointer, & 3115 bounds_start = (/1/), & 3116 bounds_end = (/n_turbines/) ) 3117 3118 output_values_1d_target = tower_diameter(1:n_turbines) 3119 output_values_1d_pointer => output_values_1d_target 3120 return_value = dom_write_var( nc_filename, & 3121 'tower_diameter', & 3122 values_realwp_1d = output_values_1d_pointer, & 3123 bounds_start = (/1/), & 3124 bounds_end = (/n_turbines/) ) 3125 3126 initial_write_coordinates = .TRUE. 3381 3127 DEALLOCATE ( output_values_1d_target ) 3382 3383 3128 ENDIF 3384 3385 END SUBROUTINE wtm_data_output 3386 3129 3130 t_ind = t_ind + 1 3131 3132 ALLOCATE ( output_values_1d_target(1:n_turbines) ) 3133 output_values_1d_target = rotor_speed(:) 3134 output_values_1d_pointer => output_values_1d_target 3135 3136 return_value = dom_write_var( nc_filename, & 3137 'rotor_speed', & 3138 values_realwp_1d = output_values_1d_pointer, & 3139 bounds_start = (/1, t_ind/), & 3140 bounds_end = (/n_turbines, t_ind /) ) 3141 3142 output_values_1d_target = generator_speed(:) 3143 output_values_1d_pointer => output_values_1d_target 3144 return_value = dom_write_var( nc_filename, & 3145 'generator_speed', & 3146 values_realwp_1d = output_values_1d_pointer, & 3147 bounds_start = (/1, t_ind/), & 3148 bounds_end = (/n_turbines, t_ind /) ) 3149 3150 output_values_1d_target = torque_gen_old(:) 3151 output_values_1d_pointer => output_values_1d_target 3152 3153 return_value = dom_write_var( nc_filename, & 3154 'generator_torque', & 3155 values_realwp_1d = output_values_1d_pointer, & 3156 bounds_start = (/1, t_ind/), & 3157 bounds_end = (/n_turbines, t_ind /) ) 3158 3159 output_values_1d_target = torque_total(:) 3160 output_values_1d_pointer => output_values_1d_target 3161 3162 return_value = dom_write_var( nc_filename, & 3163 'rotor_torque', & 3164 values_realwp_1d = output_values_1d_pointer, & 3165 bounds_start = (/1, t_ind/), & 3166 bounds_end = (/n_turbines, t_ind /) ) 3167 3168 output_values_1d_target = pitch_angle(:) 3169 output_values_1d_pointer => output_values_1d_target 3170 3171 return_value = dom_write_var( nc_filename, & 3172 'pitch_angle', & 3173 values_realwp_1d = output_values_1d_pointer, & 3174 bounds_start = (/1, t_ind/), & 3175 bounds_end = (/n_turbines, t_ind /) ) 3176 3177 output_values_1d_target = torque_gen_old(:) * generator_speed(:) * generator_efficiency 3178 output_values_1d_pointer => output_values_1d_target 3179 3180 return_value = dom_write_var( nc_filename, & 3181 'generator_power', & 3182 values_realwp_1d = output_values_1d_pointer, & 3183 bounds_start = (/1, t_ind/), & 3184 bounds_end = (/n_turbines, t_ind /) ) 3185 3186 DO inot = 1, n_turbines 3187 output_values_1d_target(inot) = torque_total(inot) * rotor_speed(inot) * air_density 3188 ENDDO 3189 output_values_1d_pointer => output_values_1d_target 3190 3191 return_value = dom_write_var( nc_filename, & 3192 'rotor_power', & 3193 values_realwp_1d = output_values_1d_pointer, & 3194 bounds_start = (/1, t_ind/), & 3195 bounds_end = (/n_turbines, t_ind /) ) 3196 3197 output_values_1d_target = thrust_rotor(:) 3198 output_values_1d_pointer => output_values_1d_target 3199 3200 return_value = dom_write_var( nc_filename, & 3201 'rotor_thrust', & 3202 values_realwp_1d = output_values_1d_pointer, & 3203 bounds_start = (/1, t_ind/), & 3204 bounds_end = (/n_turbines, t_ind /) ) 3205 3206 output_values_1d_target = wdir(:)*180.0_wp/pi 3207 output_values_1d_pointer => output_values_1d_target 3208 3209 return_value = dom_write_var( nc_filename, & 3210 'wind_direction', & 3211 values_realwp_1d = output_values_1d_pointer, & 3212 bounds_start = (/1, t_ind/), & 3213 bounds_end = (/n_turbines, t_ind /) ) 3214 3215 output_values_1d_target = (yaw_angle(:)) * 180.0_wp / pi 3216 output_values_1d_pointer => output_values_1d_target 3217 3218 return_value = dom_write_var( nc_filename, & 3219 'yaw_angle', & 3220 values_realwp_1d = output_values_1d_pointer, & 3221 bounds_start = (/1, t_ind/), & 3222 bounds_end = (/n_turbines, t_ind /) ) 3223 3224 output_values_0d_target = time_since_reference_point 3225 output_values_0d_pointer => output_values_0d_target 3226 3227 return_value = dom_write_var( nc_filename, & 3228 'time', & 3229 values_realwp_0d = output_values_0d_pointer, & 3230 bounds_start = (/t_ind/), & 3231 bounds_end = (/t_ind/) ) 3232 3233 DEALLOCATE ( output_values_1d_target ) 3234 3235 ENDIF 3236 3237 END SUBROUTINE wtm_data_output 3238 3387 3239 END MODULE wind_turbine_model_mod
Note: See TracChangeset
for help on using the changeset viewer.
