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

Last change on this file since 3046 was 3046, checked in by Giersch, 3 years ago

Remaining error messages revised, comments extended

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