umwm

UMWM Technical Reference

This file is a technical reference for the physics implemented in the University of Miami Wave Model (UMWM).

Implementation Map

Physics area Main implementation Purpose
Model domain and loop src/umwm_domain.f90 domain_type owns the model state, initializes it, and runs the time loop: source functions, source integration, advection, refraction, stress, Stokes drift, diagnostics, and output.
Output and output-time physics src/umwm_domain_output.f90 Integrated diagnostics, ocean-side stress and energy fluxes, gridded and spectrum NetCDF output.
Forcing state src/umwm_forcing.f90 forcing_type owns mutable wind, current, density, sea-ice, remapped, and interpolation-scratch fields, and the wave partitions of the last input file read.
Wave partitions src/umwm_partitions.f90 partitions_type with the wave partitions of an input file (forcing.waves), and the parametric spectra that rebuild the wave spectrum from them, for the initial state and the open boundaries.
Source functions src/umwm_source_functions.f90 Wind input, wave breaking, nonlinear downshifting, turbulence, and sea-ice attenuation.
Source integration src/umwm_physics.f90 Exponential source update, diagnostic spectral tail, physics time step limit.
Propagation, refraction, and frequency shift src/umwm_advection.f90 First-order upstream advection of the wave action in geographic, directional, and frequency space.
Stress src/umwm_stress.f90 Form/skin drag, drag coefficient, stress tail, and spectral stress integrals. The domain evaluates the atmospheric stress each time step, and the ocean-top and ocean-bottom momentum fluxes (ocean_stress in src/umwm_domain_output.f90) at gridded output times, or each time step when coupled with ESMF.
Stokes drift src/umwm_stokes.f90 stokes_type with wave-induced Stokes drift and e-folding depth.
Spectral grid src/umwm_spectrum.f90 spectrum_type, frequency and direction coordinates, and spectral spacings.
Dispersion src/umwm_dispersion.f90 Capillary-gravity dispersion relation, wavenumber solve, angular frequency, and group speed.
Initialization and precomputed factors src/umwm_domain.f90 Runtime allocation, spectral coordinates, wave kinematics, integration weights, dissipation constants, CFL limits.
Parallel environment and halo exchange src/umwm_env.f90, src/umwm_mpi.f90 Process rank and count; halo exchange and gathering of tile fields.
Command line src/umwm_cli.f90 Arguments of the umwm executable, umwm CONFIG [-o OUTPUT] [-i INPUT]: the configuration directory or main namelist to run, the input directory (CONFIG/input by default), the output directory (the current directory by default), which it creates, the copies of the namelists in it, and the preparation of the missing input files from a data source (forcing.source) with the Python package.
Data sources python/umwm/sources/ Preparation of the input files on the grid of the model from a reanalysis or forecast: the ERA5 source (era5.py), the regridding (regrid.py), the wave partitions (waves.py), and the files (prepare.py).
C API src/umwm_capi.f90 bind(c) procedures for foreign-language callers: the domain lifecycle through an opaque handle, and access to the domain fields by name as pointers into the model memory. Used by the Python package in python/.

The main program instantiates one domain_type from src/umwm_domain.f90 and calls its initialize, run, and finalize procedures, with the configuration given on the command line: a directory with main.nml, or a main namelist file. The spectrum.nml namelist is read from the directory of the main namelist, with its prefix if its name ends with _main.nml (umwmout_spectrum.nml next to umwmout_main.nml). The input files are read from the input directory, CONFIG/input or the directory given with -i, and the output and restart files are written to the output directory, given with -o or else the current directory, with copies of the namelists of the run, umwmout_main.nml and umwmout_spectrum.nml; a restart reads its restart file from the output directory. As the files the model writes are all named umwmout_*, umwmspc_* or umwmrst_*, the output directory may be the configuration or input directory. run loops over advance, which integrates one global time step dtg: the forcing update, the inner loop of source time steps, and the output due at the new time; advance and finished are public so that a driver, such as the Python package, can step the model and inspect it between steps. The domain is constructed from config_type, spectrum_type, grid_type, forcing_type, and stokes_type components, and owns the wave spectrum and all fields derived from it; there is no global model state.

The Python package in python/ (see python/README.md) drives the model through the C API of src/umwm_capi.f90, compiled into the shared library libumwm.so by make shared: a Python Domain holds a domain_type through an opaque handle, and the fields of the domain are exposed as NumPy arrays viewing the Fortran memory, so the time integration runs entirely in Fortran. config_type remains the source of forcing switches and constants, while forcing_type owns the mutable runtime forcing snapshots and remapped fields.

The source functions, integration, advection, and stress kernels operate on one spectral column (all frequencies and directions of a grid cell) with explicit-shape array arguments, so that compilers can vectorize them; the domain loops over cells and passes column slices of its state.

Within each global forcing/output time step dtg, UMWM repeatedly runs:

  1. Advance file-backed forcing snapshots, when not ESMF-coupled.
  2. Interpolate forcing fields.
  3. Compute Sin, Sds, Sice, turbulence and linear sink rates, and the largest stable physics time step of each cell.
  4. Find the dynamic physics time step dts, the minimum over all cells.
  5. Integrate source terms over dts, including Snl.
  6. Exchange MPI halos when enabled.
  7. Propagate wave action geographically.
  8. Refract wave action directionally and, with currents, shift it in frequency.
  9. Compute atmospheric stress and drag.
  10. At output times, compute Stokes drift, diagnostics, ocean stress, spectra, grid output, and restart output. When coupled with ESMF, ocean stress is computed each time step.

Steps 5 and 7 to 9 are pipelined over grid cells: a cell is propagated a fixed number of cells after its integration, at least the largest distance to a neighbor cell in the tile, so that the spectra of its neighbors are integrated and still in cache. Under MPI, the cells in the halo of neighboring tiles are integrated first, and the halo is exchanged before the pipeline starts. The halo is one grid row or column wide, and the model stops at startup if a tile is narrower than the halo of its neighbors.

State Variables and Conventions

UMWM predicts the surface elevation variance spectrum $E(\mathbf{x}, k, \phi, t)$ on a two-dimensional horizontal grid. The state is discretized in:

The main spectrum array is e(o,p,i).

Forcing fields are represented by type(forcing_type) from src/umwm_forcing.f90. It stores native-grid snapshots for backward and forward forcing time levels, interpolated native-grid scratch fields, and remapped one-dimensional fields consumed by source functions, advection/refraction, stress, diagnostics, and output. File-backed updates remain disabled under ESMF, where coupled forcing is expected to be supplied externally.

The spectral grid metadata is represented explicitly by type(spectrum_type) from src/umwm_spectrum.f90. domain % initialize constructs one spectrum object after namelist input is read, and the domain uses it for the spectral dimensions and spacings, the direction projection of the grid, restart I/O, and spectrum and grid output. The domain keeps default-real copies of the frequency and direction coordinates, f(:) and th(:), with their cosines and sines, for the model arrays and formulas.

The model uses mathematical direction convention: $\phi=0$ points in the positive $x$ direction and positive angles rotate counter-clockwise. Wind direction wdir uses the same convention.

The relationship between wave energy spectrum $E’$ and variance spectrum $E$ is

