Home

Context Navigation

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

Last change on this file since 4511 was 4497, checked in by raasch, 5 years ago
last bugfix deactivated because of compile problems, files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 143.4 KB

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