Chapter 25: Scientific Data Formats — NetCDF, HDF5, and Managing Large Datasets

"Data is a precious thing and will last longer than the systems themselves." — Tim Berners-Lee

Overview

Your solver works. After Chapter 24 it marches a temperature field forward in time on a real grid, and after Chapter 7 it can write that field to a file. So run it in earnest: a $1000 \times 1000$ plate, ten thousand timesteps, a snapshot saved every hundred steps. You come back in the morning to a directory full of text files totalling tens of gigabytes — slow to write, slower to reload, wrong in the last decimal place, and utterly silent about what they contain. A year from now, or on a colleague's machine, they are a wall of numbers with no record of their grid spacing, their units, or even whether the rows run north–south or east–west. The computation was the easy part. The data is now the problem.

This is not a beginner's problem; it is the problem, and every field that computes at scale has converged on the same answer. The climate and weather communities store their data in NetCDF. The physics, astronomy, and general large-simulation communities store theirs in HDF5. Both are self-describing, portable, binary formats: they pack the numbers efficiently as raw bytes, and they store — inside the same file — the shape, the type, the byte order, and whatever metadata you attach, so the file explains itself to any program on any machine. They are the difference between a private scratch file and a scientific record that outlives the code that wrote it. Learning to read and write them is not optional polish; it is part of being employable in computational science, and it is what this chapter teaches.

In this chapter, you will learn to:

• Say precisely why text files, and even raw Fortran binary, do not scale — in size, in speed, in precision, and in the metadata they throw away.
• Write a labelled 2D field to a NetCDF file from Fortran, defining dimensions, variables, and attributes through the nf90_* interface and checking every return status.
• Write a chunked, compressed dataset to an HDF5 file through the h5*_f interface, and organize data into a hierarchy of groups.
• Read a variable back from either format, and choose between NetCDF and HDF5 for a real job with clear reasons rather than habit.
• Attach CF-convention metadata so your output is reproducible, tool-readable, and self-documenting — the standard that makes a NetCDF file a scientific artifact.

We will be honest about one thing throughout: these are large libraries, and no code in this book is ever executed (the promise from Chapter 1). The NetCDF and HDF5 examples require their libraries to be installed and linked, which we cannot do here, so we hand-trace the pure-Fortran logic and, for the library calls, describe the file each program produces rather than claim an exact line of terminal output. Every API name is real; where a signature is one we would want a reviewer to double-check against an installed version, we say so.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — this is your daily bread if you touch climate, ocean, astrophysics, or any large simulation. Read §25.2 (NetCDF) and §25.5 (CF conventions) closely; they are how your community shares data. Skim the HDF5 internals in §25.3 unless you own the output format. - 📖 Standard — read straight through. This chapter is the scalable successor to the text I/O of Chapter 7. - 🔧 Legacy — old codes write raw Fortran unformatted files (§7.5), which do not travel between machines or compilers. §25.1 explains exactly why, and NetCDF/HDF5 are the portable fix; wrapping a legacy code's output in one of them is a common, high-value modernization. - ⚡ HPC — §25.3 (HDF5, chunking, compression) is the one that matters at scale, and the note on parallel I/O points at the MPI-IO of Chapter 34. Text I/O is a classic hidden bottleneck; this chapter removes it.

25.1 Why Text Files Don't Scale

You already know the two ways Fortran writes a number, because Chapter 7 taught both. Formatted (text) output converts each value to human-readable decimal characters. Unformatted (binary) output writes the value's exact internal bytes with no conversion. Text is the natural first choice — you can open the file in an editor — and for a few thousand numbers it is completely fine. At scientific scale it fails on four separate counts, and it is worth naming all four, because each points at something the self-describing formats fix.

Text is bulky. A real(dp) occupies exactly eight bytes in memory. Written at enough precision to recover it exactly — roughly seventeen significant digits, say the format es24.16e3 — it becomes about two dozen characters of text once you count a separator. That is a threefold inflation before you have written a single byte of actual metadata. Let us make it concrete with a small program that just counts bytes; it is pure arithmetic, so we can compute its output by hand exactly:

program storage_estimate
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: nx = 100, ny = 100      ! one modest snapshot
  integer, parameter :: text_width = 24         ! chars per value, incl. separator
  integer :: n_values, bits_per, binary_bytes, text_bytes
  real(dp) :: ratio

  n_values     = nx * ny
  bits_per     = storage_size(1.0_dp)           ! bits in one real(dp): 64
  binary_bytes = n_values * (bits_per / 8)      ! exact internal size
  text_bytes   = n_values * text_width          ! decimal-text size

  ratio = real(text_bytes, dp) / real(binary_bytes, dp)

  print '(a, i0)',      'values         : ', n_values
  print '(a, i0)',      'binary  bytes  : ', binary_bytes
  print '(a, i0)',      'text    bytes  : ', text_bytes
  print '(a, f0.2)',    'text / binary  : ', ratio
end program storage_estimate

$ gfortran -std=f2018 -Wall -O2 example-01-storage-estimate.f90 -o est && ./est
values         : 10000
binary  bytes  : 80000
text    bytes  : 240000
text / binary  : 3.00

Ten thousand values — a small $100 \times 100$ plate — is 80 kB of binary and 240 kB of text. Scale it to that overnight run: a $1000\times1000$ field is 8 MB binary per snapshot, 24 MB as text, and at a thousand snapshots you have written 8 GB versus 24 GB. The factor of three is real disk and real transfer time, and it is the least of text's problems.

Text is slow. Converting a binary double to decimal digits, and parsing decimal digits back into a double, is genuine computation — division, rounding, digit extraction — performed once per value in each direction. For a billion values it is a billion conversions each way. A profiler in Chapter 28 will happily show a naïve solver spending more time formatting output than computing physics.

