Home

Context Navigation

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

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