Outputting NetCDF?
Question:
How do I make my model write out NetCDF files that are suitable for Ferret?
Solution:
You can get the NetCDF libraries from Unidata at http://www.unidata.ucar.edu/software/netcdf>, link your software to the libraries, and make calls to the NetCDF library to open a dataset, define variables, and write data to the file. Below is a second way, using the NetCDF tools ncdump and ncgen:
The NetCDF utilities ncgen and ncdump do most of the work for you. They actually generate the FORTRAN or C code that you will need. All you have to do is to "hack it" a little bit.
You can use Ferret commands to describe what your output files should look like -- everything about the output files except the actual data and time word values.
We proceed in 4 steps:
- describe the file using Ferret commands
- use ncdump to create a "CDL" language description of the file
- use ncgen to generate FORTRAN or C code
- hack that code and insert it into our model code
1. Describe the file using Ferret commands
For this example let us suppose that we will have variables U and TEMP and that they are on staggered grids.
We will define variables like this with Ferret commands. The **values** of these variables will never be used ... but the **shape** of the variables is important. Our example doesn't use a Z axis but it could be added in a straightforward way.
The encoding of time is discussed at some length in the Ferret Users Guide
! define the axes and grids ! temperature axes define axis/x=130e:80w:2/unit="degrees" xt define axis/y=20s:20n:2/unit="degrees" yt ! velocity axes (staggered relative to temperature) define axis/x=131e:79w:2/unit="degrees" xu define axis/y=19s:21n:2/unit="degrees" yu ! time axis define axis/t="1-jan-1980:12:00":"1-jan-1990:12:00":1/unit="days" time ! define the grids define grid/x=xt/y=yt/t=time gt define grid/x=xu/y=yu/t=time gu ! define the (dummy) variables let u = SIN((x[g=gu]+y[g=gu]+t[g=gu])/100) let temp = SIN((x[g=gt]+y[g=gt]+t[g=gt])/100) set variable/title="Temperature"/units="centigrade" temp set variable/title="Zonal Velocity"/units="m/sec" u ! create the template NetCDF file save/file=template.cdf/l=1 temp,u
2. Now create a CDL description of this file
! ... note -- you can edit template.cdl with your favorite editor > ncdump template.cdf > template.cdl
3. now create FORTRAN code that would generate this NetCDF file
( use the "-c" qualifier instead of "-f" to get C code)
! ... note: you can "hack" this FORTRAN code to create the NetCDF output ! parts of your model > ncgen -f template.cdl > template.F
4. Hack the code from template.F into your model code
Usually you will
- create a subroutine to create (but not insert data into) the output file. It is executed only once.
- insert the NetCDF into your model code where you want to write out the time word and the values of your variables. It is executed once per output time step of your model
- The subroutine to create the output file is the FORTRAN code from the line
include 'netcdf.inc'
through
* leave define mode call ncendf(ncid, iret)
You will also want the blocks of code like store XT that write your models coordinate system into the output file.
Note: For maximum accuracy you may want to write the coordinates directly from your model arrays instead of from the DATA statement. For this make sure youinsert the matching data type (e.g. "NCDOUBLE") in the define variables section of the file creation code.
- The code to write the time word is hacked from store TIME. Note that the"corner" value will be "1" for time step 1, 2 for step 2, etc.As with the other coordinates make the data type of your model time word variable matches the data type specified during define variables
The code will look something like
corner(1) = istep ! the current output time step edges(1) = 1 call ncvpt(ncid, TIMEid, corner, edges, time_word, iret)
The code to write the model variables is similar:
! corner is the start point for writing -- normally "1" on each axis ! edges is the amount to write for each axis -- normally the full axis length corner(1) = 1 corner(2) = 1 corner(3) = istep edges(1) = 76 edges(2) = 21 edges(3) = 1 call ncvpt(ncid, TEMPid, corner, edges, temp_var, iret)
Last modified: Dec 9, 1996