Home

Context Navigation

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

Last change on this file since 4526 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

Rev	Line
[1820]	1	!> @file wind_turbine_model_mod.f90
[4497]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1819]	4	!
[4497]	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.
[1819]	8	!
[4497]	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.
[1819]	12	!
[4497]	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/>.
[1819]	15	!
[4481]	16	! Copyright 2009-2020 Carl von Ossietzky Universitaet Oldenburg
[4360]	17	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4497]	18	!--------------------------------------------------------------------------------------------------!
[1819]	19	!
[4497]	20	!
[1819]	21	! Current revisions:
	22	! -----------------
[1913]	23	!
[4438]	24	!
[1913]	25	! Former revisions:
	26	! -----------------
	27	! $Id: wind_turbine_model_mod.f90 4497 2020-04-15 10:20:51Z oliver.maas $
[4497]	28	! file re-formatted to follow the PALM coding standard
	29	!
	30	!
	31	! 4495 2020-04-13 20:11:20Z raasch
[4495]	32	! restart data handling with MPI-IO added
	33	!
	34	! 4481 2020-03-31 18:55:54Z maronga
[4467]	35	! ASCII output cleanup
[4497]	36	!
[4467]	37	! 4465 2020-03-20 11:35:48Z maronga
[4497]	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	!
[4465]	41	! 4460 2020-03-12 16:47:30Z oliver.maas
[4460]	42	! allow for simulating up to 10 000 wind turbines
[4497]	43	!
[4460]	44	! 4459 2020-03-12 09:35:23Z oliver.maas
[4459]	45	! avoid division by zero in tip loss correction factor calculation
[4497]	46	!
[4459]	47	! 4457 2020-03-11 14:20:43Z raasch
[4457]	48	! use statement for exchange horiz added
[4497]	49	!
[4457]	50	! 4447 2020-03-06 11:05:30Z oliver.maas
[4447]	51	! renamed wind_turbine_parameters namelist variables
[4497]	52	!
[4447]	53	! 4438 2020-03-03 20:49:28Z suehring
[4438]	54	! Bugfix: shifted netcdf preprocessor directive to correct position
[4497]	55	!
[4438]	56	! 4436 2020-03-03 12:47:02Z oliver.maas
[4434]	57	! added optional netcdf data input for wtm array input parameters
[4497]	58	!
[4434]	59	! 4426 2020-02-27 10:02:19Z oliver.maas
[4497]	60	! define time as unlimited dimension so that no maximum number of time steps has to be given for
	61	! wtm_data_output
	62	!
[4426]	63	! 4423 2020-02-25 07:17:11Z maronga
[4423]	64	! Switched to serial output as data is aggerated before anyway.
[4497]	65	!
[4423]	66	! 4420 2020-02-24 14:13:56Z maronga
[4420]	67	! Added output control for wind turbine model
[4497]	68	!
[4420]	69	! 4412 2020-02-19 14:53:13Z maronga
[4412]	70	! Bugfix: corrected character length in dimension_names
[4497]	71	!
[4412]	72	! 4411 2020-02-18 14:28:02Z maronga
[4411]	73	! Added output in NetCDF format using DOM (only netcdf4-parallel is supported).
	74	! Old ASCII output is still available at the moment.
[4497]	75	!
[4411]	76	! 4360 2020-01-07 11:25:50Z suehring
[4497]	77	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
	78	! information used in wall_flags_static_0
	79	!
[4346]	80	! 4343 2019-12-17 12:26:12Z oliver.maas
[4497]	81	! replaced <= by < in line 1464 to ensure that ialpha will not be greater than dlen
	82	!
[4343]	83	! 4329 2019-12-10 15:46:36Z motisi
[4329]	84	! Renamed wall_flags_0 to wall_flags_static_0
[4497]	85	!
[4329]	86	! 4326 2019-12-06 14:16:14Z oliver.maas
[4497]	87	! changed format of turbine control output to allow for higher torque and power values
	88	!
[4326]	89	! 4182 2019-08-22 15:20:23Z scharf
[4182]	90	! Corrected "Former revisions" section
[4497]	91	!
[4182]	92	! 4144 2019-08-06 09:11:47Z raasch
[4144]	93	! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
[4497]	94	!
[4144]	95	! 4056 2019-06-27 13:53:16Z Giersch
[4497]	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	!
[4056]	100	! 3885 2019-04-11 11:29:34Z kanani
[4497]	101	! Changes related to global restructuring of location messages and introduction of additional debug
	102	! messages
	103	!
[3885]	104	! 3875 2019-04-08 17:35:12Z knoop
[3875]	105	! Addaped wtm_tendency to fit the module actions interface
[4497]	106	!
[3875]	107	! 3832 2019-03-28 13:16:58Z raasch
[3832]	108	! instrumented with openmp directives
[4497]	109	!
[3832]	110	! 3725 2019-02-07 10:11:02Z raasch
[3725]	111	! unused variables removed
[4497]	112	!
[3725]	113	! 3685 2019-01-21 01:02:11Z knoop
[3685]	114	! Some interface calls moved to module_interface + cleanup
[4497]	115	!
[3685]	116	! 3655 2019-01-07 16:51:22Z knoop
[3593]	117	! Replace degree symbol by 'degrees'
[4497]	118	!
[4182]	119	! 1914 2016-05-26 14:44:07Z witha
	120	! Initial revision
[3593]	121	!
[4497]	122	!
[1819]	123	! Description:
	124	! ------------
[4497]	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).
[1819]	127	!> The wind turbines include the tower effect, can be yawed and tilted.
[4497]	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).
[1819]	132	!>
[3065]	133	!> @todo Replace dz(1) appropriatly to account for grid stretching
[1819]	134	!> @todo Revise code according to PALM Coding Standard
	135	!> @todo Implement ADM and ALM turbine models
	136	!> @todo Generate header information
[1917]	137	!> @todo Implement further parameter checks and error messages
[1819]	138	!> @todo Revise and add code documentation
	139	!> @todo Output turbine parameters as timeseries
	140	!> @todo Include additional output variables
[1864]	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
[1819]	144	!
[4497]	145	!--------------------------------------------------------------------------------------------------!
[1819]	146	MODULE wind_turbine_model_mod
	147
[4497]	148	USE arrays_3d, &
[2553]	149	ONLY: tend, u, v, w, zu, zw
[1819]	150
[4497]	151	USE basic_constants_and_equations_mod, &
[2553]	152	ONLY: pi
[1819]	153
[4497]	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
[1819]	160
[4497]	161	USE cpulog, &
[1819]	162	ONLY: cpu_log, log_point_s
	163
[4411]	164	USE data_output_module
[4497]	165
	166	USE grid_variables, &
[1819]	167	ONLY: ddx, dx, ddy, dy
	168
[4497]	169	USE indices, &
	170	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, &
[4346]	171	nzb, nzt, wall_flags_total_0
[1819]	172
	173	USE kinds
	174
[4497]	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, &
[4434]	186	vars_pids
[4497]	187
[1819]	188	USE pegrid
	189
[4495]	190	USE restart_data_mpi_io_mod, &
	191	ONLY: rrd_mpi_io_global_array, wrd_mpi_io_global_array
[1819]	192
[4497]	193
[1819]	194	IMPLICIT NONE
	195
[1864]	196	PRIVATE
[1819]	197
[4497]	198	CHARACTER(LEN=100) :: variable_name !< name of output variable
	199	CHARACTER(LEN=30) :: nc_filename !<
[4460]	200
[4497]	201
	202	INTEGER(iwp), PARAMETER :: n_turbines_max = 1000 !< maximum number of turbines (for array allocation)
	203
	204
[1819]	205	!
	206	!-- Variables specified in the namelist wind_turbine_par
[4497]	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
[1819]	209
[4497]	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
[4447]	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
[4497]	234	REAL(wp) :: tilt_angle = 0.0_wp !< vertical tilt_angle of the rotor [degree] ( positive = backwards )
[1819]	235
[4497]	236	!< ( clockwise, 0 = facing west )
[4460]	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
[4497]	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]
[1839]	248
[4497]	249
[1819]	250	!
	251	!-- Variables specified in the namelist for speed controller
	252	!-- Default values are from the NREL 5MW research turbine (Jonkman, 2008)
[4497]	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]
[1819]	267
[1839]	268
[1912]	269
[1819]	270	!
[4497]	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]
[1819]	275
	276	!
[4497]	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 !<
[1819]	283
[4497]	284	INTEGER(iwp), DIMENSION(1) :: lct !<
[1819]	285
[4497]	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
[1819]	294
[4497]	295	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nsegs !< number of segments per ring and turbine
[1819]	296
[1912]	297	!
[4497]	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
[1839]	305	!
[4497]	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 !<
[1819]	311
[4497]	312	REAL(wp) :: accu_cl_cd_tab = 0.1_wp !< Accuracy of the interpolation of the lift and drag coeff [deg]
[1819]	313
[4497]	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
[1912]	318	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tow_cd_surf !< 3d field of the nacelle drag coefficient
[1819]	319
	320	!
[4497]	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 !<
[3832]	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)
[1912]	335	!
[4497]	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 !<
[1912]	340	REAL(wp) :: tend_tow_y = 0.0_wp !<
[3832]	341	!$OMP THREADPRIVATE (tend_nac_x, tend_tow_x, tend_nac_y, tend_tow_y)
[1819]	342
[4497]	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 !<
[1912]	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	!
[4497]	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
[1912]	365	!
[4497]	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 !<
[1912]	370
	371	!
