Documentation for wind farm parameterisation
Natalie Theeuwes & Bert van Ulft<br> contact: Natalie.Theeuwes@knmi.nl
Introduction
We've implemented the Fitch et al. (2012) scheme in HARMONIE-AROME. First implemented in CY40 and evaluated for 1 year, 2016. Results are published in this paper: (van Stratum et al., 2022) Afterwards it has been moved to CY43 and run for 3 year (2019-2021) and evaluated.
Implementation in CY46
Switching on parameterisation
To switch the parameterisation on go to ecf/config_exp.h and set WINDFARM="yes". This will make the necessary switch to the namelist using the scr/forecast_model_settings.sh
# Fitch et al. (2012) wind turbine parametrization
export LWINDFARM=.FALSE.
if [ "$WINDFARM" = yes ]; then
LWINDFARM=.TRUE.
...
fiIn nam/harmonie_namelists.pm in the namelist in the %arome=( NAMPHY=>{ the switch for the wind farm paremeterisation should be located: 'LWINDFARM' => $ENV{LWINDFARM},
Once this switch in ecf/config_exp.h the code will look for a directory called WFP_input_files in the directory of the experiment. This should include input text files for the wind farm parameterisation.
Code
The input files are read in in src/arpifs/setup/su0yomb.F90 which calls the setup of the wind farm parameterisation src/arpifs/setup/suwindfarm.F90:
!* Initialize wind farm scheme
IF (LWINDFARM) THEN
WRITE(NULOUT,*) '---- Set up wind farms --',CLINE
CALL SUWINDFARM(YDGEOMETRY)
ENDIFThe setup routine allocates the turbines on each BLK and it's properties (i.e., radius, hub height, location, type), plus it allocates the power and thrust curves for each turbine type. src/arpifs/setup/suwindfarm.F90 uses:
src/arpifs/setup/suwindfarm_rc.F90- to read the coordinate file,src/arpifs/setup/suwindfarm_rt.F90- to read the turbine type files, andsrc/arpifs/utility/skip_comments.F90- used to skip comments in the .tab files.
The wind farm parameterisation is called in src/arpifs/phys_dmn/apl_arome.F90 :
! ------------------------------------------------------------------
! 13 - CALL WIND FARM (DRAG/TKE-prod) PARAMETERISATON
! ------------------------------------------------------------------
IF (LWINDFARM) THEN
CALL WINDFARM(YDGEOMETRY%YREGEO,KBL,KLON,KLEV,&
& PAPHIM,PAPHIFM,POROG,PUM,PVM,PDT,&
& PCWFP,PWFPA,PTENDU,PTENDV,PTENDTKE)
ENDIFThe wind farm parameterisation subroutine is located in src/arpifs/phys_dmn/windfarm.F90. The parameterisation loops over each turbine in KBL and adjust the tendencies and power production.
The U, V and TKE tendencies are changed in by this call. In addition, we output two new variables:
PCWFP- accumulated wind farm power productionPWFPA- instantaneous wind farm power production
These are passed through the CFU (accumulated) and XFU (instantaneous) thingy's to the FA file output.
Input files
There are two types of input files:
wind_turbine_coordinates.tabthere should be only one file with the location of all individual turbines, the coordinates given in WG84, turbine type, radius, and hub height. This file can look like:
# lon lat type r z
6.2031 53.4072 999 8 25
5.2345 52.7387 999 8 31
5.9171 53.3762 999 8 25
5.5033 53.1913 999 8 25
5.9415 52.9295 999 8 31
5.2454 52.7327 999 8 31
6.3057 53.0353 999 8 31
4.3082 51.5722 999 8 31
5.7044 53.2521 999 8 31
5.7689 53.1922 999 8 31
6.0345 53.3543 999 8 25
5.7798 53.3212 999 8 25
5.6054 52.7066 999 8 31
5.5258 53.2502 999 8 25
5.5189 53.2468 999 8 31
5.6110 52.7242 999 9 31
4.6530 52.7281 999 8 25
5.7255 53.1632 999 9 31
4.6861 52.7434 999 8 25
5.3789 52.3055 999 9 25
5.4726 53.1554 999 8 31
6.5278 53.2488 999 8 31
6.3006 53.3322 999 8 31
5.8826 53.3058 999 8 31wind_turbine_XXX.tabEach turbine type should have their own individual file, including the thrust and power curves, with additional info e.g. the name of the type hub height and radius. Note:: Only $C_T$ and $C_P$ are read in, the remainder of the info is for reference only. Hub height and radius are read in from thewind_turbine_coordinates.tabfile. TheXXXis the turbine type number under the columntypeinwind_turbine_coordinates.tab. Awind_turbine_XXX.tabfile can look like:
# SWT-3.6-120 (z=82 m, D=120 m)
# r (m) z (m) cT_low (-) cT_high (-)
60.0000 82.0000 0.0480 0.0480
# U (m/s) cP (-) cT (-)
3.0000 0.0000 0.0000
4.0000 0.3630 0.8580
5.0000 0.4050 0.8580
6.0000 0.4240 0.8580
7.0000 0.4320 0.8570
8.0000 0.4350 0.8580
9.0000 0.4360 0.8570
10.0000 0.4200 0.8040
11.0000 0.3690 0.6070
12.0000 0.2980 0.4180
13.0000 0.2360 0.3160
14.0000 0.1890 0.2480
15.0000 0.1540 0.2010
16.0000 0.1270 0.1650
17.0000 0.1060 0.1380
18.0000 0.0890 0.1170
19.0000 0.0760 0.1000
20.0000 0.0650 0.0870
21.0000 0.0560 0.0760
22.0000 0.0490 0.0670
23.0000 0.0430 0.0600
24.0000 0.0380 0.0530
25.0000 0.0330 0.0480The `cT_low (-) cT_high (-)` indicate the thrust coefficient used
when the turbine is switched off, in either low or very high wind
speed conditions.At the moment the directory with input files is set in Env_system:
export WINDFARM_DATA_PATH=...The coordinate file is set in ecf/config_exp.h as:
# File with coordinates and properties of wind turbines
WINDFARM_COORDS_FILE=$WINDFARM_DATA_PATH/wind_turbine_coordinates_${DOMAIN}_@YYYY@.tabThe forecastmodelsettings.sh script creates links to the input files in the working directory. In research running mode the script will abort if the input files are not found. If RUNNING_MODE==operational, if the coordinate file for that DTG is not found (e.g. to avoid crashes in operations at New Year's eve), the script searches back up to 1 year for an alternative, previous file and uses that one instead. If no suitable file is found the wind farm parameterisation will be switched off (WINDFARM=no; LWINDFARM=.FALSE.). In both cases this will be logged in the severe_warnings.txt file.
Some python code to create the coordinate files from a pandas dataframe:
# Write coordinates/types to input file HARMONIE
f = open('wind_turbine_coordinates_2022.tab', 'w')
f.write('# {0:^12s} {1:^12s} {2:^12s} {3:^12s} {4:^12s}\n'.format('lon','lat','type','r','z'))
for i in range(len(df.index)):
f.write('{0:11.4f} {1:11.4f} {2:11d} {3:12d} {4:12d}\n'.format(df.lon.iloc[i],
df.lat.iloc[i], df.type.iloc[i], int(df.r.iloc[i]), int(df.z.iloc[i])))NOTE!: When creating the wind_turbine_coordinates.tab files all turbines must be within the domain. This is not done in the code. All turbines outside the domain will accumulate in the edge cells of the domain. Easiest way to do this is in the preprocessing. Like using (with turbines in a pandas dataframe):
grid = xr.open_dataset('../common_data/Harmonie_grid.nc')
hm_lon = grid['lon'].values
hm_lat = grid['lat'].values
maxlat = hm_lat.shape[0]-1
maxlon = hm_lat.shape[1]-1
# remove turbines beyond the boundaries of the HARMONIE domain
for n in np.arange(len(df.index)):
j,i = np.unravel_index(np.sqrt((hm_lon-df.Longitude.iloc[n])**2
+ (hm_lat-df.Latitude.iloc[n])**2).argmin(), hm_lon.shape)
if ( (j == 0) | (j >= maxlat) ):
# set longitude to NaN's
df.Longitude.iloc[n] = np.nan
elif ( (i == 0) | (i>=maxlon) ):
# set longitude to NaN's
df.Longitude.iloc[n] = np.nan
df.dropna(subset=['Longitude'], how='all', inplace=True)For the creation of the turbine type files wind_turbine_XXX.tab, the following python function can be used:
def write_harmonie_turbine_input(tid, name, zhub, D, u, cp, ct):
"""
Write the turbine specifications in the correct input format for HARMONIE
"""
print('Writing turbine input file for {0} - ID = {1:2d} (zhub = {2:.0f} m,
D = {3:.0f} m)'.format(name, tid, zhub, D))
f = open('wind_turbine_{0:03d}.tab'.format(tid), 'w')
f.write('# {0} (z={1:.0f} m, D={2:.0f} m)\n'.format(name, zhub, D))
f.write('# {0:^12s} {1:^12s} {2:^12s} {3:^12s}\n'.format('r (m)','z (m)',
'cT_low (-)','cT_high (-)'))
f.write('{0:11.4f} {1:11.4f} {2:11.4f} {3:11.4f}\n'.format(D/2., zhub,
ct[-1], ct[-1]))
f.write('# {0:^12s} {1:^12s} {2:^12s}\n'.format('U (m/s)','cP (-)','cT (-)'))
for i in range(u.size):
f.write('{0:11.4f} {1:11.4f} {2:11.4f}\n'.format(u[i], cp[i], ct[i]))
f.close()Output
For the conversion to GRIB1 of the two new variables, we use:
#Wind power production, accumulated
'253211' = {
table2Version = 253 ;
indicatorOfParameter = 211 ;
timeRangeIndicator = 4 ;
}
#Wind power production, instantaneous
'253211' = {
table2Version = 253 ;
indicatorOfParameter = 211 ;
timeRangeIndicator = 0 ;
}And for the conversion to GRIB2, we use:
#Wind power production, accumulated
'253211' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 39 ;
typeOfStatisticalProcessing = 1 ;
centre = 94 ;
editionNumber = 2 ;
interpretationOfNumberOfPoints = 0 ;
subCentre = 255 ;
}
#Wind power production, instantaneous
'253211' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 39 ;
centre = 94 ;
editionNumber = 2 ;
interpretationOfNumberOfPoints = 0 ;
subCentre = 255 ;
}For both GRIB 1 and GRIB 2:
Wind power production, accumulated:
- name:
'Wind power production, accumulated' - paramId:
'253211' - shortName:
'wfpower_acc' - units:
'MJ'
- name:
Wind power production, accumulated:
- name:
'Wind power production, instantaneous' - paramId:
'253211' - shortName:
'wfpower_ins' - units:
'MW'
- name: