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

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

Error messages related to the new vertical grid stretching revised, Bugfix in IF statement before error message

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