[4497]	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
[1912]	377	!
[4497]	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 !<
[1912]	382
	383	!
[4497]	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:
[4460]	391	REAL(wp), DIMENSION(1:n_turbines_max,1:3,1:3) :: rot_coord_trans !< matrix for rotation of rotor coordinates
[1819]	392
[4497]	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
[1839]	403	!
[4497]	404	!-- Fixed variables for the speed controller:
	405	LOGICAL :: start_up = .TRUE. !<
[1819]	406
[4497]	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
[1819]	415
[4497]	416	REAL(wp), DIMENSION(n_turbines_max) :: rotor_speed_l = 0.0_wp !< local rot speed [rad/s]
[2563]	417
[1839]	418	!
[4497]	419	!-- Fixed variables for the yaw controller:
	420	INTEGER(iwp) :: wdlon !<
	421	INTEGER(iwp) :: wdsho !<
[1819]	422
[4497]	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
[1912]	431	REAL(wp), DIMENSION(:) , ALLOCATABLE :: yawdir !< direction to yaw
[4447]	432	REAL(wp), DIMENSION(:) , ALLOCATABLE :: yaw_angle_l !< local (cpu) yaw angle
[4497]	433
	434	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd2 !<
[1912]	435	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wd30 !<
[1819]	436
[4497]	437
[2563]	438	!
[4497]	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
[4460]	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
[1819]	447
[2563]	448
[1839]	449	SAVE
[1819]	450
[1839]	451
	452	INTERFACE wtm_parin
	453	MODULE PROCEDURE wtm_parin
	454	END INTERFACE wtm_parin
[2563]	455
[1912]	456	INTERFACE wtm_check_parameters
	457	MODULE PROCEDURE wtm_check_parameters
	458	END INTERFACE wtm_check_parameters
[3875]	459
[4411]	460	INTERFACE wtm_data_output
	461	MODULE PROCEDURE wtm_data_output
	462	END INTERFACE wtm_data_output
[4497]	463
[1839]	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
[2553]	471
[4411]	472	INTERFACE wtm_init_output
	473	MODULE PROCEDURE wtm_init_output
	474	END INTERFACE wtm_init_output
[4497]	475
[3875]	476	INTERFACE wtm_actions
	477	MODULE PROCEDURE wtm_actions
	478	MODULE PROCEDURE wtm_actions_ij
	479	END INTERFACE wtm_actions
[1819]	480
[3875]	481	INTERFACE wtm_rrd_global
[4495]	482	MODULE PROCEDURE wtm_rrd_global_ftn
	483	MODULE PROCEDURE wtm_rrd_global_mpi
[3875]	484	END INTERFACE wtm_rrd_global
[1819]	485
[3875]	486	INTERFACE wtm_wrd_global
	487	MODULE PROCEDURE wtm_wrd_global
	488	END INTERFACE wtm_wrd_global
[2563]	489
[3875]	490
[4497]	491	PUBLIC &
	492	dt_data_output_wtm, &
	493	time_wtm, &
[4420]	494	wind_turbine
	495
[4497]	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, &
[3875]	505	wtm_wrd_global
	506
	507
[1819]	508	CONTAINS
	509
	510
[4497]	511	!--------------------------------------------------------------------------------------------------!
[1819]	512	! Description:
	513	! ------------
[1839]	514	!> Parin for &wind_turbine_par for wind turbine model
[4497]	515	!--------------------------------------------------------------------------------------------------!
	516	SUBROUTINE wtm_parin
[1819]	517
[4497]	518	IMPLICIT NONE
[1819]	519
[4497]	520	CHARACTER(LEN=80) :: line !< dummy string that contains the current line of the parameter file
[1839]	521
[4497]	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
[4447]	535
[1819]	536	!
[4497]	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 )
[1839]	544
	545	!
[4497]	546	!-- Read user-defined namelist:
	547	READ ( 11, wind_turbine_parameters, ERR = 10, END = 12 )
[2932]	548	!
[4497]	549	!-- Set flag that indicates that the wind turbine model is switched on:
	550	wind_turbine = .TRUE.
[4465]	551
[4497]	552	GOTO 12
[2932]	553
[4497]	554	10 BACKSPACE ( 11 )
	555	READ ( 11, '(A)' ) line
	556	CALL parin_fail_message( 'wind_turbine_parameters', line )
[3246]	557
[4497]	558	12 CONTINUE ! TBD Change from continue, mit ierrn machen
[1839]	559
[4497]	560	END SUBROUTINE wtm_parin
[1839]	561
[2563]	562
[4497]	563	!--------------------------------------------------------------------------------------------------!
[2563]	564	! Description:
	565	! ------------
[2894]	566	!> This routine writes the respective restart data.
[4497]	567	!--------------------------------------------------------------------------------------------------!
	568	SUBROUTINE wtm_wrd_global
[2576]	569
	570
[4497]	571	IMPLICIT NONE
[2776]	572
[2894]	573
[4497]	574	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
[2894]	575
[4497]	576	CALL wrd_write_string( 'generator_speed' )
	577	WRITE ( 14 ) generator_speed
[2894]	578
[4497]	579	CALL wrd_write_string( 'generator_speed_f' )
	580	WRITE ( 14 ) generator_speed_f
[2894]	581
[4497]	582	CALL wrd_write_string( 'generator_speed_f_old' )
	583	WRITE ( 14 ) generator_speed_f_old
[2894]	584
[4497]	585	CALL wrd_write_string( 'generator_speed_old' )
	586	WRITE ( 14 ) generator_speed_old
[2894]	587
[4497]	588	CALL wrd_write_string( 'rotor_speed' )
	589	WRITE ( 14 ) rotor_speed
[2894]	590
[4497]	591	CALL wrd_write_string( 'yaw_angle' )
	592	WRITE ( 14 ) yaw_angle
[2894]	593
[4497]	594	CALL wrd_write_string( 'pitch_angle' )
	595	WRITE ( 14 ) pitch_angle
[2894]	596
[4497]	597	CALL wrd_write_string( 'pitch_angle_old' )
	598	WRITE ( 14 ) pitch_angle_old
[2894]	599
[4497]	600	CALL wrd_write_string( 'torque_gen' )
	601	WRITE ( 14 ) torque_gen
[4495]	602
[4497]	603	CALL wrd_write_string( 'torque_gen_old' )
	604	WRITE ( 14 ) torque_gen_old
[4495]	605
[4497]	606	ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' ) THEN
[4495]	607
[4497]	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 )
[2563]	618
[4497]	619	ENDIF
[2563]	620
[4497]	621	END SUBROUTINE wtm_wrd_global
	622
	623
	624	!--------------------------------------------------------------------------------------------------!
[2563]	625	! Description:
	626	! ------------
[4495]	627	!> Read module-specific global restart data (Fortran binary format).
[4497]	628	!--------------------------------------------------------------------------------------------------!
[4495]	629	SUBROUTINE wtm_rrd_global_ftn( found )
[2563]	630
	631
[4497]	632	USE control_parameters, &
[2894]	633	ONLY: length, restart_string
	634
	635
[2563]	636	IMPLICIT NONE
	637
[4497]	638	LOGICAL, INTENT(OUT) :: found
[2563]	639
	640
[2894]	641	found = .TRUE.
[2563]	642
	643
[2894]	644	SELECT CASE ( restart_string(1:length) )
[2563]	645
[4447]	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
[2894]	662	CASE ( 'torque_gen' )
	663	READ ( 13 ) torque_gen
	664	CASE ( 'torque_gen_old' )
	665	READ ( 13 ) torque_gen_old
[2563]	666
[2894]	667	CASE DEFAULT
[2563]	668
[2894]	669	found = .FALSE.
[2563]	670
[2894]	671	END SELECT
[2563]	672
[4497]	673
[4495]	674	END SUBROUTINE wtm_rrd_global_ftn
[2894]	675
	676
[2563]	677	!------------------------------------------------------------------------------!
	678	! Description:
	679	! ------------
[4495]	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
[4497]	698	!--------------------------------------------------------------------------------------------------!
[4495]	699	! Description:
	700	! ------------
[2563]	701	!> Check namelist parameter
[4497]	702	!--------------------------------------------------------------------------------------------------!
	703	SUBROUTINE wtm_check_parameters
[1912]	704
[4497]	705	IMPLICIT NONE
[4436]	706
[4497]	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
[4436]	712
[4497]	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
[4436]	718
[4497]	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
[4436]	722
[4497]	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 )
[1912]	725	ENDIF
[4497]	726	ENDIF
	727
	728	END SUBROUTINE wtm_check_parameters
	729	!
	730	!--------------------------------------------------------------------------------------------------!
[1839]	731	! Description:
	732	! ------------
	733	!> Allocate wind turbine model arrays
[4497]	734	!--------------------------------------------------------------------------------------------------!
	735	SUBROUTINE wtm_init_arrays
[1839]	736
[4497]	737	IMPLICIT NONE
[1839]	738
[4497]	739	REAL(wp) :: delta_r_factor !<
	740	REAL(wp) :: delta_r_init !<
