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

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

bugfix for last commit (wind turbine model)

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