Renegade Master

Back once again, pre-empting a New Year resolution to re-start my work blog.  Lots of interesting stuff happening, or about to happen, in terms of wind power forecasting for offshore wind farms.

All I will say at the moment is: we're confident we can beat ECMWF at their own game. Watch this space. (nobody else is!).

ERA-Interim

ERA-Interim data is available for free for academic and commercial work:
http://www.ecmwf.int/research/era/do/get/index

Getting ECMWF data

So, you've heard 'ECWMF is far better than  NCEP', and you want a piece of some ECMWF data action. Where do you get get it? How do you find it?.  Here's a list to help you (me) out.

There are multiple routes to ECMWF data heaven. The most flexible and powerful, but most complicated is MARS. To use MARS fully you need the client, a piece of open-source software covered by a GPL licences, but costs £100. It is used to communicate with the data servers at ECMWF, and can subset data on the server side to reduce the volume of data needed to be transferred.  Quicker way to MARS is via  Web MARS.

Some links on MARS
www.ecmwf.int/services/computing/training/material/mars.pdf
www.ecmwf.int/publications/manuals/mars/practice/
www.ecmwf.int/publications/manuals/mars/practice/
www.ecmwf.int/services/archive/

The other route is the ECMWF data server. This gives you straight access, no software needed, no (not many) questions asked. It is actually just a linux server sitting outside the ECMWF firewall running MARS, but only allowing access to the disk archive, and not the tape archive.

To grab a dataset requires a request. Some fields used in the request are:

class	ECMWF classification (od, rd, e4, …)
stream	originating forecasting system(oper, wave, enfo, seas, …)
expver	version of the experiment (01 operational, 11, aaaa)
domain	area covered by the data (Global, Mediterranean, …)
system	seasonal forecast operational system (1, 2, 3)
type	type of field (an, fc, …)
levtype	type of level (pl, ml, sfc, pt, pv)
levelist	levels for the specified levtype (off if levtype=sfc)
param	meteorological parameter (t, temperature, 130, 30.128)
number	ensemble member (1, 2, …)
channel	brightness temperature frequency band
diagnostic	iteration sensitivity forecast products
frequency	direction 2-d wave spectra products
product	section,
latitude	longitude
ocean products

An example of a MARS request is:


retrieve,
class = ei,
stream = oper,
expver = 1,
date = 20071201/to/20071231,
time = 00/06/12/18,
type = an,
levtype = sfc,
param = sd,
target = “era-int.200712.sd”


The request is used to build a MARS tree, which describes the data hierachy to work through.

Truncation is done before interpolation to save resources, and is done depending on the final grid specified. The mapping between the final grid and the resolution truncated to is:

Grid increment Truncation

2.5  ≤ Δ T63

1.5  ≤ Δ < 2.5 T106

0.6  ≤ Δ < 1.5 T213

0.4  ≤ Δ < 0.6 T319

0.3  ≤ Δ < 0.4 T511

0.15 ≤ Δ < 0.3 T799

0.09 ≤ Δ < 0.15 T1279

0.0   ≤ Δ < 0.09 T2047

To subset an area use the boundaries North/West/South/East.

Archived data is available to all registered users. Real-time data only to the big spenders!

ECMWF model data

The global gridded datasets available from central meteorological agencies are amazing. They are the distillation of thousands of global observations from all kinds of sources, pulled together into physically consistent picture of the atmosphere using immensely advanced meteorological models. They are a triumph of the scientific method and the power of cooperation.

They are also confusing as hell to work with. The main problems are:

• multiple similar datasets, often different versions derived from others 
• different vertical coordinate systems
• multiple access routes to the same/similar data 
• no up-to-date description of the current system
• information has to pieced together from multiple, contradictory sources

People assume you're familiar with the models before you start: "Oh, you need the forecast from the T159L60 model". This post is, by definition*, an idiot's guide to ECWMF data.  I'll post this in two parts: the first describing the general system, and the second describing the practicalities actually getting and working with data.

* because it was written by me

ECMWF run the Integrated Forecast System. This is a spectral model: in 2006 the deterministic operational model ran at a resolution of a resolution of T799L91. A summary of this is

Numerical scheme	TL799L91 (triangular truncation, resolving up to wavenumber 799 in spectral space, linear reduced Gaussian grid. 91 levels between the earth's surface and 80 km), Semi-Lagrangian two-time-level semi-implicit formulation.
Time-step	12 minutes
Smallest wavelength resolved	50 km
Number of grid points in model	76,757,590 in upper air and 3,373,960 in surface and sub-surface layers
Grid for computation of physical processes is the Gaussian grid, on which single level parameters are available.	The average distance between grid points is close to 25km.
Variables at each grid point	wind, temperature, humidity, cloud fraction and water/ice content (also (recalculated at each time-step) pressure at surface grid-points), ozone

The spectral resolution refers to the highest wavenumber retained in the solution. A the relation between spectral and spatial resolution is given below.

Spectral	Gaussian	Lat/lon
T63	N48	1.875
TL95	N48	1.875
T106	N80	1.125
TL159	N80	1.125
T213	N160	0.5625
TL255	N128	0.7
TL319	N160	0.5625
TL399	N200	0.450
TL511	N256	0.351
TL639	N320	0.28125
TL799	N400	0.225
TL1023	N512	0.176
TL1279	N640	0.141
TL2047	N1024	0.088

As far as I understand it, the output of the spectral model is converted to reduced gaussian grid, using some kind of magic, very likely involving Fourier transforms.  A gaussian grid plays well with spectral models on a sphere, as the grid points are evenly distributed over the sphere.  This is then converted into a regular latitude-longitude grid (aka cylindrical equidistant aka Plate Carre aka unprojected) as this is what most sane people like to work with.

So remember: the end resolution oft quoted e.g. 1.0 degrees, 0.5 degrees, 0.25 degrees  is the resolution of the end product. It is possible to derive low-resolution datasets from a high resolution model and vice-versa, although obviously going to from a low resolution model to high resolution output will not magically add information.

The resolution above is used for the deterministic forecast.  This also is used to give the initial conditions of the control run of the ensemble forecast. 

The control run (unperturbed) of the ensemble is truncated to a lower resolution from the operational analysis.  In addition, 50 perturbed members make up the full 51-member ensemble.  The ensemble model is run in two legs, with a lower resolution after 10 days.

According to this news article, the resolution of the deterministic and ensemble system were upgraded on 26th January 2010 to the following scheme, which is around 16km resolution!

	Deterministic		Ensemble Prediction System (EPS)
			Leg A		Leg B / C
	Current	Upgrade	Current	Upgrade	Current	Upgrade
Spectral	T799	T1279	T399	T639	T255	T319
Reduced Gaussian grid	N400	N640	N200	N320	N128	N160
Horizontal grid resolution	~25 km	~16 km	~50 km	~30 km	~80 km	~60 km
Dissemination (LL)	0.25	0.125	0.5	0.25	0.5	0.5
Model Level Vertical resolution	91	91	62	62	62	62

In addition, there are multiple reanalysis datasets, most notably ERA-Interim.  This uses cycle 31r2, which is the system introduced in 2006. This is T255 spectral resolution, 60 vertical levels, and a reduced gaussian grid at approximately 79km spacing.

The ERA-interim reanalysis contains analyses at 00z, 06z, 12z and 18z, and 10 day forecasts from 00z and 12z.  Pressure level and surface level data are archived every 3-hours out to 240 hours. Model level data is archived out to 12 hours. Forecast fields are not available on isentropic or potential vorticity levels.

Another product is the multi-model ensemble, where analyses from different centres: ECMWF, UK Met Office, NCEP and DWD are used to drive an ensemble.

Compiling read_wrf_nc with gfortran

Quick but simple tool for extracting time-series from WRF netcdf output. Does not do any de-stagerring of variables.

gfortran lacks iargc, so compile the following

extern int _gfortran_iargc(void);
int iargc_()
{
return _gfortran_iargc();
}

$ gfortran -c gfortran_iargs.c
$ gfortran gfortran_iargc.o read_wrf_nc.f90 -L$NETCDF/lib -lnetcdf -I/$NETCDF/include -ffree-form -o read_wrf_nc 

Installing packages with portable python on windows

Seems to be as easy as:
1. Download and unzip <package>
2. C:\Portable Python 2.7.2.1\App\Scripts\easy_install.exe <package>


Will install into the Portable Python 2.7.2.1\App\lib\site-packages

Job done!

Useful linux commands

Everything except the first line :

sed 1d 

for f in *.txt; do cat $f | sed 1d >> combined.txt; done; 

Gribmaster, cut faster.

Some notes on using the bundle of awesomeness that is the gribmaster. For my own reference, of zero interest to the outside world.

./gribmaster --dset gfs004grib2

Will fetch most recent cycles of GFS model, from servers defined in the config/gfs004grib2.conf file.

The dreaded GRIB format

I'm trying to do some processing of GFS forecast data, to try and find a nice way of interpolating wind speeds to heights above sea level (and eventually above ground level, but that may have to wait).

I thought I could use the Climate Data Operators (CDO), but my grib files appear to have duplicate records in them. Or they do according to wgrib2. To remove duplicate records you can do:

wgrib2 IN.grb -submsg 1 | unique.pl | wgrib2 -i IN.grb -GRIB OUT.grb


Note:
(a) for messages with submessages, only the inventory of the 1st message is checked. i.e. if the 1st submessage matches and the second submessage doesn't, the entire message is removed.
(b) order is preserved
(c) submessage structure is retained, see (a)
(d) note that on the command line, "GRIB" is capitalized

Where unique.pl is:

----------------------- unique.pl ------------------------
#!/usr/bin/perl -w
# print only lines where fields 3..N are different
#
while () {
chomp;
$line = $_;
$_ =~ s/^[0-9.]*:[0-9]*://;
if (! defined $inv{$_}) {
$inv{$_} = 1;
print "$line\n";
}
}
--------------------- end unique.pl ----------------------


UPDATE: This didn't help - it did seem to remove duplicate records, but CDO still didn't like the data. In any case I found out MET doesn't support grib2 format, so I need to convert my GFS files into grib1 in order to do any verification.

I did, however, find an awesome set of tools called gribmaster. Which will fetch grib data from loads of sources. This will be a godsend when I eventually try and make WRF operational, but even now it's very useful as it contains, pre-compiled, wgrib, wgrib2, and gribconv, plus some other tools for working with GEMPAK. In short it is awesome, and knows it.

So, my interim solution is to use 10m wind from GFS for validation against 80m towers (I know), but at least it lets me build the rest of the framework without wasting more time working with grib files. My strategy is this:


grib2 files ----> wgrib2 -----> 10m wind only ------> gribconv ------> grib1 files -----> met


Meteorologists were only ever interested in surface or upper air. The wind industry is interested in 80-100m. Hopefully one day they'll meet.

WRF and wind energy resources

I'll try and maintain a list of WRF and wind energy related resources here. Starting with a good set of presentations from UCAR.
Wind Energy Workshop

Hub-height wind speeds from WRF and GFS

I need to extract wind speed from an atmospheric model at 80m above ground level. That should be pretty simple right? There must be a tool which does that already? Anywhere?

Atmospheric models generally don't like height above ground as vertical coordinate. Most human beings do. This makes getting variables at specified heights a bit of a mission. Especially if the files are in GRIB format, which is a pain in itself.

One option I've have found for WRF is the Universal Post Processor, UPP. UPP will interpolate variables onto a fixed set of heights above sea level, flight levels. These are hard-coded into UPP at: 30m, 50m, 80m, 100m, 305 and a bunch of other heights up to 6000m.

A (not-documented) feature of UPP is that if you specify the vertical levels with a '2', they are interpolated to heights above ground level. Or at least above the level of the model terrain. So put

L = (22220, 0000 ....)

into you wrf_cntrl.parm file, you will get variables out at 30m, 50m, 80m and 100m, above terrain level. Additionally, if you dig into the UPP source code and edit:

src/unipost/MISCLN.f

Then you change the heights to be whatever you want. I settled on every 10m up to 120m. However, a quick check in
src/unipost/FDLVL.f reveals the interpolation is done linearly in height from the levels above and below the height you are interested in. OK if your model levels are close together. However, linearly interpolating wind speed from 150m to 80m is a bit questionable. Don't expect any magic to recreate the vertical wind speed profile.

It would be nice if UPP at least had the option to interpolate in log h. I looked at changing the source code to do this, but then backed away. It didn't look too difficult, I'm just lazy.

Oh - and UPP produces GRIB output. This is fine for me as I need GRIB format for input into the MET tool. But for most uses, GRIB is pretty painful to work with. I mean, fancy taking your lovely NetCDF files and giving you back GRIB. It's just not a deal.

GFS output proved even more of a challenge. UPP can apparently can work directly with GFS data, but there is no documentation on it. I tried it, and it attempted to read GFS files using an i/o module designed to import WRF data, so it fell over and died.

Next I tried using CDO. But I was unable to work with GFS output, as CDO seemed to think there was more than one timestep in my GRIB files (there wasn't).

The only tool left is NCL, which can't write GRIB output. I have implemented some code to which interpolate wind speeds to heights above ground level. My options now are to write that to NetCDF format, and then try and convert to GRIB so that I can input it into MET.

Surely there must be an easier way.

Building CDO with Grib2 support

I was struggling to build CDO with GRIB2 support. Some instructions here:
http://www.eamnet.eu/cms/?q=node/95

Eventually I succeeded using:

$ ./configure --prefix=$HOME/usr --with-netcdf=/cm/shared/apps/netcdf/gcc/64/4.0.1 --with-hdf5=/cm/shared/apps/hdf5/current --with-zlib=/usr --with-szlib=$HOME/usr --with-grib-api=$HOME/usr --disable-shared --with-jasper=$HOME/usr


WRF Reference Configurations

The DTC test centre have a useful set of tried-and-tested WRF configurations here:
http://www.dtcenter.org/config/

These could save me a lot of time and effort...