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

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

New vertical stretching procedure has been introduced

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