Home

Context Navigation

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

Last change on this file since 3806 was 3725, checked in by raasch, 6 years ago
modifications to avoid compiler warnings about unused variables, temperton-fft: GOTO statements replaced, file re-formatted corresponding to coding standards, ssh-calls for compilations on remote systems modified to avoid output of login messages on specific systems changed again (palmbuild, reverted as before r3549), error messages for failed restarts extended (palmrun)
Property svn:keywords set to `Id`
File size: 120.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |