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

Last change on this file since 4326 was 4326, checked in by oliver.maas, 17 months ago

changed format of turbine control output to allow for higher torque and power values

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