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

Last change on this file since 4411 was 4411, checked in by maronga, 19 months ago

Added NetCDf output for wind turbine model. Added new features to palmrungui

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