Tuesday, June 18, 2013

Outline Nested Grid Boundaries in GrADS

If you have ever attended a talk or a presentation, or read a paper for that matter, focused on modeling, you have likely encountered an image showing you the model domain with nested grid boundaries marked inside of it (see below).  This tutorial will show you how to use the 'q w2xy' command to draw these boundaries.  For this tutorial, we will show the NAM domain nested inside of the GFS.

NAM Domain outlined inside the GFS domain.


Lets first start by opening the files, and gathering our relevant information!

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130618/gfs_hd_00z'
    'set mpdset hires'
    set display color white'
    'clear'
    'set clevs 100000'
    'd capesfc'

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/nam/nam20130618/nam_00z'
    'set dfile 2'
    'q file'
    dimslin=sublin(result,5)
    xmax=subwrd(dimslin,3)
    ymax=subwrd(dimslin,6)

 
This code first opens the GFS (parent file) and then displays a dummy variable so that we have a map and a scaled environment. It then opens the NAM (nested file) and gathers the number of x and y grid-points in the domain.  This is important, as we can make use of this data to help set our boundaries.

The simplest way to draw the domain of a nested grid is to take the lower left and the upper right corners of your domain, and draw a rectangle around them.  We will set the nested model to encompass its entire domain and then, by using the 'q dims' and the 'q w2xy' commands, we will convert the lat/lon coordinates to page coordinates and draw our box.

    'set x 1 'xmax
    'set y 1 'ymax  
    'q dims'

    xline=sublin(result,2)
    yline=sublin(result,3)

    lonmin=subwrd(xline,6)
    lonmax=subwrd(xline,8)
    latmin=subwrd(yline,6)
    latmax=subwrd(yline,8)

    lonmin=lonmin+360 

    lonmax=lonmax+360

     *Bottom Left Corner
    'q w2xy 'lonmin' 'latmin
   x1=subwrd(result,3)
   y1=subwrd(result,6)

   *Top Right Corner
    'q w2xy 'lonmax' 'latmax
    x2=subwrd(result,3)
    y2=subwrd(result,6)


What we did here, was just get the minimum lat/lon coordinate, and the maximum lat/lon coordinate and convert them to page coordinates using the 'q w2xy' command.  Note: the longitude values had to be added to 360 to account for the differences between the two domains.  Then we saved each value.  Now, the last thing we need to do is draw the rectangle.  This is done by simply inputting the bottom left and top right corners into the 'draw rec' command.

    'set line 1'
    'set cthick 6'
    'draw rec 'x1' 'y1' 'x2' 'y2


I added a couple of commands to make the line thicker and black, but now you have your nested domain boundaries plotted in your parent grid (see image above).

Download Example Script Here






Friday, June 7, 2013

Tutorial: Tips on handeling more than 1 file at a time in GrADS

While it is possible to use one control file to work with multiple data or binary files by using a few options in the .ctl format (specifically for files of identical spatial dimensions but for different time periods), I find it much easier to move around numerous .ctl files in GrADS.  This tutorial is going to show you how to open more than 1 file at a time as well as address the main hiccups encountered when handling more than one file.

To start out we will need two files to work with.  For simplicity we will use two files with the same spatial dimensions and variables.  Before we start, since we are dealing with more than one file in this tutorial, we will use the 'reinit' command to close all of our current files as to not get ourselves tripped up.  Then, we'll start by opening these two files.

    'reinit'
    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_00z'
    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_12z'

Once our two files are open, we will use the 'q files' (not to be confused with q file) command, to list our open files.
 
Output:

    File 1 : GFS half degree (0.5x0.5) fcst starting from 00Z04june2013, downloaded June 04 04:32 UTC 
    Descriptor: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_h_00z
    Binary: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_00

    File 2 : GFS half degree (0.5x0.5) fcst starting from 12Z04june2013, downloaded June 04 16:33 UTC

    Descriptor: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_h_12z
    Binary: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_12

So basically, you just get a little information regarding your open files, the most relevant of which is your file number.  Since GrADS only works with one file at a time, you need to specify which file you are working with.  This is best done using the 'set dfile' command.  This changes your default file around. So:

    'set dfile 1'  - Points to file 1, in this case 00z GFS
    'set dfile 2'  - Points to file 2, in this case 12z GFS


Once your dfile is set, you treat GrADS as you normally would with just one file.  So, we will continue this example by plotting CAPE from the GFS output, we will start with file 1.

    'set dfile 1'
    'set gxout shaded'
    'set lat 22 50'
    'set lon 230 295'
    'set mpdset hires'
    'set t 5'
    'd capesfc'


Now, lets move on to file 2, and plot CAPE at the same time as a model comparison.  However, since the time of the model initialization has switched, when we work with file 2, we will point to the actual time, not the timestep as we did in file one.  So our command will be 'set time', not 'set t'  First, we will need to get our time though, and we will do that with the 'q dims' command.

By doing that, we get our time is equal to: 12Z13MAY2013.  Now we simply switch files:

    'set dfile 2'
    'set time 12Z04JUN2013'
    'set gxout contour'
    'd capesfc'

CAPE at 12Z June 4th from 00z GFS (shaded) and 12z GFS contoured

Now, we have our comparison map, since these files were identical spatially, we did not need to change the lat/lon dimensions around.  However, I will mention that in general it is not a bad practice to reset all of your space and time dimensions every time you switch from one file to the next, just to ensure that you are plotting everything correctly.

Now, this brings me into the special case:  If both of your files have the same spatial/time dimensions (for example if you were comparing two runs of the WRF, or different GFS ensemble members), you can short cut without having to go through the 'set dfile' commands.

Similar to the GrADS array syntax you can simply put a '.filenumber' at the end of your plotted variable:

   'd var.1'  - Plots var from file 1
   'd var.2'  - Plots var from file 2


Again, only use this shortcut if your model dimensions, are identical, e.g., grid spacing, timestep, etc.  The main problem you will run into handling multiple files in GrADS is an "Entire Grid Undefined" error.  This is simply because you didn't reset the dimensions to fit with your file.  If you get this error, simply do a 'q dims' command, see where x,y,z,t fall outside of your file dimensions.

So that is all for the basics of file handling in GrADS.  Hope you found it useful!

Download Example Script