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

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

Bugfix in calculation of alpha_attack and Bugfix for restarts in combination with grid stretching

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