Data Resources: PSD netCDF Support Routines: Reading netCDF data in Fortran Code


We are strongly discouraging use of "call netopen" and instead you should try to use the standard nf_open and nf_nvid calls (ersion 3 of the calls).

This document for version 1.3 of the netopen package supersedes all previous versions (see Cathy Smith in rm. 1D-106, phone: 497-6263, or e-mail cathy.smith@noaa.gov for any questions).


To link to the gfortran library on linux quartet:

gfortran -finit-local-zero -I/usr/local/gfortran/include/ -o file file.f /home/dmwork/DataWorkCommonCode/netopen/read/src/out/lnx64/gfortran/libnetopen.tmp.a -lnetcdff -ludunits
To link to the sun library on linux quartet:
f90 -o cons4x10.nocrd cons4x10.nocrd.f /home/dmwork/DataWorkCommonCode/netopen/read/src/out/lnx64/sunf90/libnetopen.a -lnetcdff -ludunits

To link to the gfortran library on linux quartet:
gfortran -finit-local-zero -I/usr/local/gfortran/include/ -o file file.f /home/dmwork/DataWorkCommonCode/netopen/read/src/out/lnx64/gfortran/libnetopen.tmp.a -lnetcdff -ludunits
For a MAC, see Dave Allured.

Jump to Changes, Using the New Version, Missing Values, Examples, Troubleshooting, Datasets, Variables and Levels, or


Latest changes 05/09/2003 by CAS

The changes to netopen should allow all COARDS-CF compatible files to be read via gridread and similar routines. In addition, you should be able to read all netCDF files that are COARDS compatible but don't have udunits times. In those cases, gridread and similar routines will return -99,-99,-99,-99 for year, month,day and hour. Routines like timeindex or dayread that rely on time will NOT work for those files.


Using the Version 1.3 Netopen routines

For those already familiar with netopen skip to linking commands, netop_init, gridread(x), dayread(x), specread, timeindex,headread, specread subroutine has been added to allow users to read spectral coefficient files.

The gridread, gridreadx, dayread, and dayreadx subroutines are set up to read rectangular or Gaussian gridded data of up to 1441x721 points and return an array starting at the South pole and proceding North regardless of how the data is stored in the actual netCDF file. The specread subroutine is capable of reading up to 4032 spectral coefficients (truncation of 124). Further, since dynamic allocation is not standard in Fortran 77, it is important that you pass in an array variable which is dimensioned according to the size of the grid (or vector of spectral coefficients) in the netCDF file. (Note: you can determine this by using the ncdump -c command on the file). If you pass an array that is too small the software will crash, and if it is too large the indexing will be wrong.

This software can be used to read grid sizes larger than 1441x720 by changing the maxlat and maxlon parameters to the appropriate grid size. Additionally, by changing the maxnumvalues parameter the spectral coefficient truncation size readable can be increased.

Basic steps to reading a netCDF file in fortran code:

  1. Put the following at the top of your program (starting in column 7 or greater):
          include '/usr/local/include/netcdf.inc'
    
  2. Call the netopen initialization routine netop_init with an empty argument list.
  3. Open one or more files. Also, retrieve identifying variable IDs for each variable in each file.
  4. If looping, find start and stop time index.
  5. Read in a grid (looping through data, if desired).

