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

Last change on this file since 4147 was 4144, checked in by raasch, 5 years ago

relational operators .EQ., .NE., etc. replaced by ==, /=, etc.

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