[1864]	741
[4436]	742	#if defined( __netcdf )
[1839]	743	!
[4497]	744	! Read wtm input file (netcdf) if it exists:
	745	IF ( input_pids_wtm ) THEN
[4434]	746
	747	!
[4497]	748	!-- Open the wtm input file:
	749	CALL open_read_file( TRIM( input_file_wtm ) // &
	750	TRIM( coupling_char ), pids_id )
[4434]	751
[4497]	752	CALL inquire_num_variables( pids_id, num_var_pids )
[4434]	753
	754	!
[4497]	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 )
[4434]	758
	759	!
[4497]	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
[4434]	764
[4497]	765	IF ( check_existence( vars_pids, 'rotor_speed' ) ) THEN
	766	CALL get_variable( pids_id, 'rotor_speed', rotor_speed(1:n_turbines) )
	767	ENDIF
[4434]	768
[4497]	769	IF ( check_existence( vars_pids, 'pitch_angle' ) ) THEN
	770	CALL get_variable( pids_id, 'pitch_angle', pitch_angle(1:n_turbines) )
	771	ENDIF
[4434]	772
[4497]	773	IF ( check_existence( vars_pids, 'yaw_angle' ) ) THEN
	774	CALL get_variable( pids_id, 'yaw_angle', yaw_angle(1:n_turbines) )
	775	ENDIF
[4434]	776
[4497]	777	IF ( check_existence( vars_pids, 'hub_x' ) ) THEN
	778	CALL get_variable( pids_id, 'hub_x', hub_x(1:n_turbines) )
	779	ENDIF
[4434]	780
[4497]	781	IF ( check_existence( vars_pids, 'hub_y' ) ) THEN
	782	CALL get_variable( pids_id, 'hub_y', hub_y(1:n_turbines) )
	783	ENDIF
[4434]	784
[4497]	785	IF ( check_existence( vars_pids, 'hub_z' ) ) THEN
	786	CALL get_variable( pids_id, 'hub_z', hub_z(1:n_turbines) )
	787	ENDIF
[4434]	788
[4497]	789	IF ( check_existence( vars_pids, 'nacelle_radius' ) ) THEN
	790	CALL get_variable( pids_id, 'nacelle_radius', nacelle_radius(1:n_turbines) )
	791	ENDIF
[4434]	792
[4497]	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
[4434]	800
[4497]	801	IF ( check_existence( vars_pids, 'tower_cd' ) ) THEN
	802	CALL get_variable( pids_id, 'tower_cd', tower_cd(1:n_turbines) )
	803	ENDIF
[4434]	804	!
[4497]	805	!-- Close wtm input file:
	806	CALL close_input_file( pids_id )
[4434]	807
[4497]	808	ENDIF
[4434]	809	#endif
	810
	811	!
[4497]	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) )
[1819]	816
[4497]	817	nrings(:) = 0
	818	delta_r(:) = 0.0_wp
[1819]	819
	820	!
[4497]	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) )
[1819]	825
[4497]	826	DO inot = 1, n_turbines
[1819]	827	!
[4497]	828	!-- Determine number of rings:
	829	nrings(inot) = NINT( rotor_radius(inot) / delta_r_init )
[1819]	830
[4497]	831	delta_r(inot) = rotor_radius(inot) / nrings(inot)
[1819]	832
[4497]	833	ENDDO
[1819]	834
[4497]	835	nrings_max = MAXVAL( nrings )
[1819]	836
[4497]	837	ALLOCATE( nsegs(1:nrings_max,1:n_turbines) )
	838	ALLOCATE( nsegs_total(1:n_turbines) )
[1819]	839
[4497]	840	nsegs(:,:) = 0
	841	nsegs_total(:) = 0
[1819]	842
	843
[4497]	844	DO inot = 1, n_turbines
	845	DO ring = 1, nrings(inot)
[1819]	846	!
[4497]	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
[1819]	851	!
[4497]	852	!-- Total sum of all rotor segments:
	853	nsegs_total(inot) = SUM( nsegs(:,inot) )
	854	ENDDO
[1819]	855
	856	!
[4497]	857	!-- Maximum number of segments per ring:
	858	nsegs_max = MAXVAL( nsegs )
[1819]	859
[1864]	860	!!
[4497]	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
[1819]	871
	872
	873	!
[4497]	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) )
[1819]	883
	884	!
[4497]	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
[1819]	893	!
[4497]	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) )
[1819]	905
	906	!
[4497]	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) )
[1819]	914
	915	!
[4497]	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) )
[1819]	920
	921	!
[4497]	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) )
[1819]	931
	932	!
[4497]	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) )
[1819]	940
	941	!
[4497]	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
[1819]	949
[4497]	950	torque_total(:) = 0.0_wp
	951	thrust_rotor(:) = 0.0_wp
[2563]	952
[4497]	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
[1819]	962
[4497]	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
[1819]	970	!
[4497]	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
[1819]	982
[4497]	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
[1819]	989
[4497]	990	rotx(:,:) = 0.0_wp
	991	roty(:,:) = 0.0_wp
	992	rotz(:,:) = 0.0_wp
[1819]	993
[4497]	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
[1819]	1002
[4497]	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
[1819]	1009
	1010
[4497]	1011	END SUBROUTINE wtm_init_arrays
[1819]	1012
	1013
[4497]	1014	!--------------------------------------------------------------------------------------------------!
[1819]	1015	! Description:
	1016	! ------------
[1839]	1017	!> Initialization of the wind turbine model
[4497]	1018	!--------------------------------------------------------------------------------------------------!
	1019	SUBROUTINE wtm_init
[1819]	1020
[4457]	1021
[4497]	1022	USE control_parameters, &
	1023	ONLY: dz_stretch_level_start
[1819]	1024
[4497]	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
[1819]	1037	!
[4497]	1038	!-- Help variables for the smearing function:
	1039	REAL(wp) :: eps_kernel !<
	1040
[1839]	1041	!
[4497]	1042	!-- Help variables for calculation of the tower drag:
	1043	INTEGER(iwp) :: tower_n !<
	1044	INTEGER(iwp) :: tower_s !<
[1819]	1045
[4497]	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
[1839]	1061	!
[4497]	1062	!--------------------------------------------------------------------------------------------------!
	1063	!-- Calculation of parameters for the regularization kernel (smearing of the forces)
	1064	!--------------------------------------------------------------------------------------------------!
[1819]	1065	!
[4497]	1066	!-- In the following, some of the required parameters for the smearing will be calculated:
[1819]	1067
[4497]	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
[1819]	1071	!
[4497]	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
[1819]	1077	!
[4497]	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
[3065]	1090	!
[4497]	1091	!-- Square of eps_min:
	1092	eps_min2 = eps_min**2
[1819]	1093	!
[4497]	1094	!-- Parameters in the polynomial function:
	1095	pol_a = 1.0_wp / eps_min**4
	1096	pol_b = 2.0_wp / eps_min**2
[1819]	1097	!
[4497]	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
[1819]	1104	!
[4497]	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
[1819]	1109
[1864]	1110
[4497]	1111	DO inot = 1, n_turbines
[1819]	1112	!
[4497]	1113	!-- Rotate the rotor coordinates in case yaw and tilt are defined:
	1114	CALL wtm_rotate_rotor( inot )
	1115
[1819]	1116	!
[4497]	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) )
[1819]	1121
	1122	!
[4497]	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) )
[1819]	1129
[4497]	1130	ENDDO
	1131
[1819]	1132	!
[4497]	1133	!-- Call the wtm_init_speed_control subroutine and calculate the local rotor_speed for the
	1134	!-- respective processor:
	1135	IF ( speed_control) THEN
[2792]	1136
[4497]	1137	CALL wtm_init_speed_control
[2792]	1138
[4497]	1139	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
[2792]	1140
[4497]	1141	DO inot = 1, n_turbines
[2792]	1142
[4497]	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
[2792]	1148
[4497]	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
[2792]	1154
[4497]	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
[2792]	1160
[4497]	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
[2792]	1166
[4497]	1167	IF ( ( nxl <= i_hub(inot) ) .AND. ( nxr >= i_hub(inot) ) ) THEN
	1168	IF ( ( nys <= j_hub(inot) ) .AND. ( nyn >= j_hub(inot) ) ) THEN
[2792]	1169
[4497]	1170	rotor_speed_l(inot) = generator_speed(inot) / gear_ratio
	1171
[2792]	1172	ENDIF
[4497]	1173	ENDIF
[2792]	1174
[4497]	1175	END DO
[2792]	1176
	1177	ENDIF
	1178
[4497]	1179	ENDIF
	1180
[2792]	1181	!
[4497]	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	!--------------------------------------------------------------------------------------------------!
[1819]	1186	!
[4497]	1187	!-- Note: so far this is only a 2D version, in that the mean flow is perpendicular to the rotor
	1188	!-- area.
[1819]	1189
	1190	!
[4497]	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
[1819]	1194
[4497]	1195	ALLOCATE( circle_points(1:2,1:upper_end) )
[1819]	1196
[4497]	1197	circle_points(:,:) = 0.0_wp
	1198
	1199
	1200	DO inot = 1, n_turbines ! loop over number of turbines
[1819]	1201	!
[4497]	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)
[1819]	1206
	1207	!
[4497]	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 )
[1819]	1211	!
[4497]	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
[1819]	1218
	1219	!
