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

Last change on this file since 3680 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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