API Reference

prepare_input_data(lat_min, lat_max, lon_min, lon_max, year, month, days;
                   download_cams=false,
                   data_dir=joinpath(@__DIR__, "..", "data")) ->
    (era5_levels_file, era5_surface_file, era5_precip_file, cams_levels_file)

Prepare local meteorology inputs for a small regional TARSA run.

This helper is the recommended setup entry point for tutorial-style workflows. It reuses existing ERA5 files from data_dir when present and otherwise downloads them for the requested domain and date range.

CAMS download/regridding through this helper is currently disabled because the old xESMF backend has been removed. Regrid CAMS fields externally and pass the resulting file to the workflow that needs it.

Arguments

• lat_min, lat_max: Latitude bounds in degrees.
• lon_min, lon_max: Longitude bounds in degrees.
• year, month, days: Requested calendar period.

Keywords

• download_cams: Must be false. Setting it to true throws an ArgumentError because built-in CAMS regridding is no longer available.
• data_dir: Directory where downloaded/reused files are stored.

Returns

• era5_levels_file: ERA5 pressure-level NetCDF path.
• era5_surface_file: ERA5 single-level NetCDF path.
• era5_precip_file: ERA5 precipitation NetCDF path.
• cams_levels_file: Currently always ""; CAMS fields must be regridded externally.

Notes

• This function is aimed at small, example-scale domains rather than large production data preparation workflows.
• Download failures raise a clear error; there is no synthetic-data fallback on this path.

load_input_data(levels_file::String, surface_file::String; verbose=true) -> InputFields

Load and prepare meteorological data with structured output.

This function wraps read_input_data() and build_dx_dy_mats_from_latlon() to provide a cleaner interface with named fields instead of an 11-tuple.

Arguments

• levels_file: Path to ERA5 pressure levels NetCDF
• surface_file: Path to ERA5 single-level NetCDF
• verbose: Print progress information

Returns

InputFields struct with all meteorological data.

Example

input = load_input_data("era5_levels.nc", "era5_surface.nc")
@info "Grid size: $(input.Nx) x $(input.Ny) x $(input.Nz)"
@info "Mean dx: $(mean(input.dx)) m"

save_results(output_file::String, input::InputFields,
             concentration::Array{Float32,4}, emissions::Array{Float32,4};
             variables=["mixing_ratio", "concentration"],
             template_file::Union{Nothing,String}=nothing,
             verbose=true)

Save simulation results to NetCDF.

Arguments

• output_file: Path to output NetCDF file
• input: InputFields struct (provides grid info)
• concentration: 4D concentration array (kg/kg) from run!()
• emissions: 4D emission array (kg/m³/s)
• variables: List of variables to save
  • "mixing_ratio": kg/kg (always saved)
  • "concentration": kg/m³
  • "mass_loading": µg/m³
• template_file: NetCDF to copy dimensions from (default: use input file)
• verbose: Print progress
• save_emissions: Save the emissions field (default: true)
• save_density: Save air density as rho (default: false)
• *_name: Optional custom variable names for multi-species outputs

Example

sim = build_simulation(input, emissions)
concentration = run!(sim)

save_results("output.nc", input, concentration, emissions;
    variables=["mixing_ratio", "mass_loading"]
)

InputFields

Structured container for meteorological input data. Replaces the 11-tuple return from read_input_data().

Fields

• u, v, w: Wind components (m/s) - each (Nx, Ny, Nz, Nt)
• temperature: Temperature (K) - (Nx, Ny, Nz, Nt)
• humidity: Specific humidity (kg/kg) - (Nx, Ny, Nz, Nt)
• geopotential: Geopotential (m²/s²) - (Nx, Ny, Nz, Nt)
• density: Air density (kg/m³) - (Nx, Ny, Nz, Nt)
• pbl_height: Planetary boundary layer height (m) - (Nx, Ny, Nt)
• X_proj, Y_proj: Projected coordinates (m) - each (Nx, Ny)
• longitude, latitude: Geographic coordinates (degrees) - vectors
• dx, dy: Grid spacing (m) - each (Nx, Ny)
• Nx, Ny, Nz, Nt: Grid dimensions

Example

input = load_input_data("era5_levels.nc", "era5_surface.nc")
# Access fields: input.u, input.longitude, input.dx, etc.

Simulation

High-level simulation object that encapsulates all inputs and configuration.

Example

input = load_input_data(levels_file, surface_file)
emissions = build_point_source_emissions(input, lon, lat; rate_kg_s=1.0)