Text is lossy. This one is subtle and dangerous. Write x with f8.2 and you keep two decimal digits; read it back and you have thrown away the rest of the double's ~15 significant figures. Even "full-precision" text is only exact if you use a round-trip-safe format and the library rounds correctly in both directions. Binary has no such worry: the bytes you write are the bytes you read.

Text is mute. This is the deepest failure, and no amount of precision fixes it. A text file of numbers carries no record of what the numbers are. Is this a $100\times100$ grid or a $10000$-element vector? Are the values Kelvin or Celsius? Which index is $x$? What was alpha? The knowledge lives in your head and in the source code, and the day either is gone, the file is noise.

💡 Intuition: A raw data file is a photograph with no caption. You know what it shows because you were there when it was taken; a stranger, or you in five years, sees only pixels. A self-describing format is a photograph with the caption, date, location, and camera settings printed on the back — the file carries its own explanation, so it means the same thing to everyone who opens it.

Now, Chapter 7 already fixed the first three problems: unformatted and stream I/O write raw bytes, so they are compact, fast, and exact. Why isn't that the end of the story? Because raw binary fixes size, speed, and precision while making the muteness worse, and adds a new failure of its own:

Definition (self-describing format). A file format that stores, alongside the data, enough metadata to interpret the data without any external information: the shape and rank of each array, its numeric type and byte order (endianness), variable and dimension names, and arbitrary user attributes. A program that has never seen your code can open a self-describing file and discover what is inside it. NetCDF and HDF5 are the two dominant self-describing formats in scientific computing; a raw unformatted or stream file is the opposite — it is exactly as mute as the bytes it contains.

A raw stream file of a million doubles is a million anonymous bytes. Its shape is not recorded (was it $1000\times1000$ or $500\times2000$?), its type is not recorded (real32 or real64?), and — the classic portability trap — its byte order is not recorded. Write it on a little-endian laptop, read it on a big-endian machine, and every number is scrambled. Traditional Fortran unformatted-sequential files add compiler-specific record markers on top, so a file one compiler writes another may not even parse. Raw binary is a private note to yourself on one machine; it is not an archive.

⚡ Performance Note: The self-describing formats give you the speed and size of binary and the portability and metadata of a real archive — and then one thing more that raw binary cannot: transparent compression. Scientific fields are usually smooth, so their bytes are highly compressible; a temperature field that is 8 MB raw may shrink to 1–3 MB on disk with no change to your reading code, because the library compresses on write and decompresses on read. At the petabyte scale of a climate archive, that ratio is the difference between affordable and impossible. We meet the mechanism — chunking and deflate — in §25.3. (Compression ratios are illustrative; smooth fields compress well, noisy ones barely at all.)

The plan for the rest of the chapter follows directly. NetCDF (§25.2) is the self-describing format the geoscience world standardized on: simplest model, best conventions. HDF5 (§25.3) is the more general, more powerful format underneath it, with hierarchy, chunking, and compression exposed directly. We read both back and compare them honestly (§25.4), and we finish with the metadata conventions (§25.5) that turn a NetCDF file from "correct bytes" into "reproducible science." These libraries were named in the ecosystem tour of Chapter 16; now we put them to work.

25.2 NetCDF: The Self-Describing Standard for Geoscience

Definition (NetCDF). NetCDF — Network Common Data Form — is a self-describing, portable, binary data format and library for array-oriented scientific data, created and maintained by Unidata (part of the U.S. university consortium UCAR). Its data model is deliberately small: a file holds named dimensions, named variables (arrays whose shape is given by dimensions), and attributes (metadata key/value pairs attached to a variable or to the file as a whole). It is the de facto standard for climate, weather, and ocean data, and it is understood out of the box by a large ecosystem of tools.

The NetCDF data model is the whole thing, and it is simpler than it sounds. A dimension is a named length: x = 100, y = 100, time = unlimited. A variable is a named array declared over some dimensions and some element type: temperature of type double over (x, y). An attribute is a labelled value hung on a variable (temperature:units = "K") or on the file itself (a global attribute like :title = "heat solver run"). That is essentially all there is — dimensions, variables, attributes — and its simplicity is exactly why the format is so interoperable.

📜 From History: NetCDF came out of Unidata around 1988–1990 to solve a mundane, universal problem: atmospheric scientists at different universities, on different hardware (and in that era, genuinely different byte orders), needed to exchange gridded data without shipping the code to read it. The format was designed portable and self-describing from day one. In 2008, NetCDF version 4 was re-engineered to store its data inside an HDF5 file — so modern NetCDF is, underneath, HDF5 with a simpler model and stronger conventions bolted on top. The two formats we are learning are cousins, not rivals.

The Fortran-90 interface

Fortran talks to NetCDF through the netcdf module, whose procedures are all named nf90_*. There is one idiom to internalize before anything else: every nf90_* function returns an integer status code, which is nf90_noerr (zero) on success and something else on failure. Ignoring it is the cardinal sin of NetCDF programming, because a failed call that goes unchecked corrupts everything after it silently. The universal pattern is a tiny helper that checks the status and stops loudly:

subroutine check(status)
  use netcdf, only: nf90_noerr, nf90_strerror
  integer, intent(in) :: status
  if (status /= nf90_noerr) then
     print '(a)', 'NetCDF error: ' // trim(nf90_strerror(status))
     error stop 1
  end if
end subroutine check

nf90_strerror turns a status code into a human-readable message — the same iostat/iomsg spirit from Chapter 7, and the error stop from Chapter 13, applied to a library. Wrap every call in check(...) and a failure becomes a clear message at the exact failing line.

Writing a NetCDF file proceeds in two phases, a design that surprises newcomers but makes sense. In define mode you declare the structure — dimensions, variables, attributes — without writing any data. Then nf90_enddef switches to data mode, where you write the actual numbers. (Classic NetCDF requires this split so it can lay out the file header; NetCDF-4 is more relaxed, but writing define-then-data keeps your code portable across both.) Here is a complete program that writes a 2D temperature field with named axes and units:

program netcdf_write
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use netcdf
  implicit none

  integer, parameter :: nx = 4, ny = 3
  real(dp) :: temperature(nx, ny)
  integer  :: ncid, x_dimid, y_dimid, temp_varid
  integer  :: i, j

  ! A field we can describe exactly: temperature(i,j) = 10*i + j.
  do j = 1, ny
     do i = 1, nx
        temperature(i, j) = 10.0_dp * i + real(j, dp)
     end do
  end do

  ! --- DEFINE MODE: declare structure, no data yet ---
  call check( nf90_create('heat.nc', NF90_CLOBBER, ncid) )      ! create/overwrite
  call check( nf90_def_dim(ncid, 'x', nx, x_dimid) )            ! a dimension of length nx
  call check( nf90_def_dim(ncid, 'y', ny, y_dimid) )
  call check( nf90_def_var(ncid, 'temperature', NF90_DOUBLE, &
                           [x_dimid, y_dimid], temp_varid) )    ! a var over (x,y)
  call check( nf90_put_att(ncid, temp_varid, 'units', 'K') )    ! variable attribute
  call check( nf90_put_att(ncid, NF90_GLOBAL, 'title', &
                           '2D heat solver snapshot') )         ! global attribute
  call check( nf90_enddef(ncid) )                               ! --> DATA MODE

  ! --- DATA MODE: write the numbers ---
  call check( nf90_put_var(ncid, temp_varid, temperature) )
  call check( nf90_close(ncid) )                                ! flush and close

  print '(a)', 'wrote heat.nc'

contains
  subroutine check(status)
    integer, intent(in) :: status
    if (status /= nf90_noerr) then
       print '(a)', 'NetCDF error: ' // trim(nf90_strerror(status))
       error stop 1
    end if
  end subroutine check
end program netcdf_write

Because NetCDF is a separate library, the compile line must find its module and link its code. The netcdf-fortran package ships a helper, nf-config, that prints the right flags, so you never hard-code paths:

$ gfortran -std=f2018 -Wall example-02-netcdf-write.f90 `nf-config --fflags --flibs` -o ncwrite
$ ./ncwrite
wrote heat.nc

Note — requires NetCDF installed. This program compiles and runs only where netcdf-fortran is present; the back-ticked nf-config supplies the include path (--fflags) and the libraries (--flibs). Without the library the compiler reports it cannot find the netcdf module — a build error, not a bug in the code. We did not run it; we reason about the file it creates.

We cannot show you terminal output for the data itself, but we can show you the file's structure, because that is fully determined by the code above. The standard way to inspect a NetCDF file is the ncdump command-line tool, which prints the file's CDL (a text description of its structure). Running ncdump -h heat.nc (the -h shows the header — structure and metadata, not the data values) would print:

netcdf heat {
dimensions:
        x = 4 ;
        y = 3 ;
variables:
        double temperature(y, x) ;
                temperature:units = "K" ;

// global attributes:
                :title = "2D heat solver snapshot" ;
}

Read that carefully, because it demonstrates the whole point of the format: the file tells you what it holds. It has two dimensions with their sizes, one double-precision variable named temperature carrying a units attribute of "K", and a global title. A stranger's program — or a plotting tool, or you next year — learns all of this from the file alone.

⚠️ Common Pitfall — dimension order looks reversed, and it is not a bug. Notice that ncdump printed temperature(y, x) even though our Fortran declared the variable over [x_dimid, y_dimid]. NetCDF's tools and its C heritage report dimensions in row-major (C) order, slowest-varying first, whereas Fortran stores column-major, fastest-varying first (the layout from Chapter 5). The library handles the translation: your temperature(i, j) with i fastest is written correctly, and ncdump simply lists the axes in the opposite convention. Do not "fix" this by reversing your nf90_def_var dimension list — that would actually transpose your data. Know which convention a tool is using, and trust the library to bridge them.

The nf90_put_var call wrote the whole array in one statement, exactly as an unformatted write(u) array would (§7.5) — but now the bytes are wrapped in a portable, labelled container. And that is the essence of NetCDF from Fortran: define the shape and metadata, end define mode, write the data, close. Six kinds of call — create, def_dim, def_var, put_att, put_var, close — cover the great majority of writing.

🔄 Check Your Understanding 1. What are the three kinds of thing in the NetCDF data model, and what is each for? 2. Why does every nf90_* call return an integer, and what should you do with it? 3. Your code declares a variable over [x_dimid, y_dimid] but ncdump shows it as var(y, x). Is your data transposed? Explain.

Answers

1. Dimensions (named lengths, e.g. x = 100), variables (named arrays declared over dimensions, with an element type), and attributes (metadata key/value pairs attached to a variable or, with NF90_GLOBAL, to the whole file). Together they make the file self-describing.
2. The integer is a status code (nf90_noerr == 0 on success). NetCDF reports failures through return values, not exceptions or crashes, so you must check every one — the idiomatic way is a check wrapper that prints nf90_strerror(status) and error stops. An unchecked failure corrupts later calls silently.
3. No. ncdump lists axes in row-major (C) order, slowest-first; Fortran is column-major, fastest-first. The library maps between them, so your temperature(i, j) is stored correctly. The reversed display is a convention difference, not a transposition — reversing your dimension list would cause a real one.

25.3 HDF5: Hierarchy, Chunking, and Compression

NetCDF's flat model — one bag of variables and dimensions — is a virtue for interoperability and a limitation for complex output. A large simulation may want to group its data: this run's fields here, that run's there; the mesh in one place, the solution in another; a nested tree of cases. For that, and for the large end of the scale, the community reaches for HDF5.

Definition (HDF5). HDF5 — Hierarchical Data Format, version 5 — is a self-describing, portable, binary format and library for large and complex scientific data, maintained by The HDF Group. Its model is a filesystem inside a file: a file is a tree of groups (like directories) that contain datasets (like files — typed, multidimensional arrays) and attributes (metadata on any group or dataset). It is the standard for large simulation output across physics, astronomy, and engineering, and it is the storage layer beneath NetCDF-4.

The two new nouns are the payload. A dataset is HDF5's array — the analogue of a NetCDF variable — a typed, multidimensional block of data with its own shape. A group is a container that holds datasets and other groups, named with filesystem-style paths: /fields/temperature, /mesh/coordinates. Every file has a root group /, and you build a tree beneath it. This hierarchy is the headline difference from NetCDF.

🚪 Threshold Concept. A self-describing file is not a stream of bytes you parse; it is a small, portable database you query. Once you stop thinking "I write bytes and remember what they mean" and start thinking "I store named, typed, documented objects and the file remembers what they mean," the whole subject reorganizes. The data is decoupled from the program that wrote it — any language, any machine, any tool that speaks the format can open it and discover its contents. That decoupling is the entire reason these formats exist, and it is why a thirty-year-old NetCDF file still reads perfectly today while the raw binary from the same era is often unrecoverable.

Before the code, the two ideas that make HDF5 fast and small at scale — and the terms this chapter owns:

Definition (chunking and compression). By default a dataset is stored contiguously: the whole array as one block in the file, in order. Chunking instead stores the array as a grid of fixed-size rectangular tiles (chunks), each written independently — for a $1000\times1000$ field you might use $100\times100$ chunks. Chunking enables two things a contiguous layout cannot: reading a small sub-region without loading the whole array (you fetch only the overlapping chunks), and compression — shrinking the stored bytes with a lossless algorithm applied per chunk. The most common HDF5 compression filter is deflate (gzip), controlled by a level from 1 (fast, less shrinkage) to 9 (slow, most shrinkage). Compression requires chunking, because the filter runs on one chunk at a time.

The Fortran interface

HDF5's Fortran API lives in the hdf5 module and uses the h5*_f naming scheme. It differs from NetCDF's in three ways worth stating up front, because they shape every program:

1. Everything is a subroutine, not a function, and the last argument is an integer error flag, conventionally hdferr (0 on success, negative on failure). NetCDF returns its status; HDF5 hands it back through the final argument.
2. You must initialize and finalize the library with h5open_f(hdferr) at the start and h5close_f(hdferr) at the end. Forget h5open_f and nothing works.
3. Objects are opaque handles of special integer kinds. File, dataspace, dataset, and property-list identifiers are integer(hid_t); array sizes are integer(hsize_t). Both kinds come from the hdf5 module. And you must close every object you open — file, dataspace, dataset, property list — or you leak handles.

Here is the same $4\times3$ field written to HDF5, this time as a chunked, deflate-compressed dataset — so the example teaches the format's signature feature, not just its plumbing:

program hdf5_write
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use hdf5
  implicit none

  integer, parameter :: nx = 4, ny = 3
  real(dp) :: temperature(nx, ny)
  integer(hid_t)   :: file_id, space_id, dset_id, dcpl
  integer(hsize_t) :: dims(2), chunk_dims(2)
  integer :: hdferr, i, j

  do j = 1, ny
     do i = 1, nx
        temperature(i, j) = 10.0_dp * i + real(j, dp)
     end do
  end do
  dims       = [int(nx, hsize_t), int(ny, hsize_t)]
  chunk_dims = [int(2,  hsize_t), int(2,  hsize_t)]   ! 2x2 tiles (tiny, for the demo)

  call h5open_f(hdferr)                                       ! init the library
  call h5fcreate_f('heat.h5', H5F_ACC_TRUNC_F, file_id, hdferr)   ! create/truncate

  ! A dataspace describes the dataset's rank and shape.
  call h5screate_simple_f(2, dims, space_id, hdferr)

  ! A dataset-creation property list turns on chunking + gzip compression.
  call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl, hdferr)
  call h5pset_chunk_f(dcpl, 2, chunk_dims, hdferr)           ! chunk the storage
  call h5pset_deflate_f(dcpl, 6, hdferr)                     ! deflate level 6

  ! Create the dataset (note dcpl passed as the optional property-list argument).
  call h5dcreate_f(file_id, 'temperature', H5T_NATIVE_DOUBLE, space_id, &
                   dset_id, hdferr, dcpl)
  call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, temperature, dims, hdferr)

  ! Close everything we opened, innermost first.
  call h5pclose_f(dcpl, hdferr)
  call h5dclose_f(dset_id, hdferr)
  call h5sclose_f(space_id, hdferr)
  call h5fclose_f(file_id, hdferr)
  call h5close_f(hdferr)                                      ! finalize the library

  print '(a)', 'wrote heat.h5'
end program hdf5_write

HDF5 ships its own compiler wrapper, h5fc, which supplies the include path and libraries the way nf-config does for NetCDF:

$ h5fc -std=f2018 -Wall example-03-hdf5-write.f90 -o h5write
$ ./h5write
wrote heat.h5

Note — requires HDF5 installed. h5fc is the HDF5 Fortran wrapper around gfortran; equivalently you link by hand with -I$HDF5_DIR/include -lhdf5_fortran -lhdf5. As with NetCDF, we hand-trace the logic and describe the resulting file; we did not execute it.

The file this produces, inspected with the h5dump -H heat.h5 tool (-H = header only), would show its structure — a compressed, chunked dataset at the root, with its shape and datatype recorded:

HDF5 "heat.h5" {
GROUP "/" {
   DATASET "temperature" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 4, 3 ) / ( 4, 3 ) }
      STORAGE_LAYOUT { CHUNKED ( 2, 2 ) }
      FILTERS { COMPRESSION DEFLATE { LEVEL 6 } }
   }
}
}

