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

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