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

Last change on this file since 4087 was 4056, checked in by Giersch, 5 years ago

Redundant example files deleted, wtm example added as a test case, bugfix in wtm

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