Every property we set is now recorded in the file: the datatype is a little-endian IEEE 64-bit float (H5T_IEEE_F64LE — HDF5 recorded the byte order for us, the endianness problem solved), the shape is $(4,3)$, storage is chunked $2\times2$, and a deflate filter at level 6 is applied. A reader does not need to be told any of this; the file states it, and the library transparently decompresses on the way back in.

Organizing data into groups adds one pair of calls. To place the field under /fields instead of at the root, create the group and hand its identifier to h5dcreate_f in place of file_id:

integer(hid_t) :: grp_id
call h5gcreate_f(file_id, 'fields', grp_id, hdferr)        ! make /fields
call h5dcreate_f(grp_id, 'temperature', H5T_NATIVE_DOUBLE, space_id, &
                 dset_id, hdferr, dcpl)                    ! dataset at /fields/temperature
! ... write ...
call h5gclose_f(grp_id, hdferr)

That is the hierarchy in miniature: groups nest, datasets live inside them, and the path /fields/temperature addresses the array — a filesystem in a file.

⚠️ Common Pitfall — leaked identifiers. Every h5*create_f and h5*open_f returns a handle you must pass to the matching h5*close_f. Unlike a Fortran allocatable, an HDF5 object is not cleaned up automatically when its variable goes out of scope — the library holds it open. Forgetting h5sclose_f on a dataspace inside a loop over a thousand timesteps leaks a thousand handles and can exhaust the library's table. The discipline is mechanical: for every open, a close, innermost object first. (This is precisely the resource-management burden that allocatable spares you in native Fortran, from Chapter 9 — HDF5 handles are not Fortran objects, so the compiler cannot help.)

🔄 Check Your Understanding 1. In HDF5, what is the difference between a group and a dataset, and which one gives HDF5 its hierarchy? 2. Why must compression be paired with chunking? What unit does the deflate filter operate on? 3. A program opens a dataspace with h5screate_simple_f inside a loop that runs once per timestep. What must it do before the next iteration, and what fails if it does not?

Answers

1. A dataset is a typed, multidimensional array (like a NetCDF variable); a group is a container that holds datasets and other groups, addressed by a filesystem-style path (/fields/temperature). Groups give HDF5 its tree — the "filesystem in a file" — which NetCDF's flat model lacks.
2. The deflate filter runs on one chunk at a time, so there must be chunks: a contiguously stored dataset has no chunk boundaries for the filter to work along. Chunking is therefore a prerequisite for compression (and also enables partial, sub-region reads).
3. It must close the dataspace with h5sclose_f (and any per-iteration dataset/property-list handle) before looping. If it does not, each iteration leaks a handle; over many timesteps the leaked identifiers accumulate and can exhaust the library's table, because — unlike an allocatable — an HDF5 object is not reclaimed when its Fortran variable is reused.

25.4 Reading Both Back, and Choosing Between Them

Writing is half the job; a format you cannot read is useless. Both libraries read with the same shape as they write — open, locate the object by name, transfer the bytes, close — and both let you discover an object's shape from the file rather than hard-coding it, which is the self-describing promise cashed.

Reading NetCDF. Open the file for reading, ask for a variable's ID by name, query the dimensions to size your array, then read:

program netcdf_read
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use netcdf
  implicit none
  real(dp), allocatable :: temperature(:,:)
  integer :: ncid, varid, x_dimid, y_dimid, nx, ny

  call check( nf90_open('heat.nc', NF90_NOWRITE, ncid) )        ! open read-only
  call check( nf90_inq_dimid(ncid, 'x', x_dimid) )              ! find dimension x
  call check( nf90_inq_dimid(ncid, 'y', y_dimid) )
  call check( nf90_inquire_dimension(ncid, x_dimid, len=nx) )   ! learn its length
  call check( nf90_inquire_dimension(ncid, y_dimid, len=ny) )
  allocate(temperature(nx, ny))                                 ! size from the FILE
  call check( nf90_inq_varid(ncid, 'temperature', varid) )      ! find the variable
  call check( nf90_get_var(ncid, varid, temperature) )          ! read all of it
  call check( nf90_close(ncid) )

  print '(a, i0, a, i0)',   'read grid : ', nx, ' x ', ny
  print '(a, f8.2)',        'sum       : ', sum(temperature)
contains
  subroutine check(status)
    integer, intent(in) :: status
    if (status /= nf90_noerr) then
       print '(a)', 'NetCDF error: ' // trim(nf90_strerror(status))
       error stop 1
    end if
  end subroutine check
end program netcdf_read

We cannot run this — it needs the library and the heat.nc that netcdf_write produced — but we can predict its result exactly, because the values are determined by the writer. That program stored temperature(i, j) = 10*i + j on a $4\times3$ grid. Summed over all twelve elements:

$$ \sum_{i=1}^{4}\sum_{j=1}^{3}(10i + j) = \underbrace{3\sum_{i=1}^{4}10i}_{3(10+20+30+40)} + \underbrace{4\sum_{j=1}^{3}j}_{4(1+2+3)} = 300 + 24 = 324 . $$

So a correct round trip — write with netcdf_write, read with netcdf_read — would report:

read grid : 4 x 3
sum       :   324.00

The grid size was discovered from the file, not assumed, and the values returned bit-for-bit, so the hand-computed sum is exact. That is the round-trip guarantee a binary format gives and a lossy text format cannot.

Reading HDF5 follows the same arc with the h5*_f verbs — h5fopen_f with H5F_ACC_RDONLY_F, h5dopen_f to open the dataset by path, h5dread_f to transfer into a correctly-sized buffer, then close each handle and h5close_f. The one wrinkle is that to size the buffer from the file you query the dataset's dataspace (h5dget_space_f, then h5sget_simple_extent_dims_f) before reading — the analogue of NetCDF's dimension inquiry. Because that adds several handle-management calls without teaching a new idea, the full reader is in code/exercise-solutions.f90; the shape is identical to the writer, run in reverse.