[4497]	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
[1819]	1224
	1225	!
[4497]	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
[1819]	1229
[4497]	1230	DO j = nys, nyn
[1819]	1231	!
[4497]	1232	!-- Loop from south to north boundary of tower:
	1233	IF ( ( j >= tower_s ) .AND. ( j <= tower_n ) ) THEN
[1839]	1234
[4497]	1235	DO k = nzb, nzt
[1819]	1236
[4497]	1237	IF ( k == k_hub(inot) ) THEN
	1238	IF ( tower_n - tower_s >= 1 ) THEN
[1839]	1239	!
[4497]	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)
[1819]	1254	!
[4497]	1255	!-- Grid boxes inbetween (where tow_cd_surf = grid box area):
[1819]	1256	ELSE
[4497]	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)
[1819]	1260	ENDIF
[1839]	1261	!
[4497]	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
[1839]	1268	!
[4497]	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
[1839]	1273	!
[4497]	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)
[1839]	1284	!
[4497]	1285	!-- Grid boxes inbetween (where tow_cd_surf = grid box area):
[1819]	1286	ELSE
[4497]	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
[1839]	1294
[4497]	1295	ENDIF ! end if k == k_hub
[1839]	1296
[4497]	1297	ENDDO ! end loop over k
[1839]	1298
[4497]	1299	ENDIF ! end if inside north and south boundary of tower
[1839]	1300
[4497]	1301	ENDDO ! end loop over j
[1839]	1302
[4497]	1303	ENDIF ! end if hub inside domain + ghostpoints
[1819]	1304
[4497]	1305
	1306	CALL exchange_horiz( tow_cd_surf, nbgp )
	1307
[1839]	1308	!
[4497]	1309	!-- Calculation of the nacelle area
	1310	!-- CAUTION: Currently disabled due to segmentation faults on the FLOW HPC cluster (Oldenburg)
[1864]	1311	!!
[4497]	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):
[1864]	1314	!
[4497]	1315	! dy_int = dy / 10000.0_wp
[1864]	1316	!
[4497]	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
[1864]	1321	!!
[4497]	1322	!!-- bottom intersection point
	1323	! circle_points(1,i_ip) = hub_z(inot) - SQRT( sqrt_arg )
[1864]	1324	!!
[4497]	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
[1864]	1331	!
	1332	!
[4497]	1333	! DO j = nys, nyn
[1864]	1334	!!
[4497]	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
[1864]	1340	!!
[4497]	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
[1864]	1346	!!
[4497]	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
[1864]	1354	!!
[4497]	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
[1819]	1385
[4497]	1386	ENDDO ! end of loop over turbines
[1819]	1387
[4497]	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
[1819]	1390
[4497]	1391	CALL wtm_read_blade_tables
[3685]	1392
[4497]	1393	IF ( debug_output ) CALL debug_message( 'wtm_init', 'end' )
[1864]	1394
[4497]	1395	END SUBROUTINE wtm_init
[1864]	1396
[4497]	1397
	1398
[4423]	1399	SUBROUTINE wtm_init_output
[4497]	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
[4411]	1411	!
[4497]	1412	!-- Create NetCDF output file:
[4465]	1413	#if defined( __netcdf4 )
	1414	nc_filename = 'DATA_1D_TS_WTM_NETCDF' // TRIM( coupling_char )
[4423]	1415	return_value = dom_def_file( nc_filename, 'netcdf4-serial' )
[4465]	1416	#else
[4497]	1417	message_string = 'Wind turbine model output requires NetCDF version 4. ' // &
[4465]	1418	'No output file will be created.'
[4497]	1419	CALL message( 'wtm_init_output', 'PA0672', 0, 1, 0, 6, 0 )
[4465]	1420	#endif
[4423]	1421	IF ( myid == 0 ) THEN
[4411]	1422	!
[4497]	1423	!-- Define dimensions in output file:
[4447]	1424	ALLOCATE( ndim(1:n_turbines) )
	1425	DO n = 1, n_turbines
[4411]	1426	ndim(n) = n
	1427	ENDDO
[4497]	1428	return_value = dom_def_dim( nc_filename, &
	1429	dimension_name = 'turbine', &
	1430	output_type = 'int32', &
	1431	bounds = (/1_iwp, n_turbines/), &
[4411]	1432	values_int32 = ndim )
	1433	DEALLOCATE( ndim )
	1434
[4497]	1435	return_value = dom_def_dim( nc_filename, &
	1436	dimension_name = 'time', &
	1437	output_type = 'real32', &
	1438	bounds = (/1_iwp/), &
[4426]	1439	values_realwp = (/0.0_wp/) )
[4497]	1440
[4411]	1441	variable_name = 'x'
[4497]	1442	return_value = dom_def_var( nc_filename, &
	1443	variable_name = variable_name, &
	1444	dimension_names = (/'turbine'/), &
[4411]	1445	output_type = 'real32' )
[4465]	1446
[4411]	1447	variable_name = 'y'
[4497]	1448	return_value = dom_def_var( nc_filename, &
	1449	variable_name = variable_name, &
	1450	dimension_names = (/'turbine'/), &
[4411]	1451	output_type = 'real32' )
	1452
	1453	variable_name = 'z'
[4497]	1454	return_value = dom_def_var( nc_filename, &
	1455	variable_name = variable_name, &
	1456	dimension_names = (/'turbine'/), &
[4411]	1457	output_type = 'real32' )
	1458
[4465]	1459
	1460	variable_name = 'rotor_diameter'
[4497]	1461	return_value = dom_def_var( nc_filename, &
	1462	variable_name = variable_name, &
	1463	dimension_names = (/'turbine'/), &
[4465]	1464	output_type = 'real32' )
	1465
	1466	variable_name = 'tower_diameter'
[4497]	1467	return_value = dom_def_var( nc_filename, &
	1468	variable_name = variable_name, &
	1469	dimension_names = (/'turbine'/), &
[4465]	1470	output_type = 'real32' )
	1471