sim = build_simulation(
    input,
    emissions;
    solver=ExplicitSolver(cfl=0.8, scheme=:dst_koren_rk3),
    lateral_bc=:neumann0,
    bottom_bc=:neumann0,
)

result = run!(sim)

Grid(dx, dy, z_center)

Numerical grid description for provider-neutral model input. dx and dy contain strictly positive horizontal cell sizes in metres with shape (Nx, Ny). z_center contains cell-centre altitude in metres AGL with shape (Nx, Ny, Nz, Nt).

Meteorology(u, v, w, temperature, humidity, density, pbl_height)

Provider-neutral meteorological state. All four-dimensional fields use shape (Nx, Ny, Nz, Nt). pbl_height is the planetary boundary layer height in metres AGL with shape (Nx, Ny, Nt).

Geography(longitude, latitude, x_proj, y_proj, crs)

Coordinate metadata kept outside numerical kernels. Longitude and latitude may be one-dimensional coordinates for a rectilinear grid or identically shaped two-dimensional coordinates for a curvilinear provider. x_proj and y_proj must be matching two-dimensional projected coordinates.

ModelInput(meteo, grid, geography)

Provider-neutral collection of meteorology, numerical grid, and coordinate metadata. This opt-in container is not accepted by the current pressure-level execution paths.

geometric_height_from_geopotential(geopotential::AbstractArray{<:Real})

Convert geopotential [m^2 s^-2] to curvature-corrected geometric height above the reference ellipsoid [m].

agl_centers_from_geopotential(level_geopotential, surface_geopotential)

Construct curvature-corrected cell-centre heights above local ground level. Level geopotential has shape (Nx, Ny, Nz, Nt) and surface geopotential may have shape (Nx, Ny) or (Nx, Ny, Nt_surface), where Nt_surface is 1 or Nt. Negative AGL values are retained.

vertical_thickness(grid::Grid, t::Integer)

Return vertical cell thickness for one time slice of a provider-neutral Grid, using its AGL cell-centre heights.

load_meteo_era5(levels_file, surface_file; verbose=true) -> Meteorology

Load ERA5 pressure-level meteorology into the provider-neutral Meteorology container. This is opt-in and does not alter legacy solver input routing.

load_grid_era5(levels_file, surface_file; verbose=true) -> Grid

Load ERA5 grid spacing and AGL cell-centre heights into a provider-neutral Grid, using surface geopotential from the single-level file.

load_geography_era5(levels_file; verbose=true) -> Geography

Load oriented ERA5 longitude/latitude metadata and projected coordinates into the provider-neutral Geography container. Surface fields are not required.

load_model_input_era5(levels_file, surface_file; verbose=true) -> ModelInput

Load ERA5 pressure-level and surface files once through the legacy reader, then compose the provider-neutral ModelInput container.

zero_emissions(input::InputFields; T=Float32) -> Array{T,4}

Allocate a zero-filled 4D emissions array on the model grid.

This is the simplest starting point for fully custom emissions or for combining multiple sources with the mutating add_...! helpers.

add_point_source!(emissions, input, source_lon, source_lat;
                  rate_kg_s=1.0,
                  injection_height=0.0,
                  sigma_z=0.0,
                  k0=1,
                  time_indices=:)

Add a point source to an existing 4D emissions array.

This mutating form is useful when combining multiple sources or when building a custom emissions field incrementally.

Arguments

• emissions: Existing (Nx, Ny, Nz, Nt) volumetric emissions array in kg/m^3/s.
• input: InputFields struct defining the grid.
• source_lon, source_lat: Source coordinates in degrees.

Keywords

• rate_kg_s: Emission rate in kg/s (scalar or vector of length Nt).
• injection_height: Injection height in m AGL (scalar or vector). Used when sigma_z > 0.
• sigma_z: Vertical Gaussian sigma for distributing emission. If 0, the emission is injected into one model level.
• k0: Vertical level index for single-level injection (sigma_z=0).
• time_indices: Which model time steps receive the source. Accepts :, a single integer, an integer range/vector, or a Boolean mask of length Nt. Vector-valued inputs may be length Nt or length length(time_indices).

add_ground_source!(emissions, input, source_lon, source_lat;
                   rate_kg_s=1.0, k0=1, time_indices=:)

Add a surface / near-surface source to an existing emissions array.

This is a convenience wrapper around add_point_source! for the common case of emissions released into the first model level.

