PSD netCDF Support Routines: Writing netCDF data in Fortran Code


This document for version 1.0 of the netwrite package supersedes all previous versions. See Cathy Smith e-mail Cathy Smith@noaa.gov) for assistance.


Using the Version 1.0 netwrite routines

The netwrite package includes functions to create a netCDF file which conforms to the COARDS netCDF standard , as well as functions for writing data to that file. The files are assumed to be grids (lat by lon) that can be either regularly or semi-regularly spaced. They can have one or more levels. The data should be real, although the routines could be adapted for data of type double or integer.

Output files are in the COARDS standard and will be readable by netopen, netread, NCL, GrADS, xgplot, FERRET, Freud, etc. Some netCDF attributes may have default values-- this will work, but you may wish to change them. A subroutine is available to make changes to the global attributes. Another will soon be available to change other less-used attributes.

It is not necessary to both create a netCDF file and store the data within the same code (as is done here). If it is more convenient, the netCDF file can be created in a separate program and then written to later. Likewise, if you are familiar with cdl's and generating netCDF files that way, you can just use the store routine (used in this code) to write the data to an "already created" netCDF file.

The subroutines are located in the file ~crdc/netopen/write/netwrite.f. Alternatively, if you do not plan on making any changes to the subroutines, you can link with the library netwrite by including:

   -lnetwrite -lnetcdf -ludunits -luduerr
at the end of your compile and load command.