[4497]	1472	return_value = dom_def_att( nc_filename, &
	1473	variable_name = 'time', &
	1474	attribute_name = 'units', &
[4411]	1475	value = 'seconds since ' // origin_date_time )
[4465]	1476
[4497]	1477	return_value = dom_def_att( nc_filename, &
	1478	variable_name = 'x', &
	1479	attribute_name = 'units', &
[4411]	1480	value = 'm' )
	1481
[4497]	1482	return_value = dom_def_att( nc_filename, &
	1483	variable_name = 'y', &
	1484	attribute_name = 'units', &
	1485	value = 'm' )
[4411]	1486
[4497]	1487	return_value = dom_def_att( nc_filename, &
	1488	variable_name = 'z', &
	1489	attribute_name = 'units', &
[4465]	1490	value = 'm' )
[4411]	1491
[4497]	1492	return_value = dom_def_att( nc_filename, &
	1493	variable_name = 'rotor_diameter', &
	1494	attribute_name = 'units', &
[4465]	1495	value = 'm' )
	1496
[4497]	1497	return_value = dom_def_att( nc_filename, &
	1498	variable_name = 'tower_diameter', &
	1499	attribute_name = 'units', &
[4465]	1500	value = 'm' )
	1501
[4497]	1502	return_value = dom_def_att( nc_filename, &
	1503	variable_name = 'x', &
	1504	attribute_name = 'long_name', &
[4411]	1505	value = 'x location of rotor center' )
	1506
[4497]	1507	return_value = dom_def_att( nc_filename, &
	1508	variable_name = 'y', &
	1509	attribute_name = 'long_name', &
	1510	value = 'y location of rotor center' )
[4411]	1511
[4497]	1512	return_value = dom_def_att( nc_filename, &
	1513	variable_name = 'z', &
	1514	attribute_name = 'long_name', &
	1515	value = 'z location of rotor center' )
[4465]	1516
[4497]	1517	return_value = dom_def_att( nc_filename, &
	1518	variable_name = 'turbine_name', &
	1519	attribute_name = 'long_name', &
[4465]	1520	value = 'turbine name')
	1521
[4497]	1522	return_value = dom_def_att( nc_filename, &
	1523	variable_name = 'rotor_diameter', &
	1524	attribute_name = 'long_name', &
[4465]	1525	value = 'rotor diameter')
	1526
[4497]	1527	return_value = dom_def_att( nc_filename, &
	1528	variable_name = 'tower_diameter', &
	1529	attribute_name = 'long_name', &
[4465]	1530	value = 'tower diameter')
	1531
[4497]	1532	return_value = dom_def_att( nc_filename, &
	1533	variable_name = 'time', &
	1534	attribute_name = 'standard_name', &
[4411]	1535	value = 'time')
	1536
[4497]	1537	return_value = dom_def_att( nc_filename, &
	1538	variable_name = 'time', &
	1539	attribute_name = 'axis', &
[4411]	1540	value = 'T')
	1541
[4497]	1542	return_value = dom_def_att( nc_filename, &
	1543	variable_name = 'x', &
	1544	attribute_name = 'axis', &
[4411]	1545	value = 'X' )
	1546
[4497]	1547	return_value = dom_def_att( nc_filename, &
	1548	variable_name = 'y', &
	1549	attribute_name = 'axis', &
[4465]	1550	value = 'Y' )
[4411]	1551
[4497]	1552
[4465]	1553	variable_name = 'generator_power'
[4497]	1554	return_value = dom_def_var( nc_filename, &
	1555	variable_name = variable_name, &
	1556	dimension_names = (/ 'turbine', 'time ' /), &
[4465]	1557	output_type = 'real32' )
[4411]	1558
[4497]	1559	return_value = dom_def_att( nc_filename, &
	1560	variable_name = variable_name, &
	1561	attribute_name = 'units', &
	1562	value = 'W' )
[4465]	1563
[4411]	1564	variable_name = 'generator_speed'
[4497]	1565	return_value = dom_def_var( nc_filename, &
	1566	variable_name = variable_name, &
	1567	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1568	output_type = 'real32' )
[4465]	1569
[4497]	1570	return_value = dom_def_att( nc_filename, &
	1571	variable_name = variable_name, &
	1572	attribute_name = 'units', &
[4411]	1573	value = 'rad/s' )
[4465]	1574
[4411]	1575	variable_name = 'generator_torque'
[4497]	1576	return_value = dom_def_var( nc_filename, &
	1577	variable_name = variable_name, &
	1578	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1579	output_type = 'real32' )
[4497]	1580
	1581	return_value = dom_def_att( nc_filename, &
	1582	variable_name = variable_name, &
	1583	attribute_name = 'units', &
[4411]	1584	value = 'Nm' )
	1585
	1586	variable_name = 'pitch_angle'
[4497]	1587	return_value = dom_def_var( nc_filename, &
	1588	variable_name = variable_name, &
	1589	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1590	output_type = 'real32' )
[4465]	1591
[4497]	1592	return_value = dom_def_att( nc_filename, &
	1593	variable_name = variable_name, &
	1594	attribute_name = 'units', &
[4465]	1595	value = 'degrees' )
	1596
	1597	variable_name = 'rotor_power'
[4497]	1598	return_value = dom_def_var( nc_filename, &
	1599	variable_name = variable_name, &
	1600	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1601	output_type = 'real32' )
[4465]	1602
[4497]	1603	return_value = dom_def_att( nc_filename, &
	1604	variable_name = variable_name, &
	1605	attribute_name = 'units', &
	1606	value = 'W' )
	1607
[4465]	1608	variable_name = 'rotor_speed'
[4497]	1609	return_value = dom_def_var( nc_filename, &
	1610	variable_name = variable_name, &
	1611	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1612	output_type = 'real32' )
[4497]	1613
	1614	return_value = dom_def_att( nc_filename, &
	1615	variable_name = variable_name, &
	1616	attribute_name = 'units', &
[4465]	1617	value = 'rad/s' )
	1618
	1619	variable_name = 'rotor_thrust'
[4497]	1620	return_value = dom_def_var( nc_filename, &
	1621	variable_name = variable_name, &
	1622	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1623	output_type = 'real32' )
[4465]	1624
[4497]	1625	return_value = dom_def_att( nc_filename, &
	1626	variable_name = variable_name, &
	1627	attribute_name = 'units', &
	1628	value = 'N' )
	1629
[4465]	1630	variable_name = 'rotor_torque'
[4497]	1631	return_value = dom_def_var( nc_filename, &
	1632	variable_name = variable_name, &
	1633	dimension_names = (/ 'turbine', 'time ' /), &
[4465]	1634	output_type = 'real32' )
[4497]	1635
	1636	return_value = dom_def_att( nc_filename, &
	1637	variable_name = variable_name, &
	1638	attribute_name = 'units', &
[4465]	1639	value = 'Nm' )
	1640
[4411]	1641	variable_name = 'wind_direction'
[4497]	1642	return_value = dom_def_var( nc_filename, &
	1643	variable_name = variable_name, &
	1644	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1645	output_type = 'real32' )
[4465]	1646
[4497]	1647	return_value = dom_def_att( nc_filename, &
	1648	variable_name = variable_name, &
	1649	attribute_name = 'units', &
	1650	value = 'degrees' )
	1651
[4411]	1652	variable_name = 'yaw_angle'
[4497]	1653	return_value = dom_def_var( nc_filename, &
	1654	variable_name = variable_name, &
	1655	dimension_names = (/ 'turbine', 'time ' /), &
[4411]	1656	output_type = 'real32' )
[4423]	1657
[4497]	1658	return_value = dom_def_att( nc_filename, &
	1659	variable_name = variable_name, &
	1660	attribute_name = 'units', &
	1661	value = 'degrees' )
	1662
[4465]	1663	ENDIF
[4423]	1664	END SUBROUTINE
[4497]	1665
	1666	!--------------------------------------------------------------------------------------------------!
[1864]	1667	! Description:
	1668	! ------------
[4497]	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	!--------------------------------------------------------------------------------------------------!
[1912]	1672	!
[4497]	1673	SUBROUTINE wtm_read_blade_tables
[1819]	1674
[1839]	1675
[4497]	1676	IMPLICIT NONE
[1864]	1677
[4497]	1678	CHARACTER(200) :: chmess !< Read in string
[1839]	1679
[4497]	1680	INTEGER(iwp) :: ii !< running index
	1681	INTEGER(iwp) :: jj !< running index
[1839]	1682
[4497]	1683	INTEGER(iwp) :: ierrn !<
[1839]	1684
[4497]	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 !<
[1864]	1695
	1696
[4497]	1697	REAL(wp) :: alpha_attack_i !<
	1698	REAL(wp) :: weight_a !<
	1699	REAL(wp) :: weight_b !<
[1864]	1700
[4497]	1701	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint1 !<
	1702	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ttoint2 !<
[1864]	1703
[4497]	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
[1864]	1709
[4497]	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 !<
[1864]	1715
[4497]	1716	ALLOCATE ( read_cl_cd(1:2 * n_airfoils + 1) )
[1864]	1717
[4497]	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 )
[1864]	1722
[4497]	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
[1864]	1727	!
[4497]	1728	!-- Read distribution table:
	1729	dlen = 0
[1864]	1730
[4497]	1731	READ ( 90, '(3/)' )
[1864]	1732
[4497]	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
[1864]	1738
[4497]	1739	ALLOCATE( trad1(1:dlen), trad2(1:dlen), ttoint1(1:dlen), ttoint2(1:dlen) )
[1864]	1740
[4497]	1741	DO jj = 1, dlen+1
	1742	BACKSPACE ( 90, IOSTAT=ierrn )
	1743	ENDDO
[1819]	1744
[4497]	1745	DO jj = 1, dlen
	1746	READ ( 90, * ) trad1(jj), trad2(jj), ttoint1(jj), ttoint2(jj)
	1747	ENDDO
	1748
[1864]	1749	!
[4497]	1750	!-- Read layout table:
	1751	dlen = 0
[1819]	1752
[4497]	1753	READ ( 90, '(3/)')
[1864]	1754
[4497]	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
[1864]	1760
[4497]	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
[1864]	1768
[4497]	1769	!
	1770	!-- Read tables (turb_cl(alpha),turb_cd(alpha) for the different profiles:
	1771	dlen = 0
[1864]	1772
[4497]	1773	READ ( 90, '(3/)' )
[1864]	1774
[4497]	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
[1864]	1780
[4497]	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)
[1839]	1794	ENDDO
[1864]	1795
[4497]	1796	ENDDO
[1864]	1797
[4497]	1798	dlenbl = dlen
[1819]	1799
[4497]	1800	CLOSE ( 90 )
	1801
[1864]	1802	!
[4497]	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) )
[1839]	1812
[4497]	1813	radres = INT( rotor_radius(1) * 10.0_wp ) + 1_iwp
	1814	dlenbl_int = INT( 360.0_wp / accu_cl_cd_tab ) + 1_iwp
[1819]	1815
[4497]	1816	ALLOCATE( turb_cl_tab(1:dlenbl_int,1:radres) )
	1817	ALLOCATE( turb_cd_tab(1:dlenbl_int,1:radres) )
[1839]	1818
[4497]	1819	DO iir = 1, radres ! loop over radius
[1864]	1820
[4497]	1821	cur_r = ( iir - 1_iwp ) * 0.1_wp
[2836]	1822	!
[4497]	1823	!-- Find position in table 1:
	1824	lct = MINLOC( ABS( trad1 - cur_r ) )
	1825	! lct(1) = lct(1)
[1864]	1826
[4497]	1827	IF ( ( trad1(lct(1)) - cur_r ) > 0.0 ) THEN
	1828	lct(1) = lct(1) - 1
	1829	ENDIF
[1864]	1830
[4497]	1831	trow = lct(1)
[2836]	1832	!
[4497]	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)
[1864]	1838
[4497]	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
[2836]	1843
[4497]	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
[2836]	1849
[4497]	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)
[2836]	1858	!
[4497]	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
[2836]	1863
[4497]	1864	DO iialpha = 2, dlenbl_int ! loop over angles
[1819]	1865
[4497]	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
[1864]	1873	!
[4497]	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) )
[3182]	1892
[4497]	1893	ENDDO ! end loop over angles of attack
[1819]	1894
[4497]	1895	ENDDO ! end loop over radius
[2836]	1896
[1819]	1897
[4497]	1898	END SUBROUTINE wtm_read_blade_tables
[1819]	1899
[4497]	1900
	1901	!--------------------------------------------------------------------------------------------------!
[1864]	1902	! Description:
	1903	! ------------