build_point_source_emissions(input::InputFields, source_lon::Real, source_lat::Real;
                              rate_kg_s::Real=1.0,
                              injection_height::Union{Real,AbstractVector{<:Real}}=0.0,
                              sigma_z::Real=0.0,
                              k0::Int=1,
                              time_indices=:) -> Array{Float32,4}

Build emission field for a point source (e.g., volcano, industrial stack).

Arguments

• input: InputFields struct
• source_lon, source_lat: Source coordinates in degrees
• rate_kg_s: Emission rate in kg/s (scalar or vector of length Nt)
• injection_height: Injection height in m AGL (scalar or vector)
• sigma_z: Vertical Gaussian sigma for distributing emission (0 = single level)
• k0: Vertical level index if sigma_z=0 (default: 1 = surface)
• time_indices: Which model time steps receive the source. Accepts :, a single integer, an integer range/vector, or a Boolean mask of length Nt.

Returns

4D emission array (Nx, Ny, Nz, Nt) in kg/m³/s.

Example

# Constant point source at surface
emissions = build_point_source_emissions(input, -155.5, 19.4;
                                          rate_kg_s=100.0, k0=1)

# Time-varying volcano with injection height
emissions = build_point_source_emissions(input, -155.5, 19.4;
                                          rate_kg_s=rate_t,  # Vector of length Nt
                                          injection_height=injh_t,  # Vector of length Nt
                                          sigma_z=800.0)  # Gaussian spread

build_ground_source_emissions(input::InputFields, source_lon::Real, source_lat::Real;
                              rate_kg_s=1.0, k0=1, time_indices=:) -> Array{Float32,4}

Build a surface / near-surface source on the model grid.

rate_kg_s can be a scalar or a vector of length Nt, which makes this the shortest high-level API for time-varying ground emissions.

build_custom_emissions(builder!, input::InputFields; T=Float32) -> Array{T,4}

Allocate a zero emissions array and hand it to builder! for in-place filling.

This is the intended public entry point for fully custom emissions fields that do not fit one of the built-in source helpers.

Example

emissions = build_custom_emissions(input) do E, input
    dz = calculate_dz_full(input.geopotential)
    i = argmin(abs.(input.longitude .- 15.0))
    j = argmin(abs.(input.latitude .- 37.75))
    rate_t = fill(2.0f0, input.Nt)

    for t in 1:input.Nt
        volume = input.dx[i, j] * input.dy[i, j] * dz[i, j, 1, t]
        E[i, j, 1, t] = rate_t[t] / volume
    end
end

build_wildfire_emissions(input::InputFields, emissions_file::String,
                         flux_var::String="bcfire"; verbose=true) -> Array{Float32,4}

Build 3D wildfire emission field from GFAS surface flux and injection height.

The emissions are distributed vertically based on injection height from the CAMS GFAS dataset, placing the flux into the nearest model level.

Arguments

• input: InputFields struct (provides grid and geopotential)
• emissions_file: Path to CAMS GFAS NetCDF file
• flux_var: Variable name for emission flux (default: "bcfire")
• verbose: Print progress

Returns

4D emission array (Nx, Ny, Nz, Nt) in kg/m³/s.

Example

input = load_input_data("era5_levels.nc", "era5_surface.nc")
emissions = build_wildfire_emissions(input, "gfas.nc", "bcfire")

AerosolInputs

Auxiliary meteorological inputs used by aerosol process parameterizations.

These fields are only needed when you enable WetDeposition(...) inside AerosolProcesses(...), or when you use parameterized DryDeposition(...).

Fields

• large_scale_precip_mmh: Large-scale precipitation rate (Nx, Ny, Nt) in mm/h
• convective_precip_mmh: Convective precipitation rate (Nx, Ny, Nt) in mm/h
• cloud_cover: Cloud cover (Nx, Ny, Nz, Nt) as unitless fraction
• cloud_total_water: Cloud liquid + ice water (Nx, Ny, Nz, Nt) in kg/kg
• friction_velocity: Surface friction velocity u_* (Nx, Ny, Nt) in m/s, or nothing
• lai: Leaf area index (Nx, Ny) as unitless field, or nothing
• land_sea_mask: Land-sea mask (Nx, Ny) as unitless fraction (≥0.5 = land), or nothing

load_aerosol_inputs(levels_file, precip_file;
                    surface_file=nothing,
                    dt_seconds=3600.0,
                    precip_is_cumulative=false,
                    verbose=true) -> AerosolInputs

