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

Last change on this file since 3698 was 3685, checked in by knoop, 5 years ago

Some interface calls moved to module_interface + cleanup

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