[4497]	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 )
[1819]	1908
[1864]	1909
[4497]	1910	IMPLICIT NONE
[1864]	1911
[4497]	1912	INTEGER(iwp) :: inot
[1864]	1913	!
[4497]	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
[1864]	1918
[4497]	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
[1864]	1922
[4497]	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
[1864]	1927	!
[4497]	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 )
[1819]	1954
[1864]	1955	!
[4497]	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 )
[1839]	1961
[4497]	1962	END SUBROUTINE wtm_rotate_rotor
[1839]	1963
[4497]	1964
	1965	!--------------------------------------------------------------------------------------------------!
[1819]	1966	! Description:
	1967	! ------------
[1839]	1968	!> Calculation of the forces generated by the wind turbine
[4497]	1969	!--------------------------------------------------------------------------------------------------!
	1970	SUBROUTINE wtm_forces
[1819]	1971
[1864]	1972
[4497]	1973	IMPLICIT NONE
[1819]	1974
	1975
[4497]	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 !<
[4467]	1982
[4497]	1983	REAL(wp) :: flag !< flag to mask topography grid points
	1984	REAL(wp) :: sin_rot, cos_rot !<
	1985	REAL(wp) :: sin_yaw, cos_yaw !<
[4467]	1986
[4497]	1987	REAL(wp) :: aa, bb, cc, dd !< interpolation distances
	1988	REAL(wp) :: gg !< interpolation volume var
[4467]	1989
[4497]	1990	REAL(wp) :: dist_u_3d, dist_v_3d, dist_w_3d !< smearing distances
	1991
	1992
[1839]	1993	!
[4497]	1994	!-- Variables for pitch control:
	1995	INTEGER(iwp), DIMENSION(1) :: lct = 0
[1839]	1996
[4497]	1997	LOGICAL :: pitch_sw = .FALSE.
[1819]	1998
[4497]	1999	REAL(wp), DIMENSION(1) :: rad_d = 0.0_wp
[1819]	2000
[4497]	2001	REAL(wp) :: tl_factor !< factor for tip loss correction
[1819]	2002
	2003
[4497]	2004	CALL cpu_log( log_point_s(61), 'wtm_forces', 'start' )
[1819]	2005
[4497]	2006
	2007	IF ( time_since_reference_point >= time_turbine_on ) THEN
	2008
[1864]	2009	!
[4497]	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
[1819]	2019	!
[4497]	2020	!-- Loop over number of turbines:
	2021	DO inot = 1, n_turbines
[1819]	2022
[4497]	2023	cos_yaw = COS(yaw_angle(inot))
	2024	sin_yaw = SIN(yaw_angle(inot))
[1819]	2025	!
[4497]	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)
[3832]	2031
[4497]	2032	thrust_seg(:) = 0.0_wp
	2033	torque_seg_y(:) = 0.0_wp
	2034	torque_seg_z(:) = 0.0_wp
[1819]	2035	!
[4497]	2036	!-- Determine distance between each ring (center) and the hub:
	2037	cur_r = (ring - 0.5_wp) * delta_r(inot)
[1819]	2038
	2039	!
[4497]	2040	!-- Loop over segments of each ring of each turbine:
	2041	DO rseg = 1, nsegs(ring,inot)
[1819]	2042	!
[4497]	2043	!-- !----------------------------------------------------------------------------------!
	2044	!-- !-- Determine coordinates of the ring segments --!
	2045	!-- !----------------------------------------------------------------------------------!
[1819]	2046	!
[4497]	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 )
[1819]	2052
[4497]	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
[1819]	2057
[4497]	2058	rote = MATMUL( rot_coord_trans(inot,:,:), re )
[1819]	2059	!
[4497]	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)
[1819]	2064
[4497]	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	!-- !----------------------------------------------------------------------------------!
[1819]	2070
[4497]	2071	u_int(inot,ring,rseg) = 0.0_wp
	2072	u_int_1_l(inot,ring,rseg) = 0.0_wp
[1819]	2073
[4497]	2074	v_int(inot,ring,rseg) = 0.0_wp
	2075	v_int_1_l(inot,ring,rseg) = 0.0_wp
[1819]	2076
[4497]	2077	w_int(inot,ring,rseg) = 0.0_wp
	2078	w_int_1_l(inot,ring,rseg) = 0.0_wp
[1819]	2079
	2080	!
[4497]	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)
[1819]	2085	!
[4497]	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
[1819]	2089
[4497]	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
[1819]	2099
[4497]	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
[1819]	2102
[4497]	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
[1819]	2105
[4497]	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 )
[1819]	2108
	2109	ELSE
	2110	u_int_1_l(inot,ring,rseg) = 0.0_wp
	2111	ENDIF
	2112
[4497]	2113	ELSE
	2114	u_int_1_l(inot,ring,rseg) = 0.0_wp
	2115	ENDIF
[1819]	2116
[4497]	2117
[1819]	2118	!
[4497]	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)
[1819]	2123	!
[4497]	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
[1819]	2127
[4497]	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
[1819]	2137
[4497]	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
[1819]	2140
[4497]	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
[1819]	2143
[4497]	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 )
[1819]	2146
	2147	ELSE
	2148	v_int_1_l(inot,ring,rseg) = 0.0_wp
	2149	ENDIF
	2150
[4497]	2151	ELSE
	2152	v_int_1_l(inot,ring,rseg) = 0.0_wp
	2153	ENDIF
[1819]	2154
[4497]	2155
[1819]	2156	!
[4497]	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)
[1819]	2161	!
[4497]	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
[1819]	2165
[4497]	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
[1819]	2175
[4497]	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
[1819]	2178
[4497]	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
[1819]	2181
[4497]	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 )
[1819]	2184	ELSE
	2185	w_int_1_l(inot,ring,rseg) = 0.0_wp
	2186	ENDIF
	2187
[4497]	2188	ELSE
	2189	w_int_1_l(inot,ring,rseg) = 0.0_wp
	2190	ENDIF
	2191
[1819]	2192	ENDDO
	2193
	2194	ENDDO
[4497]	2195	!$OMP END PARALLEL
[1819]	2196
[4497]	2197	ENDDO
	2198
[1819]	2199	!
[4497]	2200	!-- Exchange between PEs (information required on each PE):
[1819]	2201	#if defined( __parallel )
[4497]	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 )
[1819]	2208	#else
[4497]	2209	u_int = u_int_1_l
	2210	v_int = v_int_1_l
	2211	w_int = w_int_1_l
[1819]	2212	#endif
	2213
	2214
	2215	!
[4497]	2216	!-- Loop over number of turbines:
	2217	DO inot = 1, n_turbines
[1912]	2218	pit_loop: DO
[1819]	2219
[4497]	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 )
[1819]	2230	ENDIF
[4497]	2231	ENDIF
[1819]	2232
[1839]	2233	!
[4497]	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)
[1819]	2242	!
[4497]	2243	!-- Determine distance between each ring (center) and the hub:
	2244	cur_r = ( ring - 0.5_wp ) * delta_r(inot)
[1839]	2245	!
[4497]	2246	!-- Loop over segments of each ring of each turbine:
	2247	DO rseg = 1, nsegs(ring,inot)
[1819]	2248	!
[4497]	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 )
[1819]	2254	!
[4497]	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
[1819]	2259
[4497]	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
[1819]	2264
	2265	!
[4497]	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
[1819]	2272	!
[4497]	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 )
[1819]	2278	!
[4497]	2279	!-- Coordinates of the single segments (center points):
	2280	rbx(ring,rseg) = hub_x(inot) + cur_r * rote(1)
[1819]	2281
[4497]	2282	rby(ring,rseg) = hub_y(inot) + cur_r * rote(2)
[1819]	2283
[4497]	2284	rbz(ring,rseg) = hub_z(inot) + cur_r * rote(3)
[1819]	2285
	2286	!
[4497]	2287	!-- !----------------------------------------------------------------------------------!
	2288	!-- !-- Calculation of various angles and relative velocities --!
	2289	!-- !----------------------------------------------------------------------------------!
[1819]	2290	!
[4497]	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.
[1864]	2296	!
[4497]	2297	!-- Projection of the xy-velocities relative to the rotor area:
[1864]	2298	!
[4497]	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)
[1819]	2302	!
[4497]	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)
[1819]	2307	!
	2308
[4497]	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) ) )
[1819]	2312
	2313	!
[4497]	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)
[1819]	2318
[4497]	2319	IF ( cur_r == 0.0_wp ) THEN
	2320	alpha_attack(rseg) = 0.0_wp
[1819]	2321
[4497]	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
[1819]	2332	!
[4497]	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 )
[1819]	2337	!
[4497]	2338	!-- Substraction of the local pitch angle to obtain the local angle of attack:
	2339	alpha_attack(rseg) = phi_rel(rseg) - alpha_attack(rseg)
[1819]	2340	!
[4497]	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 ) )
[1819]	2344	!
[4497]	2345	!-- Correct with collective pitch angle:
	2346	alpha_attack(rseg) = alpha_attack(rseg) - pitch_angle(inot)
[1819]	2347
	2348	!
[4497]	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 )
[1819]	2351
	2352	!
[4497]	2353	!-- !----------------------------------------------------------------------------------!
	2354	!-- !-- Interpolation of chord as well as lift and drag --!
	2355	!-- !-- coefficients from tabulated values --!
	2356	!-- !----------------------------------------------------------------------------------!
