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

Last change on this file since 4457 was 4457, checked in by raasch, 3 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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