Load the auxiliary ERA5 fields needed by aerosol wet and dry deposition.

This reads:

• precipitation from precip_file
• cloud cover and cloud water from levels_file
• friction velocity u_* from surface_file if provided

The returned AerosolInputs object can be passed to build_aerosol_simulation(...; aerosol_inputs=...).

AerosolTracer(; radius_m=0.25e-6, density_kg_m3=1500.0,
                CCN=0.9, IN=0.1,
                rain_efficiency=1.0, ice_efficiency=1.0)

Describe a transported aerosol tracer using the particle properties needed by sedimentation and wet deposition.

Arguments

• radius_m: Particle radius in meters
• density_kg_m3: Particle density in kg/m^3
• CCN: Effective cloud condensation nuclei fraction
• IN: Effective ice nuclei fraction
• rain_efficiency: Below-cloud rain scavenging efficiency factor
• ice_efficiency: Below-cloud snow scavenging efficiency factor

WetDeposition(; method=:a19, scale=1.0,
                cloud_replenishment_factor=6.2,
                cloud_drop_radius_m=nothing,
                ctwc_floor_kg_m2=1e-12,
                cloud_base_threshold=0.0,
                precip_threshold_mmh=0.01,
                incloud_scale=1.0,
                belowcloud_scale=1.0,
                belowcloud_scheme=:power_law,
                boost_subgrid_precipitation=true)

Configure the wet-deposition parameterization used by build_aerosol_simulation(...).

Arguments

• method: In-cloud scavenging method, :a19 or :alt
• scale: Global multiplier applied to the final wet-deposition coefficient
• cloud_replenishment_factor: Preferred Eq. A19 parameter i_cr (dimensionless), defaulting to the FLEXPART value 6.2
• cloud_drop_radius_m: Legacy compatibility alias storing cloud_replenishment_factor / rho_water; this is kept for backward compatibility and is not a physical cloud-droplet radius
• ctwc_floor_kg_m2: Floor for precipitating cloud water in the in-cloud denominator
• cloud_base_threshold: Cloud-water threshold used to detect cloud base
• precip_threshold_mmh: Minimum precipitation rate used for below-cloud removal
• incloud_scale: Multiplier applied only to the in-cloud scavenging coefficient
• belowcloud_scale: Multiplier applied only to the below-cloud scavenging coefficient
• belowcloud_scheme: Below-cloud scavenging kernel, one of :power_law, :laakso, or :wang2014. The default :power_law uses a rain/snow precipitation-intensity law.
• boost_subgrid_precipitation: When true, use the published (I_l + I_c) / F precipitation enhancement inside the in-cloud scavenging calculation. The default true follows FLEXPART/Grythe; false keeps the grid-scale precipitation intensity unchanged.

DryDeposition(; vdep_m_s=nothing,
                lai=nothing,
                land_sea_mask=nothing,
                land_a1=4.3e-3,
                water_a1=land_a1,
                coefficient_scale_m_s=1e-2,
                fine_cutoff_m=2.5e-6,
                coarse_cutoff_m=10e-6,
                land_surface=nothing,
                water_surface=nothing)

Configure dry deposition.

Two modes are supported:

• Constant mode: pass vdep_m_s=...
• Parameterized mode: pass lai=... and load u_* through load_aerosol_inputs(...; surface_file=...)

Parameterized mode has two branches:

• Simplified fine-particle mode: omit land_surface and water_surface. TARSA uses V_d = a_1 u_* and supports D_p <= fine_cutoff_m.
• Full Zhang-type mode: pass land_surface=DryDepositionSurface(...) and/or water_surface=.... TARSA then evaluates the fine, coarse and giant branches, including LAI dependence and the super-coarse rebound correction.

The Zhang/He coefficients are tabulated in cm s^-1; TARSA converts them to m s^-1 with coefficient_scale_m_s=1e-2 in the simplified mode. In the full surface mode, the equivalent conversion is carried by each DryDepositionSurface.

If land_sea_mask is supplied, values >= 0.5 are treated as land and values < 0.5 as water. If no mask is supplied, all cells are treated as land.

DryDepositionSurface

Surface-specific coefficients for the Zhang-type aerosol dry-deposition parameterization.

