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

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

some routines instrumented with openmp directives, loop reordering for performance optimization

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