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

Last change on this file since 4431 was 4426, checked in by oliver.maas, 4 years ago

define time as unlimited dimension so that no maximum number of time steps has to be given for wtm_data_output

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