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

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

Bugfix: shifted netcdf preprocessor directive to correct position

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