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

Last change on this file since 4444 was 4438, checked in by suehring, 4 years ago

Synthetic turbulence: performance optimizations - random numbers only defined and computed locally, option to compute velocity seeds locally without need of global communication; paralell random number generator: new routine to initialize 1D random number arrays; virtual measurements: CPU-log points added

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