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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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