source: palm/trunk/SOURCE/wind_turbine_model_mod.f90

Last change on this file was 4843, checked in by raasch, 9 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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