🐛 Find the Bug. A colleague's NetCDF writer runs without error but produces a file ncdump reports as empty of data — all fill values. The relevant lines:

fortran call check( nf90_def_var(ncid, 'temperature', NF90_DOUBLE, [x_dimid, y_dimid], varid) ) call check( nf90_put_var(ncid, varid, temperature) ) ! write the data call check( nf90_enddef(ncid) ) ! end define mode call check( nf90_close(ncid) )

The bug is order: nf90_put_var is called while the file is still in define mode, before nf90_enddef. Data may only be written in data mode. In classic NetCDF the put_var returns a nonzero status (which check would catch — so the real-world version of this bug omits the check, the original sin); in some builds it is silently ignored. The fix is to move nf90_enddef above nf90_put_var: define all structure, end define mode, then write data. Define, then fill — never the reverse.

The two APIs side by side

The libraries do the same job with different spellings. This table is worth keeping; it is the practical core of the chapter:

Which one should you use?

Reach for NetCDF when your data is array-oriented and you will share it — especially in the earth sciences, where the CF conventions of §25.5 make a NetCDF file instantly readable by the whole community's toolchain (ParaView, Panoply, xarray, cdo, ncview). Its small model is a feature: less to get wrong, more that just works. Reach for HDF5 when you need hierarchy (many groups, nested cases, mesh-plus-solution), maximum control over chunking and compression, or you are producing very large simulation output whose format you own. And remember the two are not exclusive: NetCDF-4 is HDF5 underneath, so writing NetCDF-4 gives you HDF5's chunking and compression through NetCDF's friendlier model and stronger conventions — often the best of both. When in genuine doubt for shareable gridded data, default to NetCDF; when in doubt for a big private simulation archive, default to HDF5.

🐍 Python Comparison: These formats are the meeting point of the "Fortran computes, Python analyzes" workflow the book keeps returning to. The NetCDF your solver writes is read in one line by xarray (xarray.open_dataset('heat.nc')) or netCDF4, and the HDF5 by h5py (h5py.File('heat.h5')) — and because the files are self-describing, Python discovers the shapes, types, and your metadata automatically, with no agreement about byte layout of the kind raw stream files demanded in §7.5. Write the heavy field in fast Fortran, hand off a .nc file, and plot it in a notebook: the format is the clean interface between the two languages, complementing the in-memory hand-off via f2py from Chapter 15.

25.5 Metadata and CF Conventions for Reproducibility

A file can be self-describing in the mechanical sense — shapes and types recorded — and still be scientifically useless, because it does not say what the numbers mean. Are they Kelvin or Celsius? What are the physical coordinates of grid point $(i, j)$? When was the file made, by what code, from what input? Self-description is structure; reproducibility needs meaning. In the NetCDF world, meaning has a standard.

Definition (CF conventions). The CF (Climate and Forecast) conventions are a community standard for the metadata in a NetCDF file: which attributes to attach, and what their values must say, so that the data is unambiguous and machine-interpretable. CF specifies, among much else, that every variable carry a units attribute drawn from the UDUNITS system ("K", "m", "s"), an optional standard_name from a controlled vocabulary and a free-text long_name; that physical axes be given as coordinate variables; that missing data be flagged with _FillValue; and that the file carry global provenance attributes (title, institution, source, history, Conventions). A CF-compliant file is understood, without configuration, by the entire geoscience toolchain.

The keystone idea is the coordinate variable.

Definition (coordinate variable). A coordinate variable is a one-dimensional NetCDF variable with the same name as a dimension, holding the physical coordinate of each index along that axis. A dimension x of length 4 says only "there are four columns"; a coordinate variable x(x) = [0.0, 0.25, 0.5, 0.75] with units = "m" says where those columns are on the plate. Coordinate variables are what let a tool label axes in metres, place data on a map, or interpolate — they turn array indices into physical space.

Here is the temperature field again, now written with CF-style metadata: coordinate variables for the two axes with real spacing, physical units and a long_name on the field, and global provenance. The mechanics are the same nf90_* calls from §25.2 — CF is a discipline about which attributes to write, not new machinery:

! (excerpt — full program in code/, requires NetCDF installed)
call check( nf90_def_dim(ncid, 'x', nx, x_dimid) )
call check( nf90_def_dim(ncid, 'y', ny, y_dimid) )

! Coordinate variables: same name as the dimension, 1D over it.
call check( nf90_def_var(ncid, 'x', NF90_DOUBLE, x_dimid, x_varid) )
call check( nf90_put_att(ncid, x_varid, 'units', 'm') )
call check( nf90_put_att(ncid, x_varid, 'long_name', 'x coordinate') )
call check( nf90_put_att(ncid, x_varid, 'axis', 'X') )
call check( nf90_def_var(ncid, 'y', NF90_DOUBLE, y_dimid, y_varid) )
call check( nf90_put_att(ncid, y_varid, 'units', 'm') )
call check( nf90_put_att(ncid, y_varid, 'axis', 'Y') )

! The field, over (x, y), with physical units and a description.
call check( nf90_def_var(ncid, 'temperature', NF90_DOUBLE, [x_dimid, y_dimid], t_varid) )
call check( nf90_put_att(ncid, t_varid, 'units', 'K') )
call check( nf90_put_att(ncid, t_varid, 'long_name', 'plate temperature') )
call check( nf90_put_att(ncid, t_varid, '_FillValue', -999.0_dp) )