\[E'(\mathbf{x}, k, \phi) = \rho_w g E(\mathbf{x}, k, \phi).\]

For implementation, the variance spectrum is the prognostic quantity. Most integrated quantities use the discrete spectral weight

\[k \, dk \, d\phi.\]

In code, dth is $d\phi$, dwn(o,i) is $dk$, and kdk(o,i) is $k_o \, dk_o$.

Spectral Grid and Wave Kinematics

Frequencies are logarithmically spaced:

\[\Delta \ln f = \frac{\ln f_{\max} - \ln f_{\min}}{N_f - 1}, \qquad f_o = \exp\left(\ln f_{\min} + (o-1)\Delta\ln f\right).\]

This construction is implemented by spectrum_type using the DOMAIN namelist values om, pm, fmin, and fmax:

spectrum_type field Meaning
num_frequencies Number of frequency bins, initialized from om.
num_directions Number of direction bins, initialized from pm.
frequency_min, frequency_max Lowest and highest frequency bins, initialized from fmin and fmax.
frequency(:) Logarithmically spaced frequency coordinates.
direction(:) Uniformly spaced directional coordinates.
dlnf Log-frequency spacing.
dth Directional spacing.

Directions are uniformly spaced over a full circle:

\[\Delta\phi = \frac{2\pi}{N_\phi}, \qquad \phi_p = \left(p-\frac{N_\phi+1}{2}\right)\Delta\phi.\]

During initialization, the domain arrays f(:) and th(:) are assigned from spectrum % frequency and spectrum % direction.

The angular-frequency increment used to infer $dk$ is

\[d\omega_o = 2\pi \Delta \ln f \, f_o.\]

For every grid cell and frequency, the wavenumber is solved from the finite-depth capillary-gravity dispersion relation

\[\omega^2 = \left(gk + \frac{\sigma}{\rho_w} k^3\right)\tanh(kd), \qquad \omega = 2\pi f,\]

where $d$ is water depth and $\sigma$ is water surface tension (sfct). The pure elemental function wavenumber() in src/umwm_dispersion.f90 solves a nondimensional form by Newton iteration. The same module also provides angular_frequency() for the forward relation and group_speed() for the derivative used by the model.

Intrinsic phase speed and group speed are

\[c = \frac{\omega}{k},\] \[c_g = c\left[ \frac{1}{2} + \frac{\min(kd,20)}{\sinh(2\min(kd,20))} + \frac{\sigma k^2}{\rho_w g + \sigma k^2} \right].\]

The min(kd,20) cap is the implementation limit used in group_speed() to avoid hyperbolic overflow at large nondimensional depth.

The spectral wavenumber increment is approximated from the group speed:

\[dk_o = \frac{d\omega_o}{|c_{g,o}|}.\]

Precomputed factors include:

\[kdk = k\,dk,\qquad k3dk = k^3\,dk,\qquad k4 = k^4,\qquad fkovg = \frac{fk}{g}.\]

The half wavelength used by the wind input is

\[\frac{\lambda}{2} = \frac{|c|}{2f}.\]

The logarithmic wind-profile factor is capped at the field-log-layer top:

\[\log\left(\frac{\lambda/2}{z}\right) \rightarrow \log\left(\frac{20\,\mathrm{m}}{z}\right) \quad \text{when } \lambda/2 > 20\,\mathrm{m}.\]

Governing Equation

The model solves the wave action balance equation. The spectral grid is fixed in the intrinsic frequency $\sigma = 2\pi f$, and each cell has the wavenumbers $k(\sigma, d)$ of its depth. The action density per unit $\ln\sigma$ and $\phi$ is

\[A = \frac{E_\sigma}{\sigma}\,\frac{d\sigma}{d\ln\sigma} = E\,\frac{k}{c_g},\]

where $E(k,\phi)$ is the variance spectrum and $E_\sigma = E\,k\,dk/d\sigma = E\,k/c_g$ the variance per unit $\sigma$. In $(\mathbf{x}, \ln\sigma, \phi)$,

\[\frac{\partial A}{\partial t} + \nabla\cdot(\dot{\mathbf{x}}A) + \frac{\partial(\dot{\varsigma} A)}{\partial \ln\sigma} + \frac{\partial(\dot{\phi}A)}{\partial \phi} = \frac{k}{c_g}\sum_{i=1}^{n} S_i,\]

where $\dot{\varsigma} = d\ln\sigma/dt$, and $S_i$ are the source terms of $E$. The geographic propagation velocity is

\[\dot{\mathbf{x}} = \mathbf{c}_g + \mathbf{u},\]

where $\mathbf{u}=(u,v)$ is the current in the wave boundary layer. Directional-space advection is refraction, nonzero for spatially varying depth and/or inhomogeneous currents. Frequency-space advection is nonzero where a current flows across depth contours or varies along the waves; it neglects the change of depth in time.

Without currents, $\sigma$ is constant along the rays, and the action balance is the energy balance at a fixed frequency. Even then, the conserved density of a frequency bin is $E\,k/c_g$, not $E$: $k/c_g$ varies with depth at a fixed frequency, and its ratio between cells gives the shoaling of the waves. With currents, the waves exchange energy with the current, and the action, not the energy, is conserved; and the intrinsic frequency changes along the rays, with the absolute frequency $\sigma + \mathbf{k}\cdot\mathbf{u}$ conserved in a steady flow.

The factor $k/c_g$ (kovcg) of each cell and frequency is the same for all directions and fixed in time, so the prognostic variable remains $E$: the source terms act on $E$, refraction of $E$ is that of $A$, and propagation and the frequency shift convert between $E$ and $A$ by kovcg.

For the variance spectrum, the local source terms are written as signed rates or tendencies acting on $E$.

Source-Term Sign Conventions

The manual writes physical source terms as signed variance-spectrum tendencies. The implementation stores several of them as positive decay rates. This distinction is important.

Physical term Manual sign Implementation array Implementation meaning
Wind input/export positive or negative ssin(o,p,i) Signed rate multiplying e.
Breaking dissipation negative sds(o,p,i) Positive decay rate multiplying e.
Turbulent dissipation negative sdt(o,i) Positive decay rate multiplying e.
Viscous dissipation negative sdv(o,i) Positive decay rate multiplying e.
Bottom friction and percolation negative sbf(o,i) Positive combined bottom decay rate multiplying e.
Nonlinear downshifting source/sink snl(o,p,i) Absolute tendency, already multiplied by donor-bin variance.
Sea-ice attenuation negative when active sice(o,i) Signed rate, normally negative.

For prognostic bins, the source update uses

\[E_s^{n+1} = E^n \exp\left[ \Delta t \left( S_{in}^* - S_{ds}^* - S_{bottom}^* - S_{dt}^* - S_{dv}^* + S_{ice}^* \right) \right] + \Delta t\,S_{nl},\]

where stars denote the implementation rate arrays that multiply $E$. Snl is not inside the exponential because it is carried as an absolute bin-to-bin tendency.

Wind Input Sin

Wind input follows Jeffreys’ sheltering hypothesis. Define the wave-wind angle

\[\theta = \psi - \phi,\]

where $\psi$ is wind direction, and define the effective wind-wave relative speed

\[\Delta U = U_{\lambda/2}\cos\theta - c - u\cos\phi - v\sin\phi.\]

The physical source is

\[S_{in} = A_1 \Delta U |\Delta U| \frac{k\omega}{g} \frac{\rho_a}{\rho_w} E(k,\phi).\]

The corresponding implementation rate is

\[S_{in}^* = A_1 \Delta U |\Delta U| \frac{k\omega}{g} \frac{\rho_a}{\rho_w}.\]

The wind speed entering this term is evaluated at half wavelength above the surface:

\[U_{\lambda/2} = U_z + \frac{u_*}{\kappa} \log\left(\frac{\lambda/2}{z}\right),\]

with $\kappa=0.4$, hence $1/\kappa=2.5$ in the code. In coupled ESMF builds, stability-function corrections are added to this log-profile term.

The sheltering coefficient depends on sea state:

\[A_1 = \begin{cases} \text{variable}, & \Delta U > 0 \quad \text{wind sea},\\ 0.001, & 0 < U_{\lambda/2}\cos\theta < c+u\cos\phi+v\sin\phi \quad \text{swell with wind},\\ 0.1, & \cos\theta < 0 \quad \text{swell against wind}. \end{cases}\]

In the current source, the variable positive-input coefficient is computed by sheltering_coare35(wspd), then negative input is scaled by sin_diss1/sin_fac, and swell-overrunning-wind input is further scaled by sin_diss2/sin_diss1.

The sheltering_coare35 approximation is piecewise in wind speed $U$:

\[A_1(U) = \begin{cases} i_0 + m_1 U, & U \le x_1,\\ a + bU + cU^2, & x_1 < U \le x_2,\\ s_1\exp\left[-\dfrac{U-x_2}{d_c U}\right], & U > x_2, \end{cases}\]

with

\[\begin{aligned} x_1 &= 15,\quad x_2 = 33,\quad x_3 = 60,\\ y_1 &= 0.10,\quad y_2 = 0.09,\quad y_3 = 0.06,\\ i_0 &= 0.04,\quad c_r = 0.65,\quad d_c = 1.6,\\ m_1 &= \frac{y_1-i_0}{x_1},\quad m_2 = \frac{y_3-y_2}{x_3-x_2},\\ c &= c_r\frac{m_2-m_1}{x_2-x_1},\\ b &= m_1 - 2cx_1,\\ a &= y_1 - bx_1 - cx_1^2,\\ s_1 &= a + bx_2 + cx_2^2. \end{aligned}\]

Sea ice reduces wind input by multiplying ssin by $(1-f_{ice})$. For bins above the local diagnostic cutoff, negative wind input is clipped to zero.

Breaking Dissipation Sds

Breaking dissipation is nonlinear in the saturation spectrum

\[B(k,\phi) = k^4 E(k,\phi).\]

The physical dissipation source is

\[S_{ds}(k,\phi) = -A_2 \coth(0.2kd) \left(1 + A_3\overline{\chi^2}(k,\phi)\right)^2 B^n(k,\phi)\, \omega E(k,\phi).\]

The implementation stores the positive decay rate

\[S_{ds}^* = A_2 \coth(0.2kd) \left(1 + A_3\overline{\chi^2}\right)^2 B^n \omega.\]

Default coefficients are

\[A_2 = 42,\qquad A_3 = 360,\qquad n = 2.4.\]

The longer-wave mean-square-slope contribution is computed directionally:

\[\overline{\chi^2}_{o,p} = \sum_{q<o}\sum_{p'} E_{q,p'} \cos^2(\phi_{p'}-\phi_p) k_q^3\,dk_q\,d\phi.\]

Because $\cos^2 x = (1 + \cos 2x)/2$, the double sum reduces to three direction-integrated moments of the longer waves,

\[\overline{\chi^2}_{o,p} = M^{(0)}_o + \cos(2\phi_p)\,M^{(c)}_o + \sin(2\phi_p)\,M^{(s)}_o, \qquad M^{(0,c,s)}_o = \frac{d\phi}{2} \sum_{q<o} k_q^3\,dk_q \sum_{p'} E_{q,p'}\,\{1, \cos 2\phi_{p'}, \sin 2\phi_{p'}\},\]

which mss_moments evaluates in $O(N_f N_\phi)$ rather than $O(N_f N_\phi^2)$ operations. The power of the saturation spectrum is evaluated in double precision.

The code computes and stores the spilling part, before applying $\coth(0.2kd)$. This lets Snl use the spilling-only dissipation rate. The factor cothkd is applied where sds is used, to include the finite-depth/plunging-breaker enhancement in the decay rate used by the source integrator, the stress diagnostics, and the Sds output. Between output times, sds is evaluated only for the frequencies used by the integration, up to the $m+1$ bins above the prognostic range from which Snl moves energy into it (two at the default resolution).

Turbulent Dissipation Sdt

Ambient turbulence in the wave boundary layer attenuates waves:

\[S_{dt} = -A_4 u_{*w}kE(k,\phi).\]

The implementation stores the positive decay rate

\[S_{dt}^* = A_4 u_{*w} k.\]

The code approximates the water-side friction velocity from the air-side friction velocity and density ratio:

\[u_{*w} = u_*\sqrt{\frac{\rho_a}{\rho_w}}.\]

Default:

\[A_4 = 0.002.\]

Viscous Dissipation Sdv

Viscous dissipation converts wave energy to heat:

\[S_{dv} = -4\nu k^2 E(k,\phi).\]

The implementation stores

\[S_{dv}^* = 4\nu k^2,\]

where $\nu$ is nu_water. This term is usually important only for the shortest waves in clean water.

Bottom Friction and Percolation Sbf

The manual separates bottom friction and bottom percolation:

\[S_{bf} = -G_f\frac{k}{\sinh(2kd)}E(k,\phi),\] \[S_{bp} = -G_p\frac{k}{\cosh^2(kd)}E(k,\phi).\]

The implementation combines both into the positive decay-rate array sbf:

\[S_{bottom}^* = G_f\frac{k}{\sinh(2kd)} + G_p\frac{k}{\cosh^2(kd)}.\]

Typical coefficient ranges are:

\[G_f = 0.001\text{ to }0.01\,\mathrm{m\,s^{-1}}, \qquad G_p = 0.0006\text{ to }0.01\,\mathrm{m\,s^{-1}}.\]

The defaults are sbf_fac = 0.003 and sbp_fac = 0.003.

Nonlinear Downshifting Snl

UMWM represents nonlinear wave-wave transfer by a downshifting of the energy of spilling breakers: at $A_5$ times the spilling breaking rate, energy moves from each frequency $f$ to the lower frequency $f e^{-\Delta}$, by the fixed distance $\Delta$ in $\ln f$. With $\mathcal{E}(\ln f,\phi) = E\,k\,dk/d\ln f$, the variance per unit $\ln f$ and direction,

\[\left.\frac{\partial \mathcal{E}(\ln f,\phi)}{\partial t}\right|_{nl} = A_5 \left[ S_{sb}(\ln f+\Delta,\phi)\,\mathcal{E}(\ln f+\Delta,\phi) - S_{sb}(\ln f,\phi)\,\mathcal{E}(\ln f,\phi) \right],\]

where

\[S_{sb} = \frac{S_{ds}^*}{\coth(0.2kd)}\]

is the spilling-only breaking rate. The transfer conserves energy, moves it on average by exactly $\Delta$, and refers to neither the spectral resolution nor the wavenumber spacing of the bins: it is an energy flux through spectral space, the same in wavenumber and in frequency coordinates. snl_fac is $A_5$; the default is 5.0. The shift is the constant snl_shift of src/umwm_source_functions.f90,

\[\Delta = (b_1 + 2b_2)\,\Delta\ln f_0 = 0.155346, \qquad b_1 = e^{-16\Delta\ln f_0^2}, \qquad b_2 = e^{-64\Delta\ln f_0^2},\]

renormalized so $b_1+b_2=1$: the mean shift of the redistribution of Donelan et al. (2012) over the next two lower bins with the weights $b_1$ and $b_2$, at the default resolution, 37 bins from 0.0313 to 2 Hz, $\Delta\ln f_0 = \ln(2/0.0313)/36$. That redistribution moved energy by one or two bins, whatever their width, so that its mean shift depended on the resolution: 0.243 at $\Delta\ln f = 0.225$, 0.155 at 0.115, 0.118 at 0.083, and 0.057 at 0.038.

On bins of width $\Delta\ln f$, the shift is $\Delta/\Delta\ln f = m + \varphi$ bins, with $m$ whole bins (snl_bins) and the fraction $\varphi$ (snl_frac). The energy that leaves bin $j$ moves to bins $j-m$ and $j-m-1$ with the fractions $1-\varphi$ and $\varphi$: to the bins that bin $j$ overlaps when shifted by $\Delta$, a first-order conservative remapping. On bins wider than the shift ($m = 0$), this is first-order upstream finite-volume advection of energy at the speed $A_5\Delta S_{sb}$ in $\ln f$. On finer bins, energy moves over more than one bin per transfer, so the downshifting does not limit the time step more as the bins get narrower. For the variance spectrum,

\[S_{nl,o,p} = \beta_{1,o} S_{sb,o+m,p}E_{o+m,p} + \beta_{2,o} S_{sb,o+m+1,p}E_{o+m+1,p} - A_5\lambda_o S_{sb,o,p}E_{o,p},\] \[\beta_{1,o} = A_5 (1-\varphi) \frac{(k\,dk)_{o+m}}{(k\,dk)_o}, \qquad \beta_{2,o} = A_5 \varphi \frac{(k\,dk)_{o+m+1}}{(k\,dk)_o},\]

where $\lambda_o = (1-\varphi)[o>m] + \varphi[o>m+1]$ is the fraction of the energy of bin $o$ that moves: the energy that would move below the lowest bin stays where it is. The ratios of $k\,dk$ only convert between the energy of a bin, $E\,k\,dk\,d\phi$, and the spectrum, which makes the transfer conserve energy for any wavenumber spacing. Within the prognostic range the downshifting conserves energy: it changes the energy of the range only by what it moves in from the diagnostic bins $o_c+1$ to $o_c+m+1$, the downshifting through the cutoff frequency. The highest prognostic bin is therefore at most $N_f - \max(2, m+1)$, and the frequency range must be wider than the shift, $f_{\max}/f_{\min} > e^{\Delta} = 1.168$, which the model checks at startup.

At the default resolution, $m = 1$ and $\varphi = b_2$: the downshifting is the two-bin redistribution of Donelan et al. (2012), apart from the lowest two bins, which no longer lose energy below the spectrum.

In the dynamic physics time-step estimate, the breaking rate is reduced to $S_{ds}^*(1 - \beta_{1,o} - \beta_{2,o})$, which accounts for the rate at which the downshifting moves energy into a bin from bins with the same $E\,S_{sb}$ (see the time step below). As $\Delta$ is fixed, so is this estimate, apart from the interpolation of $k\,dk$ between bins.

Sea-Ice Attenuation Sice

The manual identifies sea-ice attenuation as an active source/sink term. The current source implementation follows the attenuation form cited in the source comments as Kohout et al. (2014).

If the sea-ice fraction is below fice_lth, sice = 0. If it is above fice_lth, the code computes local significant wave height

\[H_s = 4\sqrt{\sum_{o,p} E_{o,p} k_o\,dk_o\,d\phi}.\]

Then it defines an attenuation coefficient

\[\alpha = \begin{cases} C_1 H_s, & H_s < H_{th},\\ C_1 H_{th}, & H_s \ge H_{th}, \end{cases}\]

with

\[H_{th}=3\,\mathrm{m}, \qquad C_1=-5.35\times10^{-6}\,\mathrm{m^{-1}}.\]

The implemented rate is

\[S_{ice}^* = 2c_g\alpha.\]

Because $C_1<0$, this is normally a sink. When fice > fice_uth, propagation and refraction set the updated spectrum to zero in that cell. The default thresholds are fice_lth = 0.30 and fice_uth = 0.75.

Source Integration and Diagnostic Tail

The local source contribution is split from propagation and refraction:

\[\frac{\partial E}{\partial t} = \left(\frac{\partial E}{\partial t}\right)_s + \left(\frac{\partial E}{\partial t}\right)_a.\]

For source rates that multiply $E$,

\[\left(\frac{\partial E}{\partial t}\right)_s = \sum_i S_i^*E.\]

The exponential source solution is

\[E_s^{n+1} = E^n \exp\left(\sum_i S_i^*\Delta t\right).\]

The timestep is limited so the exponential factor remains bounded:

\[\exp\left(\sum_i S_i^*\Delta t\right) < r.\]

In code, explim is the logarithmic exponent limit and the diagnostic physics timestep is

\[\Delta t_{phys}(i) = \frac{\texttt{explim}} {\max_{o,p}|\Lambda_{o,p,i}|},\]

where the estimated local rate is

\[\Lambda = S_{in}^* - S_{ds}^*(1 - \beta_{1,o} - \beta_{2,o}) - S_{bottom}^* - S_{dt}^* - S_{dv}^* + S_{ice}^*,\]

with the factors $\beta$ of the downshifting Snl, which is integrated explicitly.

The actual source timestep is

\[dts = \min\left( \min_i \Delta t_{phys}(i), dtamin, dtg-sumt \right).\]

dtamin is the advective CFL limit computed at initialization:

\[dtamin = 0.98\,cfl_{lim} \frac{\min(\Delta x,\Delta y)}{\max(c_g)}.\]

The CFL factor is

\[cfl_{lim} = \begin{cases} \cos\left(\frac{\pi}{4}-\frac{\Delta\phi}{2}\right), & N_\phi \text{ divisible by } 8,\\ \frac{1}{\sqrt{2}}, & \text{otherwise}. \end{cases}\]

For the first non-restart step, dts is set to zero so only the diagnostic spectrum is initialized.

The prognostic/diagnostic cutoff frequency is

\[f_c = \texttt{fcut\_fac}\,\frac{g}{U_{10}},\]

by default with fcut_fac = 0.53, so that $f_c \approx 4f_{PM}$, and limited by fprog. Where waves are far younger than a fully developed sea, as in a laboratory tank, their spectral peak lies above this default cutoff, and a larger fcut_fac keeps them in the prognostic range. The local cutoff bin oc(i) marks the highest prognostic bin.

For bins above oc(i), UMWM assumes a wind-aligned equilibrium range. When the numerator is nonnegative, the diagnostic spectrum is assigned from the balance of wind input, turbulence, viscosity, sea ice, and breaking:

\[E_{eq} = k^{-4} \left[ \frac{ S_{in}^* - S_{dt}^* - S_{dv}^* + S_{ice}^* }{ 2\pi A_2 f \left(1+A_3\overline{\chi^2}\right)^2 \coth(0.2kd) } \right]^{1/n}.\]

After the source calculation, the implementation uses the time-split intermediate state

\[E^* = \frac{E^n + E_s^{n+1}}{2}\]

for propagation, refraction, and the frequency shift. Each computes its change from $E^*$, and the changes are added to $E_s^{n+1}$. Each is stable and keeps the spectrum positive on its own, propagation by the time step limit dtamin, and refraction and the frequency shift by their substeps; their sum is not guaranteed to keep it positive where the outflows of several of them from one bin are large. In coupled mode (ESMF), refraction and the frequency shift compute their changes from the propagated spectrum instead.

Geographic Propagation

The advection equation in Cartesian projection is

\[\frac{\partial A}{\partial t} = -\frac{\partial[(c_g\cos\phi+u)A]}{\partial x} -\frac{\partial[(c_g\sin\phi+v)A]}{\partial y},\]

with $A = E\,k/c_g$ in each cell. In terms of $E$, the flux into a cell from a neighbor carries the spectrum of the neighbor times the ratio of its $k/c_g$ to that of the cell; the flux out of the cell is unchanged.

Geographic propagation uses first-order upstream differencing. For one dimension,

\[\frac{\partial(\dot{x}A)}{\partial x} \approx \frac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x},\]

with

\[\Phi_{i+1/2} = \frac{\dot{x}_{i+1/2}+|\dot{x}_{i+1/2}|}{2}A_i + \frac{\dot{x}_{i+1/2}-|\dot{x}_{i+1/2}|}{2}A_{i+1},\] \[\Phi_{i-1/2} = \frac{\dot{x}_{i-1/2}+|\dot{x}_{i-1/2}|}{2}A_{i-1} + \frac{\dot{x}_{i-1/2}-|\dot{x}_{i-1/2}|}{2}A_i,\]

and

\[\dot{x}_{i\pm1/2} = \frac{\dot{x}_i+\dot{x}_{i\pm1}}{2}.\]

The implementation applies this form over east, west, north, and south faces using cell-face lengths and cell area. Propagation by currents is added only when the current fields are nonzero.

Open boundary conditions are applied next to land and at regional domain edges unless boundary input is supplied. Global domains use periodic east-west boundary conditions. Land cells take, at each frequency, the smallest phase and group speed over all sea cells.

Refraction

Refraction turns the wavenumber direction $\phi$ along the rays. For the dispersion relation $\omega = \sigma(k, d) + \mathbf{k}\cdot\mathbf{u}$, the ray equation $\dot{\mathbf{k}} = -\nabla\omega$ gives the rotation rate

\[\dot{\phi} = -\frac{c_g}{c}\frac{\partial c}{\partial n} -\hat{\mathbf{k}}\cdot\frac{\partial\mathbf{u}}{\partial n},\]

where $\hat{\mathbf{k}}=(\cos\phi,\sin\phi)$, $n$ is the coordinate normal to $\hat{\mathbf{k}}$, to its left, with $\partial/\partial n = -\sin\phi\,\partial/\partial x + \cos\phi\,\partial/\partial y$, and $c$ and $c_g$ are the intrinsic phase and group speeds at the fixed frequency of the bin. In components,

\[\dot{\phi} = \frac{c_g}{c}\left(\sin\phi\frac{\partial c}{\partial x} - \cos\phi\frac{\partial c}{\partial y}\right) + \sin^2\phi\frac{\partial v}{\partial x} - \cos^2\phi\frac{\partial u}{\partial y} + \sin\phi\cos\phi\left(\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}\right).\]

The depth term is Snell’s law, $d\phi/ds = -c^{-1}\,\partial c/\partial n$, along the ray $s$ that the wave energy travels at the group speed: the phase speed gradient sets the direction in which the waves turn, and the group speed how fast they cross it. Over straight and parallel depth contours it keeps $\sin\phi/c$ constant along the ray. The group speed at fixed frequency is not monotonic in depth, largest at $kd \approx 1.2$, so its gradient does not give the direction of refraction.

The current term is half the vorticity, $\frac{1}{2}(\partial v/\partial x - \partial u/\partial y)$, plus a part from the strain of the current that varies with direction. A current in solid-body rotation turns the waves of all directions at half its vorticity; a shear current $u(y)$ turns the waves along it at $-\partial u/\partial y$, and does not turn the waves across it. The ray path, in the direction of $\mathbf{c}_g + \mathbf{u}$, turns at the full vorticity to first order in $u/c_g$ (Kenyon, 1971), from $\dot\phi$ together with the advection of the waves by the current.

Positive $\dot{\phi}$ rotates wave energy counter-clockwise; negative $\dot{\phi}$ rotates it clockwise.

The gradients are differences over the sea neighbors of a cell: the phase speed gradient times $c_g/c$ (cgdlncp_dx, cgdlncp_dy) of each frequency at initialization, and the current velocity gradient (du_dx, du_dy, dv_dx, dv_dy) at the first time step and each time step with currents from file. They are centered between two sea neighbors, one-sided next to land, and zero between two land neighbors, so that the values of land cells do not enter: the current, set to zero on land, and the phase speed, the smallest of all sea cells. A uniform current over a flat bottom then neither refracts the waves nor shifts their frequency, also at the coast, while over a resolved slope next to the coast the waves turn toward the shore. The direction $\phi$ is taken in the frame of the grid, as in the propagation.

Refraction is discretized with the same upstream-flux form, replacing $x$ by $\phi$, with the rotation rate at the faces between directions the mean of the rates of the two directions. Directional space is periodic, so no directional boundary condition is required.

On its own, the scheme is stable and keeps the spectrum positive if the fraction of a bin that leaves it through both of its faces in a time step does not exceed 1:

\[\frac{\Delta t}{\Delta\phi} \left(\max(\dot{\phi}_{p+1/2}, 0) - \min(\dot{\phi}_{p-1/2}, 0)\right) \le 1.\]

Refraction splits the time step dts of each cell into as many equal substeps as this condition needs, so that every cell is refracted over the whole time step. The condition is evaluated with an upper bound from the amplitudes of the depth term, $a$, and of the strain term, $b$, which vary over direction as $\phi$ and $2\phi$, and from half the vorticity, $\zeta/2$: a bin loses at most $a + |\zeta|/2 + b$ through one face, and at most $\Delta\phi\,(a + 2b)$ through both faces where the rotation diverges, per unit time. One substep is usually enough.

Frequency Shift

With currents, the intrinsic frequency changes along the rays at

\[\dot{\varsigma} = \frac{d\ln\sigma}{dt} = \frac{1}{\sigma}\frac{\partial\sigma}{\partial d}\,\mathbf{u}\cdot\nabla d - \frac{c_g}{c}\,\hat{\mathbf{k}}\cdot(\hat{\mathbf{k}}\cdot\nabla)\mathbf{u},\]

from the ray equations for $\omega = \sigma(k, d) + \mathbf{k}\cdot\mathbf{u}$ with a steady depth. At a fixed frequency, $(\partial\sigma/\partial d)\nabla d = k\,(c_g/c)\nabla c$, so the first term is $\mathbf{u}\cdot(c_g/c)\nabla c\,/\,c$, from the same phase speed gradient as refraction. In components, $\hat{\mathbf{k}}\cdot(\hat{\mathbf{k}}\cdot\nabla)\mathbf{u} = \cos^2\phi\,\partial u/\partial x

The action $A$ is advected in $\ln\sigma$, whose bins are uniformly spaced by dlnf, with the same upstream-flux form as refraction, for the prognostic frequencies of each cell. Action leaves the lowest bin through its lower face, with none below it, and flows between the highest prognostic bin and the diagnostic tail above it. The time step of each cell is split into as many equal substeps as the upstream scheme needs to stay stable and positive, from an upper bound of $|\dot{\varsigma}|$ over frequency and direction: the strain rate $\hat{\mathbf{k}}\cdot(\hat{\mathbf{k}}\cdot\nabla)\mathbf{u}$ is half the divergence of the current plus a sinusoid in $2\phi$. The frequency shift is evaluated only when the current is nonzero.

Wind Stress and Drag

Wave energy and momentum are related by phase speed:

\[E = Mc.\]

The resolved wind form stress is computed from the wind-input source:

\[\tau_x = \rho_w g \int_{-\pi}^{\pi} \int_{k_{\min}}^{k_{\max}} \frac{S_{in}}{c} \cos\phi\,k\,dk\,d\phi,\] \[\tau_y = \rho_w g \int_{-\pi}^{\pi} \int_{k_{\min}}^{k_{\max}} \frac{S_{in}}{c} \sin\phi\,k\,dk\,d\phi.\]

The form drag coefficient is

\[C_{d,f} = \frac{\sqrt{\tau_x^2+\tau_y^2}}{\rho_a U_z^2}.\]

The smooth-flow skin drag coefficient is computed from

\[C_{d,s} = \frac{u_*^2}{U^2(z)} = \frac{\kappa^2}{\left[\log(z/z_0)\right]^2},\]

with molecular-sub-layer roughness

\[z_0 = 0.132\frac{\nu_a}{u_*}.\]

The implementation solves this relation by six fixed-point iterations. It uses wind relative to the moving surface:

\[\mathbf{U}_{rel} = \mathbf{U}_{wind} - \left(\mathbf{u}_{current}+\mathbf{u}_{Stokes}(z=0)\right).\]

The skin drag is reduced by wave sheltering:

\[C_{d,s}^{new} = \frac{C_{d,s}^{old}}{3} \left( 1 + 2\frac{C_{d,s}^{old}}{C_{d,s}^{old}+C_{d,f}} \right),\]

and capped at $10^{-2}$.

The total atmospheric stress and drag coefficient are

\[\boldsymbol{\tau} = \boldsymbol{\tau}_{form} + \boldsymbol{\tau}_{skin}, \qquad C_d = \frac{|\boldsymbol{\tau}|}{\rho_a U_z^2}.\]

The air-side friction velocity is updated from the total stress:

\[u_* = \sqrt{\frac{|\boldsymbol{\tau}|}{\rho_a}}.\]

High-Frequency Stress Tail

For accurate form stress, the manual recommends a high-frequency limit of at least 2 Hz. The implementation appends a wind-speed-dependent power-law stress tail from the highest resolved wavenumber to a capillary wavenumber:

\[p(U) = aU^2+bU+c,\]

with

\[a=0.000112,\qquad b=-0.01451,\qquad c=-1.0186.\]

The tail factor is

\[T = \frac{k_{cap}^{p+1}-k_{max}^{p+1}} {k_{max}^{p}(p+1)}, \qquad k_{cap}=10^3\,\mathrm{m^{-1}}.\]

This tail is applied in the wind direction for atmospheric stress and in the dissipative direction for ocean-top momentum flux.

Ocean Momentum Fluxes

Momentum fluxes into the ocean are defined positive downward. Ocean-top flux receives dissipation by breaking, turbulence, and viscosity, plus tail and skin stress:

\[\tau_x^{OT} = \rho_w g \int_{-\pi}^{\pi} \int_{k_{\min}}^{k_{\max}} \frac{-S_{ds}-S_{dt}-S_{dv}}{c} \cos\phi\,k\,dk\,d\phi + \tau_x(tail) + \tau_x(skin),\] \[\tau_y^{OT} = \rho_w g \int_{-\pi}^{\pi} \int_{k_{\min}}^{k_{\max}} \frac{-S_{ds}-S_{dt}-S_{dv}}{c} \sin\phi\,k\,dk\,d\phi + \tau_y(tail) + \tau_y(skin).\]

Ocean-bottom flux receives bottom friction and bottom percolation:

\[\tau_x^{OB} = \rho_w g \int_{-\pi}^{\pi} \int_{k_{\min}}^{k_{\max}} \frac{-S_{bf}-S_{bp}}{c} \cos\phi\,k\,dk\,d\phi,\] \[\tau_y^{OB} = \rho_w g \int_{-\pi}^{\pi} \int_{k_{\min}}^{k_{\max}} \frac{-S_{bf}-S_{bp}}{c} \sin\phi\,k\,dk\,d\phi.\]

The source additionally diagnoses taux_snl and tauy_snl, the momentum loss associated with nonlinear downshifting. Snl conserves energy in the implemented transfer but does not conserve momentum because phase speed changes between source and receiver bins.

Stokes Drift

When OUTPUT: stokes = .TRUE., the model reads user-provided depth levels from the STOKES namelist and computes wave-induced Stokes drift at those levels.

Depths in the namelist are positive downward. Internally they are stored as $z_l \le 0$, with $z=0$ at the surface and $z=-d$ at the bottom.

For each requested level, UMWM computes

\[\mathbf{u}_S(z_l) = \sum_{o,p} K_{o,l} E_{o,p} \begin{bmatrix} \cos\phi_p\\ \sin\phi_p \end{bmatrix},\]

where, for finite depth,

\[K_{o,l} = \omega_o k_o^2 \frac{\cosh\left[2k_o(z_l+d)\right]}{\sinh^2(k_od)} dk_o\,d\phi.\]

For large arguments that would overflow hyperbolic functions, the code uses the deep-water approximation

\[K_{o,l} = 2\omega_o k_o^2 \exp(2k_o z_l) dk_o\,d\phi.\]

If a requested level is below the local bottom, its Stokes contribution is set to zero.

The output d_stokes is the e-folding depth of the Stokes drift magnitude, found by linear interpolation to the depth where

\[|\mathbf{u}_S(z)| = \frac{|\mathbf{u}_S(0)|}{e}.\]

Integrated Diagnostics

The following diagnostics are computed from the full spectrum unless noted otherwise.

Significant wave height:

\[H_s = 4\sqrt{ \int\!\int E(k,\phi)k\,dk\,d\phi }.\]

Discrete implementation:

\[H_s = 4\sqrt{ \sum_{o,p}E_{o,p} k_o\,dk_o\,d\phi }.\]

Mean wave period:

\[\overline{T} = \sqrt{ \frac{ \int\!\int E(k,\phi)k\,dk\,d\phi }{ \int\!\int f^2E(k,\phi)k\,dk\,d\phi } }.\]

Mean wavelength:

\[\overline{L} = 2\pi \sqrt{ \frac{ \int\!\int E(k,\phi)k\,dk\,d\phi }{ \int\!\int k^2E(k,\phi)k\,dk\,d\phi } }.\]

Mean wave direction:

\[\overline{\phi} = \operatorname{atan2} \left( \sum_p M_p\sin\phi_p, \sum_p M_p\cos\phi_p \right),\]

where

\[M_p = \sum_o E_{o,p}k_o\,dk_o.\]

Mean-square slope:

\[mss = \sum_{o,p} E_{o,p}k_o^3\,dk_o\,d\phi.\]

The dominant period, wavelength, and direction are taken from the discrete spectral peak of

\[E_{o,p}k_o\,dk_o.\]

The wave momentum diagnostics are integrated over the prognostic range:

\[M_x = \rho_w g \sum_{o\le oc,p} E_{o,p} \frac{k_o\,dk_o}{c_o} \cos\phi_p\,d\phi,\] \[M_y = \rho_w g \sum_{o\le oc,p} E_{o,p} \frac{k_o\,dk_o}{c_o} \sin\phi_p\,d\phi.\]

The radiation-like momentum flux components are

\[C_{g}M_{xx} = \rho_w g \sum_{o\le oc,p} c_{g,o}E_{o,p} \frac{k_o\,dk_o}{c_o} \cos^2\phi_p\,d\phi,\] \[C_{g}M_{xy} = \rho_w g \sum_{o\le oc,p} c_{g,o}E_{o,p} \frac{k_o\,dk_o}{c_o} \cos\phi_p\sin\phi_p\,d\phi,\] \[C_{g}M_{yy} = \rho_w g \sum_{o\le oc,p} c_{g,o}E_{o,p} \frac{k_o\,dk_o}{c_o} \sin^2\phi_p\,d\phi.\]

Physics Namelist Controls

The main physics controls live in the PHYSICS namelist, with the following defaults.

Parameter Default Units Role
g 9.80665 m/s^2 Gravitational acceleration.
nu_air 1.56E-5 m^2/s Air kinematic viscosity for skin roughness.
nu_water 0.90E-6 m^2/s Water kinematic viscosity for Sdv.
sfct 0.07 N/m Surface tension in dispersion relation.
kappa 0.4 nondimensional Von Karman constant.
z 10. m Input wind height.
gustiness 0.0 nondimensional Random wind-component perturbation fraction; should be nonnegative and not larger than about 0.2.
dmin 10. m Minimum depth applied to ocean cells.
explim 0.9 nondimensional Maximum source exponent used to choose dts; 0.69 permits about 100 percent growth.
sin_fac 0.11 nondimensional Wind-sea growth factor.
sin_diss1 0.10 nondimensional Opposing-wind damping factor.
sin_diss2 0.001 nondimensional Swell-overrunning-wind damping factor.
sds_fac 42. nondimensional Breaking dissipation factor $A_2$.
sds_power 2.4 nondimensional Saturation-spectrum exponent $n$.
mss_fac 360 nondimensional Mean-square-slope enhancement factor $A_3$.
snl_fac 5.0 nondimensional Nonlinear downshifting factor $A_5$: the rate at which energy moves by snl_shift = 0.155346 in $\ln f$, relative to the spilling breaking rate.
sdt_fac 0.002 nondimensional Turbulent dissipation factor $A_4$.
sbf_fac 0.003 m/s Bottom friction coefficient $G_f$.
sbp_fac 0.003 m/s Bottom percolation coefficient $G_p$.
fcut_fac 0.53 nondimensional Factor of the prognostic cutoff frequency $f_c$ = fcut_fac$\,g/U_{10}$, about 4 times the Pierson-Moskowitz peak frequency by default.

Physics-related DOMAIN controls:

Parameter Role
om Number of frequency/wavenumber bins.
pm Number of direction bins. The manual requires divisibility by 4; the CFL optimization is best when divisible by 8.
fmin, fmax Lowest and highest frequency bins.
fprog Upper cap on the prognostic frequency range. The local cutoff is min(fcut_fac*g/U10, fprog).
dtg Global forcing/output step and maximum source-loop interval.

Sea-ice controls:

Parameter Default Role
fice0 0.0 Constant sea-ice fraction when sea ice is not read from input files.
fice_lth 0.30 Lower sea-ice fraction threshold for attenuation.
fice_uth 0.75 Upper sea-ice fraction threshold above which propagation/refraction zero the updated spectrum.

OUTPUT namelist controls:

Parameter Default Units Role
outgrid 1 hours Gridded output interval. A value of 0 disables gridded output.
outspec 0 hours Spectrum point output interval. A value of 0 disables spectrum output.
outrst 6 hours Restart output interval. A value of 0 disables restart output.
xpl 20 grid index Grid-cell x index used for stdout diagnostic reporting.
ypl 6 grid index Grid-cell y index used for stdout diagnostic reporting.
stokes .t. logical Enables Stokes drift calculation and output; when true, the STOKES namelist is read.

STOKES namelist controls:

Parameter Default Units Role
depths 0.1 0.2 ... 100 m Positive-down depth levels where u_stokes, v_stokes, and d_stokes diagnostics are computed.

GRID controls of a regular longitude-latitude grid, which the model makes when gridFromFile = .F. and dlon and dlat are positive, instead of the Cartesian grid of delx and dely:

Parameter Default Units Role
lon0, lat0 0 degrees Longitude and latitude of cell (1,1).
dlon, dlat 0 degrees Spacing in longitude and latitude; 0 for a Cartesian grid. A global grid (isGlobal) has mm*dlon = 360.

Wave partition and data source controls of the FORCING namelist:

Parameter Default Role
waves .F. Wave partitions from the input files: the initial spectrum, unless restarting, and the spectrum of the open boundaries of a regional domain.
source 'none' Data source of the input files, such as 'era5': the missing input files are prepared from it before the run. 'none' reads input files provided in the input directory.

Wave Partitions: Initial State and Open Boundaries

Reanalyses and forecasts rarely archive the wave spectrum; they describe the wave field by integrated parameters of its partitions, the wind sea and a few swell systems. With waves = .T. in the FORCING namelist, the model reads the partitions from the input files, and rebuilds its spectrum from them: at the start time for the initial state, instead of starting from rest, and at every forcing time at the open boundaries of a regional domain.

Partitions in the input files

The partitions are variables of the input files umwmin_<time>.nc (or umwmin.nc with steady), on the grid like the other forcing fields, one set per partition $p$ = 0 to 9 (by convention, 0 is the wind sea):

Variable Units Meaning
phs<p> m Significant wave height $H_p$; 0 where the partition is absent, and negative (-1) where there are no wave data, such as on land or sea ice in the dataset.
ptm<p> s Mean period $T_{m-1,0} = m_{-1}/m_0$, the energy period; or, instead,
ptp<p> s the peak period $T_p$.
pdir<p> rad Mean direction of propagation $\bar\phi_p$, counterclockwise from east, in the direction convention of the model.
pspr<p> rad Directional spread $\sigma_p = \sqrt{2(1-M_1)}$ (Kuik et al., 1988); optional, 0.5 rad by default.
pgam<p> 1 Peak enhancement $\gamma_p$ of the JONSWAP spectrum; optional, 3.3 by default.

The values must be finite, as the model is compiled with fast math, which may not tell a NaN; the Python package writes missing values as -1 for phs<p> and 0 for the others.

Parametric spectrum of a partition

The spectrum of a partition is a JONSWAP spectrum in frequency times a $\cos^{2s}$ distribution in direction (Longuet-Higgins et al., 1963):

\[S_p(f,\phi) \propto f^{-5}\exp\left[-\frac{5}{4}\left(\frac{f_p}{f}\right)^4\right] \gamma_p^{\exp\left[-(f-f_p)^2/(2\sigma_J^2 f_p^2)\right]} \cos^{2s}\left(\frac{\phi-\bar\phi_p}{2}\right),\]

with $\sigma_J$ = 0.07 below and 0.09 above the peak frequency $f_p$. The first circular moment of the distribution in direction is $M_1 = s/(s+1)$, so that the spread $\sigma_p$ gives

\[s = \frac{2}{\sigma_p^2} - 1,\]

from 0 (isotropic, $\sigma_p = \sqrt 2$) to 20000 ($\sigma_p$ = 0.01 rad).

The spectrum is built on the discrete frequencies $f_o$ of the model: with a mean period, the peak frequency is found by fixed-point iteration, $f_p \leftarrow f_p\,T(f_p)/T_{m-1,0}$, where $T(f_p)$ is the mean period $m_{-1}/m_0$ of the discrete spectrum, from $f_p = 0.9/T_{m-1,0}$, within $[f_1/2, f_{om}]$, so that the discrete spectrum has the given mean period whatever the frequency resolution. The variance of the partition is then distributed over the bins as

\[E_{o,q}\,k_o\,dk_o\,d\phi = \frac{H_p^2}{16}\, \frac{f_o S_p(f_o)}{\sum_{o'} f_{o'} S_p(f_{o'})}\, \frac{D(\phi_q)}{\sum_{q'} D(\phi_{q'})},\]

since the frequency bins have the width $f_o\,\Delta\ln f$, in the weights $k\,dk\,d\phi$ of the significant wave height of the model, with the wavenumbers of the depth of the cell. The discrete spectrum thus has the significant wave height $H_p$ exactly, in deep and shallow water, its mean period, and its mean direction and spread to the resolution of the directions. The partitions add up: the spectrum of a cell is the sum of those of its partitions, and is tiny where it has none.

src/umwm_partitions.f90 implements the partitions and their spectra, and tests/test_partitions.f90 tests the significant wave height, mean period, mean direction, and spread of the discrete spectra.

Initial state

Unless restarting, the model starts from the spectrum of the partitions of the input file of the start time, which must have them. The first time step is then, as from rest, of zero length: it sets the diagnostic range of the spectrum to the equilibrium of the model, and the stress and friction velocity of the initial spectrum, before the output of the start time. A restart reads the spectrum of the restart file instead, and the partitions serve the open boundaries.

The partitions of a dataset describe the wave field of its own wave model, and the model adjusts it to its own physics in the first hours. From ERA5, the global wave energy decreases by 7% in the first hour and 12% in six hours, and then holds: the model damps the swell that travels across or against the local wind (sin_diss1) much faster than the wave model of ERA5, which removes the components of the swell more than 90 degrees from the wind in a few hours, most in the trade winds that the swell of the Southern Ocean crosses.

Open boundaries

The edges of a domain are land. In a regional domain (not isGlobal), the sea cells next to its edges where the input file of the start time has wave data (phs<p> not negative) are open boundaries: their spectrum is that of the partitions of the input files, which they propagate into the domain. At each forcing time, the model rebuilds the spectrum of the boundary cells from the partitions of the new input file, which must have them, and after each source time step it sets their spectrum to

\[E = (1-\alpha)\,E_b + \alpha\,E_f,\qquad \alpha = \frac{t - t_b}{t_f - t_b},\]

between the spectra $E_b$ and $E_f$ of the backward and forward forcing times, as the forcing fields are interpolated. The source functions of the boundary cells are not integrated, since their spectrum is prescribed, and they do not limit the source time step: the spectrum of a dataset is out of equilibrium with the source functions of the model, and its fast relaxation would otherwise hold the time step of the whole domain to a few seconds. In the output, the wind input and nonlinear transfer of the boundary cells are therefore zero, and their spectrum keeps the parametric shape of the partitions up to the highest frequency, rather than the equilibrium range that the source functions set in the other cells. A global domain has no open boundaries. In ESMF coupled mode, the input files are not read after the start time, and the boundaries keep the spectrum of the start time.

Data Sources of the Input Files

The model reads its forcing, grid, and wave partitions from NetCDF files in its input directory. A data source prepares these files from a reanalysis or a forecast, on the grid of the model, so that a run from such a dataset needs no other program: with source = 'era5' in the FORCING namelist, the input files of the configuration that are missing are prepared before the run, by

They are prepared in the input directory of the run, CONFIG/input by default.

The namelist says which files the model reads, as without a source: the grid of lon0, lat0, dlon and dlat in GRID (or a grid file with gridFromFile), the bathymetry with topoFromFile, and the forcing fields with winds, air_density, seaice and waves in FORCING. The files that exist are kept, such as a bathymetry from a finer dataset, unless they are for another grid, such as those of another case prepared in the same input directory: fields of another size than mm by nm, or, with lon and lat, first and last cells at other coordinates than those of the grid (within 0.001 degrees). The files for another grid that the preparation wrote (their history attribute starts with prepared by umwm) are prepared again; the others stop the preparation with an error, rather than being replaced. The model also stops with an error at an input file whose fields are not of mm by nm cells, whatever its source, rather than read part of it.

ERA5

The era5 source reads ERA5 (Hersbach et al., 2020), the hourly reanalysis of ECMWF from 1940, from the analysis-ready, cloud-optimized Zarr store of Google Cloud Public Datasets (Carver and Merose, 2023), on a regular 0.25 degree grid, anonymously over HTTPS, where each variable is a chunk of the whole globe per hour (UMWM_ERA5_URL sets another store, such as a local copy). It provides:

Input ERA5 variables
Bathymetry model_bathymetry of the wave model, sea where at least half of the weight of the bilinear interpolation is on its sea
winds 10m_u_component_of_wind, 10m_v_component_of_wind
air_density air_density_over_the_oceans
seaice sea_ice_cover
waves, partition 0 wind sea: significant_height_of_wind_waves, mean_period_of_wind_waves, mean_direction_of_wind_waves, wave_spectral_directional_width_for_wind_waves
waves, partitions 1-3 first, second and third swell: significant_wave_height_of_*_swell_partition, mean_wave_period_of_*_swell_partition, mean_wave_direction_of_*_swell_partition, and wave_spectral_directional_width_for_swell
waves, partition 4 the residual of the whole wave field: significant_height_of_combined_wind_waves_and_swell, mean_wave_period, mean_wave_direction, wave_spectral_directional_width

The mean periods of ERA5 are $T_{m-1,0}$, its directions are where the waves come from, clockwise from north, and its directional spreads are $\sqrt{2(1-M_1)}$ at the spectral peak (Bidlot, 2020). ERA5 does not provide currents or water density, which keep the constant values of FORCING_CONSTANT.

Preparation

The preparation, common to all sources (python/umwm/sources/prepare.py):

  1. reads the fields of the source at the times it needs, a block of times at a time, restricted to the latitudes of the grid of the model;
  2. interpolates them bilinearly to the cells of the model, periodic in longitude, leaving out the corners without data, and in time, linearly, between the times of the source; the wave partitions are those of the nearest time;
  3. converts the wave partitions: directions to the convention of the model, and, where the partitions have less energy than the whole wave field (they are interpolated separately in the dataset), adds the difference as a residual partition with the period, direction and spread of the whole wave field, or scales them down where they have more, so that the model starts from the significant wave height of the dataset; the period, direction and spread of a partition are interpolated weighted by its energy at the corners;
  4. writes umwm.gridtopo and the forcing files, with the wave partitions at the start time, and at every time if the domain is regional with sea cells next to its edges, its open boundaries.

Adding a data source

A source is a class of python/umwm/sources that reads its dataset in the standard names and conventions of python/umwm/sources/base.py (u10, v10, rhoa, fice, depth, and the partitions hs<p>, tm<p> or tp<p>, dir<p> and spr<p>, with hs, tm, dir and spr of the whole wave field), declared in SOURCES of python/umwm/sources/__init__.py with the inputs it provides and its dependencies, as a pip extra of the package; the preparation, the namelist, and the command-line tools then serve it as they serve ERA5. For example, an NCEP GFS source would read the 10 m wind, the air density from the pressure, temperature and humidity, and the sea ice of the GFS atmosphere, and the wave partitions of its WAVEWATCH III wave model (GFS-Wave: the wind sea WVHGT, WVPER, WVDIR and the swell SWELL, SWPER, SWDIR of three partitions), from GRIB2 files with byte ranges of their index files, from NOMADS or from the NOAA Open Data Dissemination buckets, with the peak or mean period of WAVEWATCH III (tp<p> or tm<p>), and the default spread where it has none; the forecast cycle is that of the start time.

Output Fields Relevant to Physics

Gridded output includes forcing, integrated wave quantities, stress, and diagnostics:

Field Meaning
wspd, wdir Wind speed and mathematical wind direction.
uc, vc Ocean current components.
rhoa, rhow Air and water density.
fice Sea-ice fraction.
taux_form, tauy_form Wind form drag components.
taux_skin, tauy_skin Skin drag components.
taux_ocn, tauy_ocn Momentum flux into ocean top.
taux_bot, tauy_bot Momentum flux into ocean bottom.
taux_snl, tauy_snl Momentum loss due to nonlinear downshifting.
cd Total atmospheric drag coefficient.
ust Air-side friction velocity.
shelt Variable sheltering coefficient.
swh, mwp, mwl, mwd Significant height, mean period, mean wavelength, mean direction.
mss Mean-square slope.
dwp, dwl, dwd Dominant period, wavelength, and direction from the discrete spectral peak.
momx, momy Wave momentum components.
cgmxx, cgmxy, cgmyy Group-velocity-weighted momentum flux components.
epsx_atm, epsy_atm Wave energy growth flux from air.
epsx_ocn, epsy_ocn Wave energy dissipation flux to ocean.
physics_time_step Local source-limited physics timestep diagnostic.
u_stokes, v_stokes, d_stokes Stokes drift components and e-folding depth, when enabled.

Spectrum point output includes:

Field Meaning
F Surface elevation variance spectrum.
Sin Wind input/export implementation rate.
Sds Wave breaking implementation decay rate.
Sdt Turbulent implementation decay rate.
Sdv Viscous implementation decay rate.
Sbf Combined bottom friction and bottom percolation implementation decay rate.
Snl Nonlinear wave-wave interaction absolute source/sink tendency.

The current code writes the implementation arrays directly in spectrum point output. Thus Sin, Sds, Sdt, Sdv, and Sbf are rates that must be multiplied by F to obtain variance-spectrum tendencies, while Snl is already an absolute tendency.

Tests

The focused unit tests live in tests/ and are built by tests/Makefile. They can be run either from the top level

make test

or directly from the test directory

make --directory=tests test

The current focused tests are:

Test Coverage
test_dispersion Wavenumber/frequency round-trip accuracy and positive finite group speed for the capillary-gravity dispersion module.
test_spectrum Spectrum constructor metadata, endpoint frequencies, legacy log-space frequency equivalence, direction-grid equivalence, dlnf, and dth.
test_source_functions Mean-squared slope moments against the direct sum over directions, the vectorized power function against the intrinsic, and the wind input against a step-by-step evaluation.
test_advection Current-induced ray rotation against the ray equation and in a solid-body eddy, depth-induced ray rotation against Snell’s law along a ray over a slope, the conservation, positivity, and rotation of the spectrum by refraction in several substeps, the conservation of action and the shift of the spectrum by the frequency shift, steady shoaling against linear theory, and the steady variance of waves against a strengthening current against ray theory.
test_domain One-hour integration of the idealized hurricane case by domain_type with gridded and spectrum output and Stokes drift, compared with the mean and maximum of output fields, ocean stress, and Stokes drift from UMWM before the domain refactor without the refraction by its spurious current gradients against land; the uniform current over the flat bottom neither refracts the waves nor shifts their frequency.

The Python package has its own test suite in python/tests, run with python -m pytest python/tests, which includes the same idealized hurricane case driven from Python and compared with the same reference values, and python/tests/mpi_smoke.py for a library built with MPI.

make clean at the top level recurses into tests/ and removes test-local binaries, object files, and module files.

References

The physics basis and manual equations are from UMWM manual version 2, sections 1 through 4 and 9. The cited references are:

This site is open source. Improve this page.