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

Last change on this file since 4360 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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