The specific routines that accomplish the above are as follows:

  1. Call the initialization routine.

    call netop_init() - Note: you need only call this routine once in your program

  2. Locating and opening files and returning variable IDs for each variable.

    To obtain file and variable id for each file, use the following:

    inet=ncvid("filepath/name",0,ierr)
    ivar=ncvid(inet,"variable",ierr2)

    The filepathname should be the full name and might be something like /Datasets/cmap/enh/precip.mon.mean.nc". ierr and ierr2 are integers which can hold the return error code. variable is the netCDF variable such as "sst" or "hgt". You will use inet for the file id and ivar for the variable id for a particaular file.

  3. Read the data.

    To read a (previously opened) file and retrieve a rectangular or Gaussian grid (depending on the file that was opened) for one time period, you can use either gridread or gridreadx. The grid returned will be ordered South to North regardless of the order of data in the netCDF file. These routines have identical interfaces except that gridreadx returns an extra longitude vector in the grid (a copy of the longitude vector in column 1 is placed in the last column e.g. 145x73 for 2.5 degree data). Here is the details of calling gridread:

    call gridread(iadd,itime,x,flevel,netid,ivarid,ny,nm,nd,nh)

    iadd   = integer- indicates whether subr. will automatically advance 
    	 time step
    itime  = input time index for that file -integer
             itime is 1-based.  The first time step in the file is itime = 1.
    x      = returned real*4 array, S-N dimensioned according to netCDF 
    	 file info. You MUST dimension this array correctly!
    flevel = input floating point level (e.g. 850.0) or 0.0 if no
    	 level.
    netid  = input integer file ID (gotten from netopen)
    ivarid = input integer variable ID (gotten from netopen)
    ny     = returned integer year (2 digits if 1901 - 1999, 4 otherwise)
    nm     = returned integer month
    nd     = returned integer day
    nh     = returned integer hour
    
    If you wish to have the subroutine increment the time index for you automatically, set iadd=1 for gridread. The subroutine will automatically advance the time increment for you at the START of the subroutine. So, be sure to set itime to 1 less than the first time index if you use that option.

    To read more than one netCDF data variable or level at a time using the same time index program variable, you should set iadd=0 for the calls to retrieve the additional variables or levels. Time index values are specific to a dataset. So if you are reading one dataset that starts on Jan 1, 1985 and another dataset that starts on Dec 1 1984, the time index value for the same time will NOT be the same. Therefore it is advisable to use a different program variable for each of the datasets.

    gridread and gridreadx will advance to the next file if the read extends beyond the last record of the current file. (Like netread2 in the previous version). To read more than one data variable at the same time, set iadd=0 for the other data variables. Alternatively, use a separate time index variable for each data variable.

    All functions in this library return grids from South to North regardless of how the data is stored in the netCDF file.

  4. Read spectral coefficient data.

    To read a (previously opened) file and retrieve a vector of spectral coefficients for one time period, you can use either specread. This routine has an identical interface to that of the gridread routine. Here is a description of the use of specread:

    call specread(iadd,itime,x,flevel,netid,ivarid,ny,nm,nd,nh)

    iadd   = integer- indicates whether subroutine will automatically advance 
    	 time step
    itime  = input time index for that file -integer
    x      = returned complex array dimensioned according to the number of
             spectral coefficients in the file. You MUST dimension this 
             array correctly!
    flevel = input floating point level (e.g. 0.995 sigma) or 0.0 if no
    	 level.  Surface is 2000.0 
    netid  = input integer file ID (gotten from netopen)
    ivarid = input integer variable ID (gotten from netopen)
    ny     = returned integer year (2 digits if 1901 - 1999, 4 otherwise)
    nm     = returned integer month
    nd     = returned integer day
    nh     = returned integer hour
    
    If you wish to have the subroutine increment the time index for you automatically, set iadd=1 for gridread. The subroutine will automatically advance the time increment for you at the START of the subroutine. So, be sure to set itime to 1 less than the first time index if you use that option.

    To read more than one netCDF data variable or level at a time using the same time index program variable, you should set iadd=0 for the calls to retrieve the additional variables or levels. Time index values are specific to a dataset. So if you are reading one dataset that starts on Jan 1, 1985 and another dataset that starts on Dec 1 1984, the time index value for the same time will NOT be the same. Therefore it is advisable to use a different program variable for each of the datasets.

    gridread and gridreadx will advance to the next file if the read extends beyond the last record of the current file.


    To find a time index associated with a certain file and date, you can call:

    call timeindex(netid,ny,nm,nd,nh,index)

    netid = file ID- integer
    ny    = year (4 digit default or 2 digit for 1901-1999) -integer
    nm    = month (1-12) -integer
    nd    = day -integer
    nh    = hour (0-24) -integer
    index = returned time index -integer
    

    To read a grid for a particular time, one can use:

    call dayread(iadd,itime,x,flevel,netid,ivarid,ny,nm,nd,nh)

    Where ny,nm,nd and nh are specified (iadd and itime are ignored on input, itime is returned). Or if you would like an extra longitude vector (like gridreadx) call dayreadx which has the same parameter list as dayread. The grid returned will be ordered South to North regardless of the order of the data in the netCDF file.


    Some of the data (the NCEP once daily, EC and NCEP twice daily) contain a variable that indicates whether the field was interpolated. To read this field, use:

    call nf_vid(netid,"head",headvarid)

    call headread(iadd,itime,nhead,flevel,netid,headvarid,ny,nm,nd,nh)

    which has the same parameters as netread2 with the exception that nhead is a single integer. nhead will take on one of four values indicating the following interpolation information for the grid represented:

    	0) no interpolation
    	1) both hemispheres interpolated
    	2) just NH
    	3) just SH
    


    Then use the automatic increment feature of gridread(x) and when the record beyond the last record in dataset 1 is read, the same variable and statistic (if existent) will be opened in dataset 2 and the netid and variable ids will be set properly.

    Currently the following datasets are linkable in this manner:


Missing Values:

Whenever any of the netopen routines encounters a value in the netCDF file which has been defined as missing, it will replace that value with the value 1e30.

Examples:

Let's say you wanted to read 3 years (1985 - 1987) of NCEP u and v winds at 200mb. Your program could look something like this:

      real u(144,73)
      real v(144,73)

      call netop_init()

      do iyear=1985,1987
          write(cy,222)iy
222      format(i4.4)
          file1=,'/Datasets/ncep.reanalysis.dailyavgs/pressure/uwnd.'//cyear//.nc'
          file2=,'/Datasets/ncep.reanalysis.dailyavgs/pressure/vwnd.'//cyear//.nc'
      call nf_open(netidu,file1,ierr)
      call  nf_open(netidv,file2,ierr)
      call nf_vid(netidu,'uwnd',ierr1)
      call nf_vid(netidv,'vwnd',ierr1)

         do  iday=1,365
            call gridread(0,iday,u,200.0,netidu,ivaridu,ny,nm,nd,nh)
            call gridread(0,iday,v,200.0,netidv,ivaridv,ny,nm,nd,nh)

C           User functionality with this timestep u and v winds at 200mb.

         enddo
      enddo
If you wished your program to resemble the sequential type read that you may have used in previous programs, (where you needed the extra longitude vector for IDL or something like that) your code could look something like this:

      real u(144,73)
      real v(144,73)

      call netop_init()

          file1=,'/Datasets/ncep.reanalysis.dailyavgs/pressure/uwnd.'//cyear//.nc'
          file2=,'/Datasets/ncep.reanalysis.dailyavgs/pressure/vwnd.'//cyear//.nc'
      call nf_open(netidu,file1,ierr)
      call nf_open(netidv,file2,ierr)
      call nf_vid(netidu,'uwnd',ierr1)
      call nf_vid(netidv,'vwnd',ierr1)

100   call gridreadx(1,iday,u,200.0,netidu,ivaridu,ny,nm,nd,nh)
      call gridreadx(0,iday,v,200.0,netidv,ivaridv,ny,nm,nd,nh)

C	 ...
C        ...process u,v any way you wish...
C	 ...

      if((ny.eq.87).and.(nm.eq.12).and.(nd.eq.31))iend=1
      if(iend.ne.1)goto 100
To open the coads sst, monthly mean data and read Jan 1975-Feb 1976:
      real sst(180,90)     ! Array will be returned S-N

      call netop_init()

      call timeindex(netid,1975,1,0,0,itime1)
      call timeindex(netid,1976,2,0,0,itime2)
      do it=itime1,itime2
          call gridread(0,it,sst,0.0,netid,ivarid,ny,nm,nd,nh)
      enddo
