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

Last change on this file since 3875 was 3875, checked in by knoop, 2 years ago

Implemented wtm_actions and moved respective module code into it

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