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

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

added steering for NetCDF output for wind turbine model; minor fix in palmrungui

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