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

Last change on this file since 4329 was 4329, checked in by motisi, 16 months ago

Renamed wall_flags_0 to wall_flags_static_0

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