Monthly Archives: January 2012

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

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.

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