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

Last change on this file since 3259 was 3248, checked in by sward, 6 years ago

Minor format changes

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