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

Last change on this file since 3069 was 3069, checked in by witha, 3 years ago

Initialization of some arrays changed

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