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

Last change on this file since 4465 was 4465, checked in by maronga, 2 years ago

added output of rotor and tower diameters and removed ASCII output in wind turbine model

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