[1819]	2357
	2358	!
[4497]	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
[1819]	2363
[4497]	2364	ELSE IF ( cur_r >= lrd( size(crd) ) ) THEN
	2365	chord(rseg) = ( crd( size(crd) ) + ard( size(crd) - 1 ) ) / 2.0_wp
[1819]	2366
[4497]	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
[1819]	2374	!
[4497]	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
[1819]	2380	!
[4497]	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 )
[1819]	2384	!
[4497]	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)
[1819]	2389
	2390	!
[4497]	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 )
[1819]	2393
[4497]	2394	IF ( tip_loss_correction ) THEN
[4465]	2395	!
[4497]	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
[4465]	2402
[4497]	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
[4465]	2408
[4497]	2409	turb_cl(rseg) = tl_factor * turb_cl(rseg)
	2410
	2411	ENDIF
[1819]	2412	!
[4497]	2413	!-- !----------------------------------------------------------------------------------!
	2414	!-- !-- Calculation of the forces --!
	2415	!-- !----------------------------------------------------------------------------------!
[1819]	2416
	2417	!
[4497]	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)
[1819]	2421
	2422	!
[4497]	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) ) )
[1819]	2426
	2427	!
[4497]	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) ) )
[1819]	2433	!
[4497]	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
[1819]	2438
[1912]	2439	!
[4497]	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)
[1912]	2443
[4465]	2444
[4497]	2445	torque_total(inot) = torque_total(inot) + (torque_seg * cur_r)
	2446	!$OMP END CRITICAL
[1819]	2447
[4497]	2448	ENDDO !-- end of loop over ring segments
[1819]	2449
	2450	!
[4497]	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(:)
[1819]	2455
	2456
[4497]	2457	ENDDO !-- end of loop over rings
	2458	!$OMP END PARALLEL
[1819]	2459
	2460
[4497]	2461	CALL cpu_log( log_point_s(62), 'wtm_controller', 'start' )
[1819]	2462
[4465]	2463
[4497]	2464	IF ( speed_control ) THEN
[1912]	2465	!
[4497]	2466	!-- Calculation of the current generator speed for rotor speed control:
[4465]	2467
[4497]	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
[1819]	2474
[1912]	2475	!
[4497]	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 )
[4465]	2478
[4497]	2479	ENDIF
[4465]	2480
[4497]	2481	IF ( pitch_control ) THEN
[1819]	2482
[1912]	2483	!
[4497]	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.
[1912]	2491	!
[4497]	2492	!-- Go back to beginning of pit_loop:
	2493	CYCLE pit_loop
	2494	ENDIF
[4465]	2495
[1912]	2496	!
[4497]	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
[1819]	2503
[1864]	2504
[1819]	2505	!
[4497]	2506	!-- Call the rotor speed controller:
	2507	IF ( speed_control ) THEN
[1819]	2508	!
[4497]	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 )
[1912]	2513	ENDIF
[1819]	2514	ENDIF
	2515
[4497]	2516	ENDIF
[1819]	2517
	2518
[4497]	2519	CALL cpu_log( log_point_s(62), 'wtm_controller', 'stop' )
[1819]	2520
[4497]	2521	CALL cpu_log( log_point_s(63), 'wtm_smearing', 'start' )
[1819]	2522
[4497]	2523
	2524	!-- !----------------------------------------------------------------------------------------!
	2525	!-- !-- Regularization kernel --!
	2526	!-- !-- Smearing of the forces and interpolation to cartesian grid --!
	2527	!-- !----------------------------------------------------------------------------------------!
[1819]	2528	!
[4497]	2529	!-- The aerodynamic blade forces need to be distributed smoothly on several mesh points in
	2530	!-- order to avoid singular behaviour.
[1819]	2531	!
[4497]	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.
[1819]	2535	!
[4497]	2536	!-- To save computing time, apply smearing only for the relevant part of the model domain:
[1819]	2537	!
	2538	!--
[4497]	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 )
[1819]	2542
[4497]	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)
[1819]	2552	!
[4497]	2553	!-- Predetermine flag to mask topography:
	2554	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[4465]	2555
[2232]	2556	!
[4497]	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
[1819]	2562
[4497]	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
[1819]	2571	!
[4497]	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
[1819]	2581
[4497]	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
[4465]	2587
[4497]	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
[1819]	2597	!
[4497]	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
[4465]	2602
[4497]	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
[4465]	2606
[4497]	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
[1819]	2610
[4497]	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
[1819]	2615
[4497]	2616	CALL cpu_log( log_point_s(63), 'wtm_smearing', 'stop' )
[4465]	2617
[4497]	2618	ENDDO !-- end of loop over turbines
[1819]	2619
[4465]	2620
[4497]	2621	IF ( yaw_control ) THEN
[1819]	2622	!
[4497]	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) )
[4465]	2629
[4497]	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
[1819]	2636
	2637	!
[4497]	2638	!-- Calculate the inflow wind speed:
[1819]	2639	!--
[4497]	2640	!-- Loop over number of turbines:
	2641	DO inot = 1, n_turbines
[1819]	2642	!
[4497]	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)) ) )
[1864]	2652
[4497]	2653	CALL wtm_yawcontrol( inot )
[1864]	2654
[4497]	2655	yaw_angle_l(inot) = yaw_angle(inot)
[1864]	2656
[1819]	2657	ENDIF
[4497]	2658	ENDIF
[4465]	2659
[4497]	2660	ENDDO ! end of loop over turbines
[1864]	2661
	2662	!
[4497]	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 )
[1929]	2671	#else
[4497]	2672	u_inflow = u_inflow_l
	2673	wdir = wdir_l
	2674	yaw_angle = yaw_angle_l
[4465]	2675
[1929]	2676	#endif
[4497]	2677	DO inot = 1, n_turbines
[4465]	2678	!
[4497]	2679	!-- Update rotor orientation:
	2680	CALL wtm_rotate_rotor( inot )
[1819]	2681
[4497]	2682	ENDDO ! End of loop over turbines
[4465]	2683
[4497]	2684	ENDIF ! end of yaw control
[4465]	2685
[4497]	2686	IF ( speed_control ) THEN
[1864]	2687	!
[4497]	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 )
[1929]	2698	#else
[4497]	2699	torque_gen_old = torque_gen
	2700	rotor_speed = rotor_speed_l
	2701	generator_speed_f_old = generator_speed_f
[1929]	2702	#endif
[4465]	2703
	2704	ENDIF
[1819]	2705
[4497]	2706	ENDIF
[1819]	2707
[4497]	2708	CALL cpu_log( log_point_s(61), 'wtm_forces', 'stop' )
[1819]	2709
	2710
[4497]	2711	END SUBROUTINE wtm_forces
[4465]	2712
[4497]	2713
	2714	!--------------------------------------------------------------------------------------------------!
[1839]	2715	! Description:
	2716	! ------------
	2717	!> Yaw controller for the wind turbine model
[4497]	2718	!--------------------------------------------------------------------------------------------------!
	2719	SUBROUTINE wtm_yawcontrol( inot )
[4465]	2720
[4497]	2721	USE kinds
[4465]	2722
[4497]	2723	IMPLICIT NONE
[4465]	2724
[4497]	2725	INTEGER(iwp) :: inot
	2726	INTEGER(iwp) :: i_wd_30
[1819]	2727
[4497]	2728	REAL(wp) :: missal
[1819]	2729
[4497]	2730	i_wd_30 = 0_iwp
	2731
[1864]	2732	!
[4497]	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!
[1843]	2741	!
[4497]	2742	!-- Check if turbine is not yawing:
	2743	IF ( .NOT. doyaw(inot) ) THEN
[1843]	2744	!
[4497]	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
[1819]	2752
[4497]	2753	missal = SUM( wd30_l ) / SIZE( wd30_l ) - yaw_angle(inot)
[1819]	2754	!
[4497]	2755	!-- Check if turbine is missaligned by more than yaw_misalignment_max:
	2756	IF ( ABS( missal ) > yaw_misalignment_max ) THEN
[1843]	2757	!
[4497]	2758	!-- Check in which direction to yaw:
	2759	yawdir(inot) = SIGN( 1.0_wp, missal )
[1819]	2760	!
[4497]	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
[1819]	2765	ENDIF
[4497]	2766	ENDIF
[4465]	2767
[4497]	2768	wd30(inot,:) = wd30_l
[1819]	2769
[4465]	2770	!
[4497]	2771	!-- If turbine is already yawing:
	2772	!-- Initialize 2 s running mean and yaw until the missalignment is smaller than
	2773	!-- yaw_misalignment_min
[1819]	2774
[4497]	2775	ELSE
[1819]	2776	!
[4497]	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)
[4465]	2782	!
[4497]	2783	!-- Check if array is full ( no more dummies ):
	2784	IF ( .NOT. ANY( wd2_l == 999.0_wp ) ) THEN
[1864]	2785	!
[4497]	2786	!-- Calculate missalignment of turbine:
	2787	missal = SUM( wd2_l - yaw_angle(inot) ) / SIZE( wd2_l )
[1864]	2788	!
[4497]	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
[1864]	2793	!
[4497]	2794	!-- Continue yawing:
	2795	yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d
	2796
