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

Last change on this file since 4459 was 4459, checked in by oliver.maas, 3 years ago

avoid division by zero in tip loss correction factor calculation

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