A Python-based implementation of the Model for Prediction Across Scales (MPAS) atmospheric dynamical core with tangent linear model (TLM) and adjoint model capabilities.
pyMPAS-Adjoint provides a Python interface to MPAS atmosphere dynamics with support for forward model simulation, tangent linear model (TLM), and adjoint model (ADJ) integration. This implementation uses f2py to wrap Fortran code, enabling efficient numerical computation while maintaining Python's ease of use for setup, analysis, and visualization.
The system implements the non-hydrostatic dynamical core on unstructured Voronoi meshes, supporting global atmospheric simulations with advanced numerical schemes including Runge-Kutta time integration (SRK3) and various horizontal mixing parameterizations.
- Forward Model: Full non-hydrostatic atmospheric dynamical core
- Tangent Linear Model (TLM): Linearized dynamics for sensitivity analysis
- Adjoint Model (ADJ): Backward integration for gradient computation and data assimilation
- Flexible Grid Support: Unstructured Voronoi meshes with variable resolution
- Advanced Numerics:
- 3rd-order Runge-Kutta time integration (SRK3)
- Multiple horizontal mixing schemes (2D Smagorinsky, fixed diffusion)
- Split-explicit time stepping with configurable sub-steps
- Scalar advection with monotonicity and positive-definite options
- Python-Fortran Integration: High-performance Fortran kernels accessible via Python
- NetCDF/HDF5 I/O: Standard atmospheric data formats for input/output
pyMPAS-Adjoint/
├── src/
│ ├── mpas_module.f90 # Main MPAS Fortran module with state variables
│ ├── mpas_atm_dyn.f90 # Forward model dynamics subroutines
│ ├── mpas_atm_dyn_tlm.f90 # Tangent linear model dynamics
│ ├── mpas_atm_dyn_adj.f90 # Adjoint model dynamics
│ ├── mpas_atm_core_subs.f90 # Core atmospheric physics subroutines
│ ├── mpas_operators.f90 # Mesh operators and interpolation
│ ├── mpas_driver.py # Python driver routines for model I/O
│ ├── mpas_namelist.py # Namelist parsing utilities
│ └── Makefile # Build configuration for f2py
├── run-fwd.py # Forward model execution script
├── run-tlm.py # Tangent linear model execution script
├── run-adj.py # Adjoint model execution script
├── namelist.atmosphere # Configuration file for model parameters
└── README.md # This file
- Python: 3.x
- Python Libraries:
- numpy
- xarray
- netCDF4
- h5py
- f90nml (for namelist parsing)
- pandas (for time utilities)
- Fortran Compiler: gfortran or Intel Fortran
- f2py: Included with numpy (>= 1.x)
Navigate to the src/ directory and compile the Fortran code:
cd src/
makeThis creates module_mpas.so, a Python-importable shared library containing the MPAS dynamical core.
Note: The Makefile uses gfortran by default with -fdefault-real-8 flag for double precision. For Intel Fortran, modify the FLAG variable in the Makefile to -real-size 64.
Model configuration is controlled through namelist.atmosphere. Key parameters include:
config_dt = 1800.0 ! Time step (seconds)
config_start_time = '2020-02-02_00:00:00'
config_run_duration = '00_02:00:00' ! Format: DD_HH:MM:SS
config_time_integration_order = 2 ! Time integration order
config_number_of_sub_steps = 6 ! Sub-steps per time stepconfig_split_dynamics_transport = false
config_horiz_mixing = '2d_fixed' ! Options: '2d_smagorinsky', '2d_fixed'
config_len_disp = 480000. ! Diffusion length scale (m)
config_scalar_advection = true
config_positive_definite = true
config_monotonic = trueinput = 'x1.10242.init.nc' ! Input file (initial conditions)
output = 'x1.10242.fwdState.nc' ! Output file
output_interval = '01:00:00' ! Output frequencyThe model requires an initialization file (NetCDF format) containing:
- State variables: u (horizontal velocity), w (vertical velocity), theta (potential temperature), rho (density), scalars (tracers)
- Mesh information: cell/edge/vertex coordinates, connectivity arrays, mesh metrics
- Base state: theta_base, rho_base for perturbation calculations
- Grid: zgrid (vertical coordinate), various interpolation coefficients
python run-fwd.pyWhat it does:
- Reads initial conditions from the file specified in
namelist.atmosphere - Integrates the non-hydrostatic equations forward in time using SRK3
- Outputs state variables (u, w, rho, theta, pressure, scalars) at specified intervals
- Writes results to NetCDF file with xarray
Output: NetCDF file containing time series of 3D atmospheric fields on the unstructured mesh.
python run-tlm.pyWhat it does:
- Initializes the forward model state
- Allocates TLM variables with
mpas.var_allocation_tlm() - Applies small perturbations to initial state (configurable, typically 10^-5)
- Integrates both forward and TLM equations simultaneously
- The TLM evolution is: δX(t+dt) = M'[X(t)] · δX(t) where M' is the tangent linear operator
Use cases:
- Sensitivity analysis
- Singular vector calculation
- Forecast error growth studies
- Verification of adjoint model correctness
python run-adj.pyWhat it does:
- Loads forward model trajectory from saved state file (
FWD-States.h5) - Initializes adjoint forcing (typically from final TLM state for verification)
- Integrates adjoint equations backward in time
- At each time step, uses saved forward state to compute adjoint sensitivities
- Performs adjoint identity check: ⟨M'δX, δX⟩ = ⟨δX₀, M*δX⟩
Use cases:
- Gradient-based optimization
- Adjoint sensitivity analysis
- Data assimilation (4D-Var)
- Observation impact studies
Adjoint Verification: The script computes LHS and RHS of the adjoint identity and reports relative differences to verify correctness.
Core driver functions:
read_configs(nml): Parse and set model configuration from namelistread_dims(nml): Read mesh dimensions from NetCDF inputread_vars(nml): Load state variables and mesh dataclock_namelist(nml): Parse time settings and compute simulation duration
Namelist parsing with f90nml:
- Reads Fortran namelist format
- Provides defaults for unspecified parameters
- Handles time format conversion (HH:MM:SS to seconds)
Main module containing:
- State variables: u, w, rho, theta, scalars (forward)
- TLM variables: u_tl, w_tl, rho_tl, theta_tl, scalars_tl
- ADJ variables: u_ad, w_ad, rho_ad, theta_ad, scalars_ad
- Mesh arrays: cell/edge/vertex coordinates, connectivity, metrics
- Configuration parameters mapped from namelist
- Allocation routines:
var_allocation(),var_allocation_tlm(),var_allocation_adj() - Time-stepping:
atm_srk3(),atm_srk3_tlm(),atm_srk3_adj()
Forward model dynamics:
atm_compute_dyn_tend(): Computes tendencies for momentum (u), vertical velocity (w), density (rho), and potential temperature (theta)- Includes advection, pressure gradient force, Coriolis, mixing
- Implements horizontal mixing schemes (Smagorinsky, fixed diffusion)
- Scalar transport with various advection orders
Tangent linear model:
- Linearized versions of all forward model operators
atm_compute_dyn_tend_tlm(): TLM tendency calculations- Preserves structure of forward model for verification
- Used in sensitivity studies and Hessian calculations
Adjoint model:
- Adjoint of all forward/TLM operators
atm_compute_dyn_tend_adj(): Adjoint tendency calculations- Time-reversed integration (from final time to initial time)
- Requires stored forward trajectory for nonlinear terms
Geometric operators on unstructured mesh:
- Divergence and curl operators
- Edge-to-cell and cell-to-edge interpolation
- Reconstruction operators for velocity components
- Vertical interpolation
Time Integration: 3rd-order Strong Stability Preserving Runge-Kutta (SRK3)
- Explicit RK3 scheme for non-hydrostatic equations
- Configurable sub-stepping for acoustic/gravity wave terms
- Split-explicit option for dynamics and transport
Spatial Discretization:
- C-grid staggering on Voronoi mesh
- Normal velocity (u) at cell edges
- Vertical velocity (w), scalars at cell centers
- Vorticity at cell vertices
- 3rd-order advection schemes with flux limiters
Mixing Parameterization:
- 2D Smagorinsky turbulence closure
- Configurable horizontal/vertical eddy viscosity
- Monotonicity preservation for scalars
The adjoint code was developed using algorithmic differentiation principles:
- Forward code structured for AD compatibility
- TLM generated by linearizing forward operators
- Adjoint generated by transposing TLM operators and reversing time loop
- Careful treatment of in-place operations and overwriting
Verification Strategy:
- Tangent Linear Test: Compare TLM with finite differences
- Adjoint Identity Test: Verify ⟨M'δX, δX⟩ = ⟨δX₀, M*δX⟩ to machine precision
- Gradient Check: Validate cost function gradients for optimization problems
MPAS uses Centroidal Voronoi Tessellation (CVT):
- Cells: Polygonal cells (primarily hexagons) where scalars are defined
- Edges: Connect cell centers, where normal velocity is defined
- Vertices: Cell corners, where vorticity is computed
- Dual Grid: Triangle mesh formed by connecting cell centers
This structure provides natural conservation properties and eliminates pole singularities in global models.
- Initial Conditions: NetCDF following MPAS conventions
- State:
u,w,theta,rho,qv, etc. - Mesh:
latCell,lonCell,xCell, connectivity arrays - Metrics:
areaCell,dvEdge,dcEdge, interpolation weights
- State:
- Forward Model: NetCDF (xarray) with dimensions
(nTime, nLev, nCell/nEdge) - TLM/Adjoint: HDF5 for efficient state storage and checkpointing
- Forward trajectory: saved at each time step for adjoint
- TLM perturbations: saved for adjoint verification
# Edit namelist.atmosphere to set duration, time step, initial conditions
python run-fwd.py
# Output: x1.10242.fwdState.nc containing full 4D fields# Run TLM with initial perturbation
# Modify run-tlm.py to set perturbation pattern
python run-tlm.py
# Analyze growth of perturbations over time# First, run forward model and save trajectory
# Then initialize adjoint forcing based on cost function
python run-adj.py
# Obtain gradients with respect to initial conditions-
Skamarock, W. C., et al. (2012): "A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and C-Grid Staggering", Monthly Weather Review, 140(9), 3090-3105.
-
Errico, R. M. (1997): "What is an adjoint model?", Bulletin of the American Meteorological Society, 78(11), 2577-2591.
-
Thuburn, J., et al. (2009): "Numerical representation of geostrophic modes on arbitrarily structured C-grids", Journal of Computational Physics, 228(22), 8321-8335.
-
Tian, X., and Zou, X. (2020): "Development of the Tangent Linear and Adjoint Models of the MPAS-Atmosphere Dynamic Core and Applications in Adjoint Relative Sensitivity Studies", Tellus A: Dynamic Meteorology and Oceanography, 72(1), 1-20, doi:10.1080/16000870.2020.1814602.
-
Tian, X., and Zou, X. (2021): "Validation of a Prototype Global 4D-Var Data Assimilation System for the MPAS-Atmosphere Model", Monthly Weather Review, 149(8), 2803-2817, doi:10.1175/MWR-D-20-0408.1.
-
MPAS Documentation: https://mpas-dev.github.io/
MIT License - see LICENSE file for details.
Xiaoxu Tian (2024)
This implementation builds upon the MPAS-Atmosphere dynamical core developed by the MPAS development team at NCAR and collaborating institutions.
- For large-scale simulations, consider running the forward model separately and saving checkpoints for adjoint integration
- Memory requirements scale with mesh resolution and number of vertical levels
- The adjoint model requires stored forward states, increasing storage demands
- Current implementation focuses on dry dynamics; moist physics adjoint is under development (commented out in source)
Build Issues:
- Ensure f2py version matches numpy version
- Check Fortran compiler flags for precision (real8 vs real16)
- Verify all Fortran source files are in
src/directory
Runtime Issues:
- Verify input file path in namelist.atmosphere
- Check mesh dimensions match between code and input file
- Ensure sufficient memory for mesh size and time integration
- For adjoint runs, confirm forward state files exist and are accessible
Numerical Issues:
- Reduce time step if simulation becomes unstable
- Increase horizontal mixing for better stability
- Check CFL condition: dt < dx/(max wind speed)
- Enable monotonicity and positive-definite options for tracers