Home

Context Navigation

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

Last change on this file since 4327 was 4326, checked in by oliver.maas, 5 years ago
changed format of turbine control output to allow for higher torque and power values
Property svn:keywords set to `Id`
File size: 118.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |