Home

Context Navigation

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

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