[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: fortran code for station data
Dear Tim,
Below is a bunch of various emails that deal with setting up
GrADS for station contour plots, etc... (along with fortran code)
Hopefully, this is enough to get you started... good luck!
Cheers,
Adam
P.S. there is a nifty trick given at the end of these emails for doing
these plots!
============================================================================
Hi Greg,
Well, you are certainly right - GrADS cannot plot station data. What
you need to do is to interpolate the station data onto a grid and then display
that gridded field...
Fortunately, GrADS does all this for you. If you look on Page 83 of
the V1.5.1.12 postscript documentation that has just recently been released,
there is a function "oacres" (Cressman objective analysis) that is described.
This is the function you want. Basically, you just need to set up a gridded
template (has to be linear in lat-lon) which GrADS will write the interpolated
field onto. You can then do contour,shaded,etc. plots with this field.
The documentation seems to do a fairly good job of describing the
syntax. However, should you run into problems don't hesitate ask for
further help...
Good luck,
Adam
--
============================================================================
* * * * *
C. Adam Schlosser o () o * * *
GFDL/NOAA * / \/||\/ \ Think *
Forrestal Campus, US Rt. 1 / || \ * snow....
Princeton, NJ 08542 * / \\\\ * \ * *
x \// x
Phone: (609) 987-5054 * * \\ * *
Fax: (609) 987-5063 \\ * *
email: cas@gfdl.gov * * ^^
============================================================================
> Adam,
> Thanks for the reply. What is the format of the gridded
> template that GrADS needs to interpolate the field onto? Do you have
> an example?
>
> --
> Thanks for your time,
>
> Greg
The format of the template is just like any other GrADS gridded
data set you would construct. Usually, I just make a template grid with
each value at each grid point having a value of zero (so, you just have
a constant value grid field of zeros). Just remember that the gridded
field has to be LINEAR in lat-lon. Also, you only have to make your
gridded template cover the region that your station data covers. In
other words, you don't need to make a GLOBAL grid template for station
data that only covers say the U.S.
I guess one more crucial element to the template set-up is your
resolution. You'd like to make your resolution enough so that
they're aren't alot of station points within one gridpoint domain
especially if the station data shows alot of spatial variability (i.e.
if you have alot of station points within one gridpoint in your template,
oacres will not pick up the spatial variability of the station data)...
So just set-up your resolution of your template so that it avoids having
multiple station pts. within a gridpoint domain.... BUT... BUT.... don't
make your template resolution alot higher than your station density -
the oacres function will then tend to put unrealistic extremma values
for a grid point that may lie between two station points.
As for an example, well - say I needed to perform oacres for
station data covering a region that roughy spans a 20x20 degree domain -
and I need a 1x1 degree resolution. In FORTRAN, I would set-up the
template as follows
real template(20,20)
.... open a unit 15 for GrADS formatting, then
.... just set the grid values to zero
do 10 i=1,20
do 10 j=1,20
template(i,j)=0.0
10 continue
write (15,rec=irec) template
To save diskspace, you can really just write the template for
one timestep. Oh yeah, you want to make sure that your time domains of
your template and station data match up (if they don't, I'm pretty sure
you get an error). So just make sure that the first (and only) timestep
of your template matches the first timestep of your station data (this
would be your TDEF line in the descriptor files). Also, set-up your
XDEF and YDEF lines so that they cover the region you're looking at.
So then, say I have surface temp for the station data (GrADS variable
"tair") - and my variable name for the GrADS template variable is "templ".
To peform oacres I just type:
d oacres(templ.2(t=1),tair) (Note: I have opened two descriptor files
the first being the station data and
the second being the template)
I included the "(t=1)" in the templ.2 variable to note that
if your station data has many timesteps - you can do oacres for any
timestep of the station data so long as you specify that your grid
template is at t=1. (remember? I suggested that you really only need to
write your template for one timestep...)
I would strongly suggest then displaying your station data to
compare with the gridded field you have just plotted. Note again that
you can do contours, shaded, and other GrADS type functions with the
oacres created field. However, in some cases (large station data or high
station density) oacres can take a while. What I would suggest is that
once you have checked out the oacres field against the station data and
the analysis looks reasonable, use the "set gxout fwrite" option and
write the grid interpolation into a GrADS file. After that, you can
just read the interpolated grid field that you've just written to
"grads.fwrite" (or whatever file name you choose!).
Well, this should be enough to get you started. Good Luck!
------------------ HERE'S THE TRICK ----------------------
Instead of going through the trouble of actually creating the
template file try this:
------------- cut here for sample template.ctl file ------------
dset /tmp/%y2
undef 1e20
title dummy template grid data for oacres
options template
xdef 360 linear 0 1
ydef 181 linear -90 1
zdef 1 levels 1013
tdef 1 linear 1jan1901 1mo
vars 1
dum 0 0 dummy grid
endvars
----------- end of template.ctl ---------------
dset points to a non existant data file so when you,
d dum
in GrADS, a grid is CREATED with undefined data points. Thus,
d oacres(templ.2(t=1),tair) (Note: I have opened two descriptor files
the first being the station data and
the second being the template)
can become,
d oacres(dum.2(t=1),tair)
AND, most importantly, the resolution, etc. of the template grid
can be changed by simply changing the xdef and ydef cards in
template.ctl. Thus, YOU HAVE AVOIDED CREATING FICTIOUS DATA!!!
Good luck and CC4N,
Mike
--
============================================================================
* * * * *
C. Adam Schlosser o () o * * *
GFDL/NOAA * / \/||\/ \ Think *
Forrestal Campus, US Rt. 1 / || \ * snow....
Princeton, NJ 08542 * / \\\\ * \ * *
x \// x
Phone: (609) 987-5054 * * \\ * *
Fax: (609) 987-5063 \\ * *
email: cas@gfdl.gov * * ^^
============================================================================