! Global provenance — who, what, when, and to which standard.
call check( nf90_put_att(ncid, NF90_GLOBAL, 'title', '2D heat solver snapshot') )
call check( nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Your Lab') )
call check( nf90_put_att(ncid, NF90_GLOBAL, 'source', 'heat-solver v1.0, explicit FD') )
call check( nf90_put_att(ncid, NF90_GLOBAL, 'history', '2026-07-22: created by heat_solver') )
call check( nf90_put_att(ncid, NF90_GLOBAL, 'Conventions', 'CF-1.11') )
call check( nf90_enddef(ncid) )

call check( nf90_put_var(ncid, x_varid, [(real(i-1, dp)*dx, i = 1, nx)]) )   ! 0, dx, 2dx, ...
call check( nf90_put_var(ncid, y_varid, [(real(j-1, dp)*dy, j = 1, ny)]) )
call check( nf90_put_var(ncid, t_varid, temperature) )

The resulting header, via ncdump -h, is now a complete scientific label:

netcdf heat_cf {
dimensions:
        x = 4 ;
        y = 3 ;
variables:
        double x(x) ;
                x:units = "m" ;
                x:long_name = "x coordinate" ;
                x:axis = "X" ;
        double y(y) ;
                y:units = "m" ;
                y:axis = "Y" ;
        double temperature(y, x) ;
                temperature:units = "K" ;
                temperature:long_name = "plate temperature" ;
                temperature:_FillValue = -999. ;

// global attributes:
                :title = "2D heat solver snapshot" ;
                :institution = "Your Lab" ;
                :source = "heat-solver v1.0, explicit FD" ;
                :history = "2026-07-22: created by heat_solver" ;
                :Conventions = "CF-1.11" ;
}

Nothing about the data changed — the same twelve temperatures — but the file is transformed. It now states its units, locates every value in physical space, flags its missing-data sentinel, names the code and date that made it, and declares which standard it follows. Hand it to ParaView and the axes come up in metres; hand it to xarray and the coordinates are attached automatically; hand it to a reviewer and it answers the questions a reviewer asks.

🔗 Connection: CF metadata is what makes the leap to visualization (Chapter 26) effortless — ParaView, VisIt, Panoply, and ncview read CF-compliant NetCDF directly and lay it out in physical space with no configuration, because the file told them how. It is also the concrete face of reproducibility, the software-engineering discipline of Chapter 37: the history and source attributes record how the data was made, so a result can be traced back to the code, inputs, and date that produced it. A figure in a paper should trace to a CF file, and that file should trace to a commit. Appendix H's library reference lists the NetCDF, HDF5, and CF documentation to keep at hand: Appendix H.

The reproducibility habit is cheap and it compounds. Writing units, source, and history costs three lines at write time and saves an afternoon of forensic guessing every time someone — including future you — opens the file. The history attribute in particular is conventionally append-only: each tool that modifies the data adds a line, so the attribute becomes a provenance log of everything that ever happened to the file. That is what "the data outlives the systems" looks like in practice: a file that still explains itself, and its own history, long after the program that wrote it is gone.

⚠️ Common Pitfall — inventing standard_name values. CF's standard_name is drawn from a controlled vocabulary — a published table of exact strings like air_temperature or sea_water_salinity. You cannot make one up; a name not in the table is not CF-compliant and tools may reject or ignore it. When no table entry fits (as for a generic teaching plate), omit standard_name and use the free-text long_name instead, which has no vocabulary restriction. Better a correct long_name than a fabricated standard_name.

🔄 Check Your Understanding 1. What makes a NetCDF variable a coordinate variable, and what can a tool do with the data once it has one that it cannot do with a bare dimension? 2. Which CF global attributes record a file's provenance, and why does the history attribute earn its reproducibility value by being append-only? 3. You want to label a temperature field. Which of units, long_name, and standard_name is unrestricted free text, and which is drawn from a controlled vocabulary you must not invent?

Answers

1. A coordinate variable is a 1D variable with the same name as a dimension (x(x)), holding the physical position of each index along that axis. With it, a tool can label axes in real units, place the data in physical space (e.g. on a map), and interpolate — none of which a bare dimension (just a length) permits.
2. title, institution, source, and history (plus Conventions) carry provenance. history is append-only — each tool that touches the data adds a line — so the attribute becomes a running log of everything that was ever done to the file, letting a result be traced back to the code, inputs, and steps that produced it.
3. long_name is unrestricted free text; units must come from UDUNITS; standard_name is drawn from the CF controlled vocabulary and must not be invented — omit it and use long_name when no table entry fits.

Project Checkpoint

Through Chapter 7 your solver's write_field wrote the temperature field to a text file — perfect for a $4\times4$ demonstration, ruinous for the $1000\times1000$, ten-thousand-step production run this part exists to enable. This checkpoint gives heat_io a scalable sibling: write_field_netcdf, which writes the field of the field_t type from Chapter 9 to a compact, self-describing, CF-labelled NetCDF file. It is the optional-but-recommended alternative to text output for large runs; the text writer stays for small cases and quick eyeballing.

Recall field_t (Ch. 9): it bundles nx, ny, dx, dy and the allocatable field u(:,:). The new routine takes one field_t, plus the timestep number and simulation time for the metadata, and produces one .nc file per snapshot:

subroutine write_field_netcdf(field, filename, step, time)
  use netcdf
  type(field_t),    intent(in) :: field           ! the Ch.9 derived type
  character(len=*), intent(in) :: filename
  integer,          intent(in) :: step
  real(dp),         intent(in) :: time
  integer :: ncid, x_dimid, y_dimid, t_varid

  call check( nf90_create(filename, NF90_NETCDF4, ncid) )    ! HDF5-backed NetCDF-4
  call check( nf90_def_dim(ncid, 'x', field%nx, x_dimid) )
  call check( nf90_def_dim(ncid, 'y', field%ny, y_dimid) )
  call check( nf90_def_var(ncid, 'temperature', NF90_DOUBLE, &
                           [x_dimid, y_dimid], t_varid) )
  call check( nf90_put_att(ncid, t_varid, 'units', 'K') )
  call check( nf90_put_att(ncid, t_varid, 'long_name', 'plate temperature') )
  call check( nf90_put_att(ncid, NF90_GLOBAL, 'Conventions', 'CF-1.11') )
  call check( nf90_put_att(ncid, NF90_GLOBAL, 'source', 'heat-solver, explicit FD') )
  call check( nf90_put_att(ncid, NF90_GLOBAL, 'step', step) )       ! which snapshot
  call check( nf90_put_att(ncid, NF90_GLOBAL, 'time', time) )       ! simulated seconds
  call check( nf90_enddef(ncid) )
  call check( nf90_put_var(ncid, t_varid, field%u) )               ! the whole field
  call check( nf90_close(ncid) )
end subroutine write_field_netcdf

Requires NetCDF installed; compile the solver with `nf-config --fflags --flibs` appended to the usual command, exactly as in §25.2. The library calls are shown as representative and are not run here; the surrounding module structure compiles as ordinary modern Fortran. The full version, with the check helper and a field_t definition so it builds on its own, is code/project-checkpoint.f90.

The design mirrors write_field's canonical signature — one bundled object plus a filename — so the two writers are interchangeable at the call site: a run switches between text and NetCDF output by swapping one line. Three lessons from this chapter are baked in. The output is self-describing (shape, type, and units travel with the data), it is CF-labelled (so Chapter 26 can open it in ParaView with no configuration), and using NF90_NETCDF4 means the field is stored in HDF5 underneath, so a one-line addition later — the chunksizes/deflate_level options on nf90_def_var — buys transparent compression for the big runs. This is the output format the Chapter 38 capstone writes and the figures in its "paper" are made from.

Summary

This chapter replaced the solver's text output with the self-describing, portable, compressible formats that scientific computing actually runs on.

The two things to memorize. First, self-describing means the file carries its own metadata — shape, type, byte order, units — so it means the same thing to every program and outlives the code that wrote it; that single property is why these formats exist and why raw binary is not enough. Second, check every library call: nf90_* through its returned status (a check wrapper), h5*_f through its trailing hdferr, and close every HDF5 handle you open — the library will not clean up after you the way allocatable does.

Spaced Review

Four questions reaching back to Chapter 7 (I/O) and Chapter 9 (derived types) — the two chapters this one stands on: Chapter 7's binary I/O is the foundation these formats refine, and Chapter 9's field_t is what the project checkpoint serializes.

1. (Ch. 7) In §7.5 you wrote a field with unformatted write(u) array and it round-tripped bit-for-bit, where text with f8.2 did not. Which of those two properties does NetCDF preserve, and what does NetCDF add that neither raw form has?
   

   Answer

   NetCDF preserves the **exact, bit-for-bit** round trip of unformatted binary (it stores raw bytes, not decimal text), so it is not lossy. What it adds over *both* text and raw binary is that it is **self-describing and portable**: the file records the array's shape, type, byte order, and your metadata, so a program on any machine can read it correctly — where a raw unformatted file records none of that and breaks across endianness and compilers.

2. (Ch. 7) The check subroutine in this chapter prints nf90_strerror(status) and calls error stop. Which Chapter 7 mechanism is this the library analogue of, and what are the two things it does that letting the call fail silently would not?
   

   Answer

   It is the library analogue of the **`iostat`/`iomsg`** pattern (§7.6): ask whether the operation succeeded rather than let it crash or corrupt silently. The two things it does are (1) turn a numeric status into a **human-readable message** (`nf90_strerror`, like `iomsg`), and (2) **stop at the failing line** with a nonzero exit code (`error stop`), instead of continuing past a corrupt file — so a failure is a clear report, not a mystery downstream.

3. (Ch. 9) The project's write_field_netcdf takes one type(field_t) argument instead of nx, ny, dx, dy, u as five separate ones. Name two concrete benefits of passing the derived type, per Chapter 9.
   

   Answer

   (1) The arguments **cannot be passed in the wrong order** — one object with named components replaces five loose scalars/arrays that are easy to transpose. (2) The interface is **stable and self-documenting**: `write_field_netcdf(field, ...)` says what it operates on, and adding a component to `field_t` (say `dz`) does not change the routine's signature. (Also acceptable: the field's `u` is an **allocatable component** that carries its own shape, so `field%nx`/`field%ny` and `size(field%u)` travel with the data.)

