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

Last change on this file since 2894 was 2894, checked in by Giersch, 4 years ago

Reading/Writing? data in case of restart runs revised

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