The fields a1, b1:b3, c1:c3, d1:d3, and e1:e3 follow Zhang and He (2014) and are tabulated in cm s^-1; velocity_scale_m_s=1e-2 converts the regressed velocities to m s^-1 by default. lai_max and collector_radius_m are only used for vegetated surfaces. smooth_surface=true disables the LAI adjustment and uses the smooth-surface Stokes-number branch for rebound.

AerosolProcesses(; sedimentation=true,
                   wet_deposition=nothing,
                   dry_deposition=nothing)

Select which linear aerosol processes are applied in build_aerosol_simulation(...).

• sedimentation=true subtracts the settling velocity from the vertical wind
• wet_deposition=WetDeposition(...) adds in-cloud and below-cloud wet removal
• dry_deposition=DryDeposition(...) adds a first-order loss in the bottom layer

build_aerosol_simulation(input, emissions;
                         aerosol=AerosolTracer(),
                         aerosol_inputs=nothing,
                         processes=AerosolProcesses(),
                         solver=ExplicitSolver(),
                         lateral_bc=:neumann0,
                         bottom_bc=:neumann0,
                         initial_bc=nothing) -> Simulation

Build a high-level aerosol simulation with optional sedimentation, wet deposition, and dry deposition.

Use this when you want the public API to assemble aerosol process fields for you, instead of building lambda arrays or modified vertical winds manually.

Example

input = TARSA.load_input_data(levels_file, surface_file)
aerosol_inputs = TARSA.load_aerosol_inputs(levels_file, precip_file; surface_file=surface_file)
emissions = TARSA.zero_emissions(input)
TARSA.add_ground_source!(emissions, input, 15.0, 37.75; rate_kg_s=5.0)

sim = TARSA.build_aerosol_simulation(
    input,
    emissions;
    aerosol=TARSA.AerosolTracer(radius_m=0.11e-6, density_kg_m3=1700.0),
    aerosol_inputs=aerosol_inputs,
    processes=TARSA.AerosolProcesses(
        sedimentation=true,
        wet_deposition=TARSA.WetDeposition(),
        dry_deposition=TARSA.DryDeposition(lai=fill(2.0f0, input.Nx, input.Ny)),
    ),
    solver=TARSA.ExplicitSolver(cfl=0.7),
)

SO2ChemistryInputs

Auxiliary inputs used by the public SO2 chemistry API.

These are separate from InputFields because simplified SO2 chemistry needs metadata and optional fields that are not part of the transport state:

• cloud_cover: ERA5 cloud cover (Nx, Ny, Nz, Nt)
• pressure_levels_pa: pressure-level coordinate (Nz) in Pa
• valid_time: model timestamps (Nt)
• oh_field: optional OH field (Nx, Ny, Nz, Nt) in units compatible with SulfurOxidation.k1_s_per_ppb

load_so2_chemistry_inputs(levels_file; oh_field_name="cams_oh", verbose=true)

Load the auxiliary ERA5/CAMS fields needed by the public SO2 chemistry API.

This reads cloud cover, pressure levels, timestamps, and an optional OH field from the ERA5-levels file.

OHOxidation(; source=:none, scale=1.0, ppb0=NaN,
              sza_power=1.0, alt_p_ref_hpa=1000.0,
              alt_p_top_hpa=200.0, alt_min=1.0)

Configure the simplified gas-phase OH pathway for SO2.

• source=:none disables gas-phase oxidation
• source=:field uses chemistry_inputs.oh_field
• source=:estimated uses a daylight/altitude OH estimate with baseline ppb0
• source=:auto uses chemistry_inputs.oh_field when available, otherwise the estimated OH path if ppb0 > 0

CloudOxidation(; rh_percent=90.0, scale=1.0)

Configure the simplified cloud/aqueous SO2 oxidation term k2(c, RH) used by the public SO2 chemistry API.

SO2Chemistry(; model=:sulfur_k1k2, lifetime_hours=0.0,
               oh=OHOxidation(), cloud=nothing)

Describe the simplified SO2 chemistry applied by build_so2_simulation(...).

build_so2_simulation(input, emissions;
                     chemistry=SO2Chemistry(),
                     chemistry_inputs=nothing,
                     solver=ExplicitSolver(),
                     lateral_bc=:neumann0,
                     bottom_bc=:neumann0,
                     initial_bc=nothing) -> Simulation

Build a high-level SO2 simulation with simplified chemistry.

Use this when you want TARSA to assemble the first-order SO2 loss field from:

• optional lifetime loss
• optional gas-phase OH oxidation
• optional cloud/aqueous oxidation

Example

