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

Last change on this file since 3226 was 3182, checked in by suehring, 6 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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