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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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