4. (Ch. 9) field_t holds u as real(dp), allocatable :: u(:,:). When you pass a field_t into write_field_netcdf with intent(in) and read field%u, why is there no risk of a dangling handle of the kind §25.3 warned about for HDF5 dataspaces?
   

   Answer

   An `allocatable` component is a **native Fortran object**: the compiler tracks its lifetime and deallocates it automatically when the `field_t` goes out of scope, and `intent(in)` guarantees the routine only reads it. An HDF5 dataspace or dataset is an **opaque library handle**, not a Fortran object, so the compiler cannot track it — you must `h5*close_f` it by hand. The hazard exists for HDF5 handles precisely because they are outside the language's automatic memory management, which `allocatable` (Ch. 9) provides for free.

What's Next

Your solver now writes compact, self-describing, CF-labelled data — but a .nc file full of numbers is still not a picture, and the whole point of a simulation is usually to see it. Chapter 26 closes Part VI by turning output into images: the VTK format that ParaView and VisIt animate, simpler formats for gnuplot and matplotlib, and how to write a time series your solver can produce one snapshot per step. The CF-compliant NetCDF you just learned to write is one input ParaView reads directly; VTK is the other, purpose-built for unstructured and structured visualization. Between them, your heat equation stops being a column of temperatures and starts being a movie of heat spreading across a plate — the payoff the whole book has been building toward.