Basic steps to writing a netCDF file from FORTRAN code:

  1. Create real arrays large enough to store all needed values for each of longitude, latitude, and levels.

  2. Create an integer(4) array to represent the timestep of the data. The first array value represents the yearly time step, the second represents the monthly time step, the third value represents the daily time step, and the last value represents the hourly time step. For daily data the array would have values (0,0,1,0). Twice daily data would have values (0,0,0,12).

  3. Declare a real variable called xmiss and assign it a value to be interpreted as missing. Its value must not be the same as any of your legitimate data values. (If you don't set a missing value, the value will be 0.0!) If you have no missing data, use something like -9999 or 1.e30.

  4. Create a netCDF file, and save the returned handles for the file and variable id.

    IMPORTANT:If you call netcreate with an existing file name, all of the data in that file will be destroyed.

    call netcreate(prefix, short_varname, long_varname, dataset_name, varunits, delta_time, nlon, xlon, nlat, xlat, nlev, xlev, levunits, xmiss, fileid, varid)

    TYPE      NAME                         DESCRIPTION
    =====  ============    ==============================================
    char    prefix          output filename prefix (.nc will be appended)
    char    short_varname   short variable name (eg, sst or uwnd)
    char    long_varname    descriptive name (eg, zonal wind component)
    char    dataset_name    name of data set (eg, NMor Cathy's Data)
    char    varunits        variable units (eg, m/s)
    int     delta_time      integer array of time step (see step 2).
    int     nlon            number of longitudes
    real    xlon            array of longitude values
    int     nlat            number of latitudes
    real    xlat            array of latitude values
    int     nlev            number of levels (for no levels use 1)
    real    xlev            array of levels  (for no levels use something
                                             like 0.0, or 2000. for sfc.)
    char    levunits        level units (eg, millibar or sigma_level)
    real    xmiss           missing value (see step 3)
    int     fileid          returned file id for future netCDF calls
    int     varid           returned variable id for future netCDF calls
    
    To create a packed data file, the call is similar but a scale and offset must be included.

    call netcreatpack(prefix, short_varname, long_varname, dataset_name, varunits, delta_time, nlon, xlon, nlat, xlat, nlev, xlev, levunits, xmiss, xscale,xoffset,fileid, varid)

    The scale and offset are used to pack a real*4 array into an integer*2 array. You can look at other files to get a scale and offset (e.g. the NMC reanalysis) or you can calculate them with the following subroutine:

    call offsetscale(xlow,xhigh,iprec,xscale,xadd)

    where xlow and xhigh are the valid ranges of the data
    iprec is the required precision (number of decimal places)
       if iprec=-99, then the highest precision possible is used
       xscale and xadd are the returned scale and offset
    
  5. Optional: Set the values of the global attributes.

    call createopt1(fileid, history, description, title, platform)

    TYPE      NAME                         DESCRIPTION
    =====  ============    ==============================================
    int     fileid          returned file id from netcreate subroutine call
    char    history         may contain date of creation and person 
                           creating file
    char    description     long description of the dataset 
                           (caveats, supplier, etc.)
    char    title           title of the datafile (eg, Reynold's Blended SST)
    char    platform        model, analysis, observations, etc.
     
    
  6. Optional: Add additional variables to file

    To add additional unpacked variables to a file, use the following:

    call subroutine varadd(varname,vardes, & vunit,dataset,xmiss,fileid,ivarid2)

    TYPE      NAME                         DESCRIPTION
    =====  ============    ==============================================
    char    varname        Second variable shortname
    char    vardes         Long variable description
    char    vunit          Variable units
    char    dataset        Name of dataset
    real    xmiss          Missing value indicator
    int     fileid         File id
    int     ivarid2        Variable id (different than first variable id!)
    
    To store additional packed variables, call:

    call subroutine varaddpack(varname,vardes, & vunit,dataset,xmiss,fileid,ivarid2,xscale,xoffset)

  7. Loop through times and levels for each grid and store data in the netCDF file.

    call storeunpack(data, iyr, imon, iday, ihr, itime_step, nlon, nlat, nlevel, xmiss, fileid, varid)

    TYPE      NAME                         DESCRIPTION
    =====  ============    ==============================================
    real    data            2-d array of data being stored
    int     iyr             the year   (1900 will be added, or you may specify
                                       a four digit year.  If data represent
                                       long term climatological values, set iyr = -1
                                       and the metadata will contain year=1)
                           
    int     imon            the month  (1-12)
    int     iday            the day    (1-31)
    int     ihr             the hour   (0-24)
    int     itime_step      time index (Increase by 1 as you loop through
                                       time steps.  If you don't, the  
                                       data will be overwritten!)
    int     nlon            number of longitudes
    int     nlat            number of latitudes
    int     nlevel          number of level being stored
    real    xmiss           missing value in the array (Don't set to value in
                                       the array unless it is missing!)
    int     fileid          returned from netcreate, ncopn, or netopen
    int     varid           variable id (from netcreate, ncopn, or netopen)
    
    To store packed data, a similar call is used but with a integer*2 work array of the same dimension as data in the call.

    call storepack(data, work,iyr, imon, iday, ihr, itime_step, nlon, nlat, nlevel, xmiss, fileid, varid)

    Data are packed to short integer as follows:

    		short_value = (float_value - add_offset)/scale_factor
    
    When the data are read by the application program, the floating point value is recovered from the packed data as follows:
    		float_value = scale_factor*short_value + add_offset
    

  8. An optional subroutine is available to calculate the next year, month, day and hour given the current values and the timestep. The values for year, month, day, hour and delta_time should be initialized before the first call to this routine, and the routine should be called after each time step is written. The call is:

    call nexttime(iyr, imon, iday, ihr, delta_time)

    TYPE      NAME                         DESCRIPTION
    =====  ============    ==============================================
    int     iyr             the year   (current year value; see notes above)
    int     imon            the month  (current month value in the range 1-12)
    int     iday            the day    (current day value in the range 1-31)
    int     ihr             the hour   (current hour value in the range 0-24)
    int(4)  delta_time	the array specifying the time step
    

    IMPORTANT: close the file after you have written it, or the results will not be saved.

    Call ncclos(fileid, icode)

    TYPE      NAME                         DESCRIPTION
    =====  ============    ==============================================
    int     fileid          file id used in calls to write the file
    int	icode		error code (returned)
    

  9. The netwrite subroutines are in ~crdc/netopen/write/netwrite.f. You may either link with libnetwrite.a, or compile your code with the netwrite source:

    f77 -o netwrite netwrite.f netwrite.o -lnetcdf -ludunits -luduerr

  10. Run cdcncdump prefix.nc to check the output file.

Example:

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  This program is a "driver" for Cathy Smith's FORTRAN code for
C  generating COARDS compliant netCDF files.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  
      implicit none      
      real  data(144,73)     ! Part of the example program

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  REQUIRED ITEMS  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  DECLARATIONS
      integer      nlon, nlat, nlev 
      parameter    (nlon=144, nlat=73, nlev=1)   
      integer      i, j, nlevel, icode, delta_time(4), fileid, varid
      integer      iyr, imon, iday, ihr, itime_step
      character*80 prefix, long_varname, dataset_name
      character*30 short_varname, varunits, levunits
      real         xmiss, xlat(nlat), xlon(nlon), xlev(nlev)

C  OPTIONAL CALL #1 *** for specifying global attributes ***
      character*132 history, title, description    ! Global attributes
      character*20  platform                       

C  ASSIGN VALUES.      
      prefix        = 'testcdfdata'                ! Output file - .nc
      short_varname = 'vwnd'                       ! Short variable name
      long_varname  = 'meridional wind component'  ! Variable description
      dataset_name  = 'NMC daily analysis.'        ! Data source
      varunits      = 'm/s'                        ! Variable units
      xmiss         = -9999.                       ! Missing value
      delta_time(1) = 0                            ! Delta time Years
      delta_time(2) = 0                            ! Delta time Months
      delta_time(3) = 0                            ! Delta time Days
      delta_time(4) = 12                           ! Delta time Hours
      
      do i = 1, nlon
         xlon(i) = 0.  + float(i-1)*2.5            ! Longitude(s)
      end do

      do i = 1, nlat
         xlat(i) = 90. - float(i-1)*2.5            ! Latitude(s)
      end do

      xlev(1)       = 300.                         ! Vertical level(s) 
      levunits      = 'millibar'                   ! Level units
      
C  OPTIONAL CALL #1 *** set values of the global attributes ***
      history       = 'Created by Cathy Smith 12/12/95 using netwrite.f'
      description   = 'Data were artifically generated and represent '
     &              //'nothing in particular.'
      title         = 'Cathys fake data.'
      platform      = 'Observational data.'

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  SAMPLE PROGRAM  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
C  Create netCDF file template.
      call netcreate(prefix, short_varname, long_varname, dataset_name, 
     &               varunits, delta_time, nlon, xlon, nlat, xlat, nlev,
     &               xlev, levunits, xmiss, fileid, varid)

C  OPTIONAL CALL #1 for setting global attributes.
      call createopt1(fileid, history, description, title, platform)

C  Create fake data.  These represent a 2-d slab at some level.
      do i = 1,144
         do j = 1,73
            data(i,j) = float(i+j)
         enddo 
      enddo

C  **** Write data to netCDF file ****
C  Loop through levels then time.  
C  Update iyr, imon, iday, ihr, and increment itime_step by one (1).
      itime_step =  1    ! Initial value.
      iyr        = 95    ! Set initial time parameters before 
      imon       = 12    ! writing data to the file.
      iday       = 25
      ihr        =  0

C  Begin looping.
      do i = 1, 32                   !(Begin time loop, 32 time steps)

         do nlevel = 1, nlev         !(Begin level loop)
C  HERE ==> One would normally read 2-D grid (data) at level=nlevel
            call storeunpack( data, iyr, imon, iday, ihr, 
     &           itime_step, nlon, nlat, nlevel, xmiss, fileid, varid)
         end do  !*** End level loop

         call nexttime(iyr, imon, iday, ihr, delta_time)
         itime_step = itime_step + 1
      end do     !*** End time loop


C  Close ouput data file.
      call ncclos(fileid, icode)     ! Output file

      stop
      end
Example program courtesy of Andrew F. Loughe (andrew.f.loughe@noaa.gov)