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

Last change on this file since 4434 was 4434, checked in by oliver.maas, 21 months ago

added optional netcdf data input for wtm array input parameters

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