input = TARSA.load_input_data(levels_file, surface_file)
chem_inputs = TARSA.load_so2_chemistry_inputs(levels_file; oh_field_name="cams_oh")
emissions = TARSA.zero_emissions(input)
TARSA.add_point_source!(emissions, input, -155.287, 19.421;
    rate_kg_s=rate_t,
    injection_height=injh_t,
    sigma_z=800.0,
)

sim = TARSA.build_so2_simulation(
    input,
    emissions;
    chemistry_inputs=chem_inputs,
    chemistry=TARSA.SO2Chemistry(
        lifetime_hours=18.0,
        oh=TARSA.OHOxidation(source=:auto, scale=1.0, ppb0=0.05),
        cloud=TARSA.CloudOxidation(rh_percent=90.0, scale=1.0),
    ),
)

ExplicitSolver(; cfl=0.8f0, scheme=:dst_koren_rk3, kzz_scale=1.0f0, dt=3600.0f0,
                mass_consistent=true)

Create solver configuration for the high-level explicit IMEX solver.

Options

• cfl: CFL number (0.5-0.9 recommended)
• scheme: Advection scheme
  • :dst_koren_rk3 - SSP-RK3 time + Koren limiter (default forward scheme, most accurate)
  • :dst_koren_dimsplit - Koren limiter with directional (dimensional) split; this is the differentiable / inversion operator. (Deprecated alias: :dst_koren_strang.)
  • :upwind - First-order upwind (fastest, more diffusive)
  • :upwind_dimsplit - First-order upwind with directional split. (Deprecated alias: :upwind_strang.)
  ("Strang" alone is reserved for the engine's diffusion–advection–diffusion process split.)
• kzz_scale: Multiplier for turbulent diffusion coefficient
• dt: Model timestep in seconds (usually 3600 for hourly)
• mass_consistent: Whether to use density-weighted (mass-consistent) advection for the :dst_koren family (default true). Set to false for the legacy volume-flux form (e.g. when an AD-consistent inversion forward run is required, since the Enzyme adjoint is not yet implemented for the MC path).

ImplicitSolver(; dt=3600.0f0, kzz_scale=1.0f0, collect_matrices=false)

Create solver configuration for implicit solver.

Options

• dt: Model timestep in seconds
• kzz_scale: Multiplier for turbulent diffusion coefficient
• collect_matrices: Return (mixing_ratio, A_store, volumes_store, dz_array) from run!(sim) for matrix diagnostics

build_simulation(input::InputFields, emissions::AbstractArray{<:Real,4};
                 solver=ExplicitSolver(),
                 lateral_bc::Symbol=:neumann0,
                 bottom_bc::Symbol=:neumann0,
                 initial_bc=nothing,
                 scavenging=nothing,
                 chemistry::Union{Nothing,NamedTuple}=nothing) -> Simulation

Build a Simulation object from input fields and emissions.

Arguments

• input: InputFields struct from loadinputdata()
• emissions: 4D emission array (Nx, Ny, Nz, Nt) in kg/m^3/s. Real-valued arrays are converted to Float32 internally.
• solver: Solver configuration from ExplicitSolver() or ImplicitSolver()
• lateral_bc: :neumann0 (zero flux) or :dirichlet
• bottom_bc: :neumann0 (zero flux) or :dirichlet
• initial_bc: Initial/boundary conditions (default: zeros)
• scavenging: First-order loss coefficient in s^-1 (default: zeros)
• chemistry: Legacy SO2 chemistry tuple. Prefer build_so2_simulation(...) for new SO2 workflows.

Returns

Configured Simulation object ready to run.

Example

input = load_input_data("era5_levels.nc", "era5_surface.nc")
emissions = build_wildfire_emissions(input, "gfas.nc", "bcfire")

sim = build_simulation(input, emissions;
    solver = ExplicitSolver(cfl=0.8),
    lateral_bc = :neumann0,
    bottom_bc = :neumann0
)

result = run!(sim)

run!(sim::Simulation; verbose=true)

Run a simulation and return the simulated mixing-ratio field.

Returns

For normal runs, returns a 4D array (Nx, Ny, Nz, Nt) with mixing ratio in kg/kg. For ImplicitSolver(collect_matrices=true), returns (mixing_ratio, A_store, volumes_store, dz_array).

Example

sim = build_simulation(input, emissions)
mixing_ratio = run!(sim)

# Save output
save_results("output.nc", input, mixing_ratio, emissions)