[1819]	2797	ELSE
[1864]	2798	!
[4497]	2799	!-- Stop yawing:
	2800	doyaw(inot) = .FALSE.
	2801	wd2_l = 999.0_wp ! fill with dummies again
[1819]	2802	ENDIF
[4497]	2803	ELSE
	2804	!
	2805	!-- Continue yawing:
	2806	yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d
	2807	ENDIF
[4465]	2808
[4497]	2809	wd2(inot,:) = wd2_l
[4465]	2810
[4497]	2811	ENDIF
[4465]	2812
[4497]	2813	END SUBROUTINE wtm_yawcontrol
[1819]	2814
[1864]	2815
[4497]	2816	!--------------------------------------------------------------------------------------------------!
[1819]	2817	! Description:
	2818	! ------------
[1864]	2819	!> Initialization of the speed control
[4497]	2820	!--------------------------------------------------------------------------------------------------!
	2821	SUBROUTINE wtm_init_speed_control
[1864]	2822
	2823
[4497]	2824	IMPLICIT NONE
[1864]	2825
	2826	!
[4497]	2827	!-- If speed control is set, remaining variables and control_parameters for the control algorithm
	2828	!-- are calculated:
[1864]	2829	!
[4497]	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 )
[1864]	2833	!
[4497]	2834	!-- Calculate upper limit of slipage region:
	2835	vs_sysp = generator_speed_rated / 1.1_wp
[1864]	2836	!
[4497]	2837	!-- Calculate slope of slipage region:
	2838	region_25_slope = ( generator_power_rated / generator_speed_rated ) / &
	2839	( generator_speed_rated - vs_sysp )
[1864]	2840	!
[4497]	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 )
[1864]	2845	!
[4497]	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
[1864]	2853
[4497]	2854	END SUBROUTINE wtm_init_speed_control
[1864]	2855
[4497]	2856
	2857	!--------------------------------------------------------------------------------------------------!
[1864]	2858	! Description:
	2859	! ------------
	2860	!> Simple controller for the regulation of the rotor speed
[4497]	2861	!--------------------------------------------------------------------------------------------------!
	2862	SUBROUTINE wtm_speed_control( inot )
[1864]	2863
	2864
[4497]	2865	IMPLICIT NONE
[1864]	2866
[4497]	2867	INTEGER(iwp):: inot
[1864]	2868
[4465]	2869
[1864]	2870	!
[4497]	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"
[1864]	2873
	2874	!
[4497]	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)
[1864]	2879
[4497]	2880	IF ( generator_speed_f(inot) <= region_15_min ) THEN
[4465]	2881	!
[4497]	2882	!-- Region 1: Generator torque is set to zero to accelerate the rotor:
	2883	torque_gen(inot) = 0
[4465]	2884
[4497]	2885	ELSEIF ( generator_speed_f(inot) <= region_2_min ) THEN
[4465]	2886	!
[4497]	2887	!-- Region 1.5: Generator torque is increasing linearly with rotor speed:
	2888	torque_gen(inot) = slope15 * ( generator_speed_f(inot) - region_15_min )
[4465]	2889
[4497]	2890	ELSEIF ( generator_speed_f(inot) <= region_25_min ) THEN
[1864]	2891	!
[4497]	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)
[4465]	2895
[4497]	2896	ELSEIF ( generator_speed_f(inot) < generator_speed_rated ) THEN
[4465]	2897	!
[4497]	2898	!-- Region 2.5: Slipage region between 2 and 3:
	2899	torque_gen(inot) = region_25_slope * ( generator_speed_f(inot) - vs_sysp )
[4465]	2900
[4497]	2901	ELSE
[4465]	2902	!
[4497]	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)
[4465]	2906
[4497]	2907	ENDIF
[4465]	2908	!
[4497]	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 )
[4465]	2912	!
[4497]	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 )
[1864]	2916	!
[4497]	2917	!-- Overwrite values for next timestep:
	2918	rotor_speed_l(inot) = generator_speed(inot) / gear_ratio
[1864]	2919
[4465]	2920
[4497]	2921	END SUBROUTINE wtm_speed_control
[1864]	2922
	2923
[4497]	2924	!--------------------------------------------------------------------------------------------------!
[1864]	2925	! Description:
	2926	! ------------
[4497]	2927	!> Application of the additional forces generated by the wind turbine on the flow components
	2928	!> (tendency terms)
[1839]	2929	!> Call for all grid points
[4497]	2930	!--------------------------------------------------------------------------------------------------!
	2931	SUBROUTINE wtm_actions( location )
[1819]	2932
	2933
[4497]	2934	CHARACTER(LEN=*) :: location !<
[3875]	2935
[4497]	2936	INTEGER(iwp) :: i !< running index
	2937	INTEGER(iwp) :: j !< running index
	2938	INTEGER(iwp) :: k !< running index
[1819]	2939
	2940
[4497]	2941	SELECT CASE ( location )
[1819]	2942
[4497]	2943	CASE ( 'before_timestep' )
[3875]	2944
[4497]	2945	CALL wtm_forces
[3875]	2946
[4497]	2947	CASE ( 'u-tendency' )
[1819]	2948	!
[4497]	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 )
[1819]	2955	!
[4497]	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 ) )
[1819]	2961	ENDDO
	2962	ENDDO
[4497]	2963	ENDDO
	2964	ENDIF
[1819]	2965
[4497]	2966	CASE ( 'v-tendency' )
[1819]	2967	!
[4497]	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 ) )
[1819]	2977	ENDDO
	2978	ENDDO
[4497]	2979	ENDDO
	2980	ENDIF
[1819]	2981
[4497]	2982	CASE ( 'w-tendency' )
[1819]	2983	!
[4497]	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 ) )
[1819]	2991	ENDDO
	2992	ENDDO
[4497]	2993	ENDDO
	2994	ENDIF
[1819]	2995
	2996
[4497]	2997	CASE DEFAULT
	2998	CONTINUE
[1819]	2999
[4497]	3000	END SELECT
[1819]	3001
	3002
[4497]	3003	END SUBROUTINE wtm_actions
[1819]	3004
	3005
[4497]	3006	!--------------------------------------------------------------------------------------------------!
[1819]	3007	! Description:
	3008	! ------------
[4497]	3009	!> Application of the additional forces generated by the wind turbine on the flow components
	3010	!> (tendency terms)
[1839]	3011	!> Call for grid point i,j
[4497]	3012	!--------------------------------------------------------------------------------------------------!
	3013	SUBROUTINE wtm_actions_ij( i, j, location )
[1819]	3014
	3015
[4497]	3016	CHARACTER (LEN=*) :: location !<
[1819]	3017
[4497]	3018	INTEGER(iwp) :: i !< running index
	3019	INTEGER(iwp) :: j !< running index
	3020	INTEGER(iwp) :: k !< running index
[1819]	3021
[4497]	3022	SELECT CASE ( location )
[3875]	3023
[4497]	3024	CASE ( 'before_timestep' )
[3875]	3025
[4497]	3026	CALL wtm_forces
	3027
	3028	CASE ( 'u-tendency' )
[1819]	3029	!
[4497]	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)
[1819]	3033	!
[4497]	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
[1819]	3041
[4497]	3042	CASE ( 'v-tendency' )
[1819]	3043	!
[4497]	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
[1819]	3053
[4497]	3054	CASE ( 'w-tendency' )
[1819]	3055	!
[4497]	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
[1819]	3063
	3064
[4497]	3065	CASE DEFAULT
	3066	CONTINUE
[1819]	3067
[4497]	3068	END SELECT
[1819]	3069
	3070
[4497]	3071	END SUBROUTINE wtm_actions_ij
[1819]	3072
[4411]	3073
[4497]	3074	SUBROUTINE wtm_data_output
[4411]	3075
[4465]	3076
[4497]	3077	INTEGER(iwp) :: return_value !< returned status value of called function
	3078	INTEGER(iwp) :: t_ind = 0 !< time index
[4465]	3079
[4497]	3080	IF ( myid == 0 ) THEN
[4465]	3081
[4497]	3082	!
	3083	!-- At the first call of this routine write the spatial coordinates:
	3084	IF ( .NOT. initial_write_coordinates ) THEN
[4447]	3085	ALLOCATE ( output_values_1d_target(1:n_turbines) )
[4497]	3086	output_values_1d_target = hub_x(1:n_turbines)
[4423]	3087	output_values_1d_pointer => output_values_1d_target
[4497]	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/) )
[4411]	3093
[4497]	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/) )
[4411]	3101
[4497]	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/) )
[4411]	3109
[4497]	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/) )
[4411]	3117
[4497]	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/) )
[4411]	3125
[4497]	3126	initial_write_coordinates = .TRUE.
	3127	DEALLOCATE ( output_values_1d_target )
	3128	ENDIF
[4411]	3129
[4497]	3130	t_ind = t_ind + 1
[4411]	3131
[4497]	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
[4411]	3135
[4497]	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 /) )
[4411]	3141
[4497]	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 /) )
[4411]	3149
[4497]	3150	output_values_1d_target = torque_gen_old(:)
	3151	output_values_1d_pointer => output_values_1d_target
[4411]	3152
[4497]	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 /) )
[4411]	3158
[4497]	3159	output_values_1d_target = torque_total(:)
	3160	output_values_1d_pointer => output_values_1d_target
[4411]	3161
[4497]	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
[1819]	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 |