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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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