Case Study 1: Auditing a Legacy Text-Output Pipeline

"The bytes were fine. Nobody could tell us what they meant." — a data manager, on inheriting a decade of simulation output

Executive Summary

A research group runs a working 2D diffusion solver that has, for years, dumped its temperature field to plain-text .dat files — one file per snapshot, a grid of numbers per file. The science is sound and the code is trusted, so nobody wants to touch the solver. But the output has become a liability: the archive is enormous, painfully slow to reload for analysis, and — the real problem — undocumented. A new collaborator sends her own data as CF-compliant NetCDF and asks the group to match it. This case study does the audit: we read the legacy writer to see exactly what it records and omits, quantify what the text format costs, read the collaborator's NetCDF file to understand the target, and write a small Fortran reader that validates it. We change nothing yet — the port itself is Case Study 2 — but by the end we know precisely what has to change and why.

Skills applied - Reading legacy text-I/O code and cataloguing its metadata gaps (§25.1) - Quantifying the size and speed cost of text at scale (§25.1) - Reading a NetCDF file's structure from its CDL header (§25.2, §25.4) - Mapping an ad-hoc format onto the NetCDF data model — dimensions, variables, attributes (§25.2) - Reading and validating a real NetCDF file from Fortran (§25.4), with CF metadata (§25.5)

Background

The solver writes a snapshot with a routine much like the write_field you built in Chapter 7 — assumed-shape array in, one text row per grid row out. The production version looks like this:

subroutine write_dat(u, filename)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp),         intent(in) :: u(:,:)
  character(len=*), intent(in) :: filename
  integer :: unit, i, j, nrow, ncol
  nrow = size(u, 1);  ncol = size(u, 2)
  open(newunit=unit, file=filename, status='replace', action='write')
  do i = 1, nrow
     write(unit, '(*(es24.16e3))') (u(i, j), j = 1, ncol)   ! full precision, so "safe"
  end do
  close(unit)
end subroutine write_dat

It is careful code: full-precision es24.16e3, so no digits are lost, and an unlimited-repeat format so it handles any width. And it is still the wrong format for an archive. The audit says why in numbers.

Phase 1 — Read the Writer: What Does It Record?

Before measuring anything, we catalogue what a .dat file actually contains versus what a reader needs. Walk the code and list every piece of information, present or absent:

The verdict is immediate: the file records the data and essentially none of the metadata. Everything a reader needs to interpret it lives in the group's collective memory and the solver's source. That is the "mute" failure of §25.1 in the wild — and it is why the collaborator, handed a .dat file, cannot use it.

Finding. The format is not lossy (full precision) but it is mute. The engineering fix is not "more decimals"; it is metadata that travels with the data.

Phase 2 — Measure the Cost

Now quantify the two costs that make the archive painful, so the case for change is a number, not an opinion. The production grid is $2048 \times 2048$; the run saves 2,000 snapshots. A short program does the arithmetic (pure Fortran, so the output is exact):

program audit_cost
  use, intrinsic :: iso_fortran_env, only: dp => real64, int64
  implicit none
  integer(int64), parameter :: nx = 2048, ny = 2048, n_snaps = 2000
  integer(int64) :: values, bin_bytes, txt_bytes
  real(dp) :: bin_gb, txt_gb

  values    = nx * ny                          ! per snapshot
  bin_bytes = values * 8_int64 * n_snaps       ! real64 binary
  txt_bytes = values * 25_int64 * n_snaps      ! ~25 chars per value (es24.16e3 + sep)
  bin_gb    = real(bin_bytes, dp) / 1.0e9_dp
  txt_gb    = real(txt_bytes, dp) / 1.0e9_dp

  print '(a, i0)',   'values per snapshot : ', values
  print '(a, f0.1)', 'binary archive (GB) : ', bin_gb
  print '(a, f0.1)', 'text   archive (GB) : ', txt_gb
end program audit_cost

values per snapshot : 4194304
binary archive (GB) : 67.1
text   archive (GB) : 209.7

Hand-check the magnitudes: $2048^2 = 4{,}194{,}304$ values; at 8 bytes over 2,000 snapshots that is $\approx 6.71 \times 10^{10}$ bytes $= 67.1$ GB binary, and at ~25 characters, $\approx 2.10 \times 10^{11}$ bytes $= 209.7$ GB of text. The text archive is roughly three times the size — 140 GB of pure overhead — and every byte of it was produced by a decimal conversion the CPU had to perform, and must perform again in reverse every time the data is reloaded.

⚡ The speed cost, illustratively. Reloading 200 GB of text means parsing ~8 billion numbers from decimal. At a Tier-2 illustrative rate of tens of millions of parses per second, that is minutes to tens of minutes of pure conversion per full-archive pass — before any analysis. The same data as compressed NetCDF is a fraction of the size and reads with no per-value conversion. (Order-of-magnitude figures, not a benchmark.)

Phase 3 — Read the Target: The Collaborator's NetCDF

The collaborator's file, reference.nc, is what the group must match. We inspect it with ncdump -h (header only) to learn the target structure — reading, not writing:

netcdf reference {
dimensions:
        x = 2048 ;
        y = 2048 ;
variables:
        double temperature(y, x) ;
                temperature:units = "K" ;
                temperature:long_name = "temperature" ;
                temperature:_FillValue = -999. ;
        double x(x) ;
                x:units = "m" ;
                x:axis = "X" ;
        double y(y) ;
                y:units = "m" ;
                y:axis = "Y" ;
// global attributes:
                :title = "diffusion benchmark" ;
                :Conventions = "CF-1.11" ;
                :history = "2026-05-02: created by diffuse v2.1" ;
}

Everything the .dat file omitted is here, in the file itself: the axes are named and located in metres (coordinate variables x(x) and y(y)), the field carries its units and a fill value, and the global attributes record the title, the convention, and the provenance. This is what "self-describing and reproducible" looks like as an artifact.

Phase 4 — Map the Ad-hoc Format onto the Model

With source and target in hand, the port becomes a mapping problem. Each piece of the group's world lands somewhere in the NetCDF data model:

The mapping exposes the one input the solver must start passing through: the grid spacing and units it already knows but currently throws away. Nothing about the physics changes; the port is about preserving metadata the code already has. That is the whole insight of the audit.

Phase 5 — Validate the Target File from Fortran

Finally, confirm the group can read the collaborator's format — a reader is the other half of any format decision. This program opens reference.nc, discovers its shape from the file, reads the field, checks the units attribute is what the group expects, and reports a sanity statistic. (Requires NetCDF installed; not run here.)

program validate_reference
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use netcdf
  implicit none
  real(dp), allocatable :: temperature(:,:)
  character(len=32) :: units
  integer :: ncid, varid, xid, yid, nx, ny

  call check( nf90_open('reference.nc', NF90_NOWRITE, ncid) )
  call check( nf90_inq_dimid(ncid, 'x', xid) )
  call check( nf90_inq_dimid(ncid, 'y', yid) )
  call check( nf90_inquire_dimension(ncid, xid, len=nx) )
  call check( nf90_inquire_dimension(ncid, yid, len=ny) )
  allocate(temperature(nx, ny))
  call check( nf90_inq_varid(ncid, 'temperature', varid) )
  call check( nf90_get_var(ncid, varid, temperature) )
  units = ''
  call check( nf90_get_att(ncid, varid, 'units', units) )   ! read the attribute
  call check( nf90_close(ncid) )

  print '(a, i0, a, i0)', 'grid  : ', nx, ' x ', ny
  print '(a, a)',         'units : ', trim(units)
  print '(a, f0.2, a, f0.2)', 'range : ', minval(temperature), ' to ', maxval(temperature)
  if (trim(units) /= 'K') print '(a)', 'WARNING: units are not Kelvin!'
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 validate_reference

The structure of the expected output is fixed even though the numbers depend on the file:

grid  : 2048 x 2048
units : K
range : <min> to <max>

The units check is the point: because the file declares its units, the reader can verify them and warn if a file arrives in Celsius by mistake. That safety is impossible with the .dat files, which never said what unit they were in — the very failure this audit set out to document.

Sanity check. A correct read reports the same $2048 \times 2048$ shape the header advertised, units = K, and a physically plausible temperature range. A shape mismatch would mean a corrupt file; a units mismatch would fire the warning — both caught because the file describes itself.

Discussion Questions

1. The legacy writer used full-precision es24.16e3, so it is not lossy. Why is "not lossy" insufficient for an archive? Which of §25.1's four failures does full precision not fix?
2. The grid spacing and units the port needs are values the solver already computes and then discards. What does that tell you about the cost of the port versus the cost of the science?
3. The collaborator's file encodes axes as coordinate variables in metres. What can she do with her data that the group currently cannot do with theirs, purely because of that metadata?
4. The reader validates units == 'K' and warns otherwise. Name two other CF attributes a defensive reader might check before trusting a file, and what each protects against.

Your Turn: Extensions

• Option A (analysis). Extend validate_reference to also read and print the global history and Conventions attributes, and to error stop if Conventions does not start with "CF-". What does requiring a convention buy a downstream pipeline?
• Option B (measurement). Reproduce audit_cost for your target grid and snapshot count, then add a fourth line: the compressed NetCDF size at an assumed 3× deflate ratio. How does the compressed archive compare to the original text archive — one number that would justify the port to a skeptical PI.
• Option C (mapping). Write the mapping table of Phase 4 as a checklist for your code's output. Which metadata does your solver already know and discard? That list is your port's scope.

Key Takeaways

• A working format can still be a broken archive. The legacy .dat files were correct to the last bit and useless to a stranger, because they were mute: no shape, no units, no provenance. Full precision fixes lossiness; it does nothing for muteness.
• Quantify the cost before you argue for change. 210 GB of text versus 67 GB of binary — a measured 3×, plus the reload conversion tax — turns "we should modernize the output" from a preference into a decision.
• The target teaches the port. Reading the collaborator's CDL header showed exactly which metadata was missing and where each piece belongs in the NetCDF model. Most of it — grid spacing, units — the solver already knows and throws away.
• A self-describing file can be validated, not just read. Because reference.nc declares its units, a reader can check them and refuse a wrong file. That verification is the reproducibility payoff, and it is impossible with an ad-hoc text dump. The build — actually writing the CF NetCDF — is Case Study 2.