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

Last change on this file since 3049 was 3049, checked in by Giersch, 6 years ago

Revision history corrected

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