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

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

replaced <= by < in line 1464 to ensure that ialpha will not be greater than dlen

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