To read NCEP/NCAR Reanalysis heights at 50mb through all of 1986:
      real x(144,73)
      call netop_init()


      it1 = 0

      do it = 1,1500

         call gridread(1,it1,x,50.0,net2,ivar2,ny1,nm1,nd1,nh1)

         if (ny1.ne.1986) then

C            User functionality with this timestep 50mb height

         endif
      enddo
Reading NCEP/NCAR Reanalysis surface Gaussian air temperatures (note: this set of calls will read all of the 86 data - 1460 timesteps, and into the 1987 data - remaining 40 timesteps).

      real x(192,94)

      call netop_init()

      call netopen(1986,'nmcrean','gaus','air','',2000.0,net2,ivar2)
CC here 
      it1 = 0

      do it = 1,1500    

         call gridread(1,it1,x,2000.0,net2,ivar2,ny1,nm1,nd1,nh1)

C        User functionality with this timestep surface air temperature
C        Gaussian grid.

      enddo
And finally, for reading a spectral coefficient file from NCEP/NCAR Reanalysis looking at the specific humidity variable at 0.995 sigma.
      complex x(2016)

      call netop_init()

CC      call netopen(1986,'nmcrean','spec','shum','',0.0,net2,ivar2)

      it1 = 0

      do it = 1,1460

         specread(1,it1,x,0.995,net2,ivar2,ny1,nm1,nd1,nh1)

C        User functionality with this timestep specific humidity spectral 
C        coefficients at 0.995 sigma

      enddo
A sample working program is in example.f

Possible source of problems:

If you run into any problems when you run your netCDF reading program, see if any of the below might apply. If not, I will be happy to look at your code. If you discover any bugs or non-robustness of any of the subroutines or if you have any suggestions, please bring them to my attention A.S.A.P..

  1. If you run into problems because your code has too many files open, you can use the NCCLOS command to close any files that you are finished with. The syntax is CALL NCCLOS(nc_fileid, RCODE).
  2. Make sure that the size of the grid or vector (spectral coefficients) is the same as that in the netCDF file (use ncdump -h to learn about the dimensions)
  3. Be sure that all subroutine parameters are the right type (i.e. integer,real, character, etc). NetCDF is very picky about such things. Especially the new floating point level value!
  4. Don't ask for data that does not exist. For example, there are no relative humidities above 300mb and you can't get slp at 1000mb for the NCEP data. Also, make sure you have spelled all dataset and variable names correctly (see below).
  5. Make sure that you are using the appropriate variable IDs and file IDs.
  6. If you choose to have the subroutine keep track of the time index, don't have it advanced for other data variables you wish to read at the same time (unless you use different names for the time index variable).

Datasets Available

The data descriptions have more detailed descriptions of what is in a particular dataset (e.g. dates, variable names, levels, statistics), references for that dataset and caveats on the use of the data. We also provide a summary of PSD datasets. ASCII versions of the dataset descriptions are in the README files located in the ~ftp/Datasets subdirectories.

Variable and Level Names

Levels in netopen are needed when files in the datasets are divided by levels. This applies to the nmc2d, nmcltm, and cddb datasets. Levels are needed in the gridread(x) and dayread(x) commands when a particular file is divided into different levels. This applies to the nmc, nmcrean, nmcp, and nmcmm. Levels should be specified as floating point values (1000.0, 850.0, 200.0, etc.).

There are several "special" levels in the netopen call. The NCEP/NCAR and ECMWF surface is represented as 2000.0 The NCEP/NCAR Reanalysis tropopause level is represented as 3000.0. Use 0.0 as the level for the NCEP/NCAR Reanalysis pressure level data and the NCEP (Leetmaa) Ocean data and other data sets with multi-level files.

When referring to variable names, the following are applicable in netopen - Note that not all datasets have all variables (in fact as you'll note below some datasets use different variable names for the same physical variable) so you should check the PSD Data Directory Dataset List for the variable(s)/dataset(s) in which you are interested. These are sorted by netopen variable id for those variables that are in more than one dataset.