source: palm/trunk/SOURCE/wind_turbine_model_mod.f90 @ 4598

Last change on this file since 4598 was 4537, checked in by oliver.maas, 11 months ago

reset n_turbines_max to 1E4 (10 000), because it was set to 1 000 in r4497 by mistake

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