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

Last change on this file since 3352 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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