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

Last change on this file since 4423 was 4423, checked in by maronga, 16 months ago

Switched back to serial NetCDF output for wind turbine output

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