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

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