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

Last change on this file since 3241 was 3241, checked in by raasch, 3 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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