First Forward Run: Tracer Transport with Real Meteorology

This tutorial demonstrates an end-to-end forward simulation using AtmosTransport.jl with real GEOS-FP meteorological data downloaded directly via OPeNDAP — no authentication required.

By the end, you will have:

Downloaded a coarsened GEOS-FP snapshot from NASA's OPeNDAP server
Built a model grid with 72 vertical levels
Initialized a CO₂-like tracer with a localized perturbation
Run 24 hours of advective transport using Slopes advection
Written the output to a NetCDF file

Prerequisites

Only the main AtmosTransport package and NCDatasets (a standard dependency) are needed. No API keys, no credentials, no extra packages.

using AtmosTransport
using AtmosTransport.Architectures
using AtmosTransport.Grids
using AtmosTransport.Advection
using NCDatasets
using Dates

Step 1: Download GEOS-FP test data

NASA's GEOS Forward Processing system provides near-real-time assimilated meteorological fields at 0.25° × 0.3125° resolution, freely accessible via OPeNDAP. We subsample by a stride of 10 to get a manageable ~3° grid.

const STRIDE = 10
const TEST_DATE = Date(2024, 12, 1)
const OUTDIR = joinpath(homedir(), "data", "metDrivers", "geosfp", "test")
const FT = Float64

Float64

The download function opens the OPeNDAP endpoint, reads a single 12Z snapshot with spatial subsampling, and writes it as a local NetCDF file. Once downloaded, subsequent runs use the cached file automatically.

function download_geosfp(date, outdir; stride=STRIDE)
    mkpath(outdir)
    outfile = joinpath(outdir, "geosfp_asm_$(Dates.format(date, "yyyymmdd")).nc")
    isfile(outfile) && filesize(outfile) > 1000 && return outfile

    url = "https://opendap.nccs.nasa.gov/dods/GEOS-5/fp/0.25_deg/assim/inst3_3d_asm_Nv"
    ds = NCDataset(url)
    try
        lons, lats, levs = ds["lon"][:], ds["lat"][:], ds["lev"][:]
        times = ds["time"][:]
        target = DateTime(date) + Hour(12)
        tidx = if eltype(times) <: DateTime
            argmin([abs(Dates.value(t - target)) for t in times])
        else
            argmin(abs.(times .- (Dates.value(date - Date(1,1,1)) + 0.5)))
        end
        li, la = 1:stride:length(lons), 1:stride:length(lats)

        NCDataset(outfile, "c") do out
            defDim(out, "lon", length(li))
            defDim(out, "lat", length(la))
            defDim(out, "lev", length(levs))
            defVar(out, "lon", Float64, ("lon",))[:] = lons[li]
            defVar(out, "lat", Float64, ("lat",))[:] = lats[la]
            defVar(out, "lev", Float64, ("lev",))[:] = levs
            for v in ["u","v","omega","delp"]
                defVar(out, v, Float32, ("lon","lat","lev"))[:] = Float32.(ds[v][li,la,:,tidx])
            end
        end
    finally
        close(ds)
    end
    return outfile
end

metfile = download_geosfp(TEST_DATE, OUTDIR)

"/home/runner/data/metDrivers/geosfp/test/geosfp_asm_20241201.nc"

Step 2: Load met data and build staggered velocities

The advection scheme requires velocities at cell faces:

u at longitude faces: (Nx+1, Ny, Nz), periodic
v at latitude faces: (Nx, Ny+1, Nz), zero at poles
w at level interfaces: (Nx, Ny, Nz+1), zero at top/surface

We compute these by averaging adjacent cell-center values and apply a CFL limiter to ensure numerical stability.

ds = NCDataset(metfile)
lons, lats = ds["lon"][:], ds["lat"][:]
Nx, Ny, Nz = length(lons), length(lats), length(ds["lev"][:])
u_cc   = FT.(ds["u"][:,:,:])
v_cc   = FT.(ds["v"][:,:,:])
omega  = FT.(ds["omega"][:,:,:])
delp   = FT.(ds["delp"][:,:,:])
close(ds)

closed Dataset

# Column-mean pressure thickness → vertical coordinate
mean_dp = [sum(delp[:,:,k])/(Nx*Ny) for k in 1:Nz]
p_edges = cumsum(vcat(0.0, mean_dp))
Ps = p_edges[end]

# Stagger winds to faces with CFL limiting
Δt = 300.0  # 5-minute time step
cfl = 0.4
Δx_min = STRIDE * 0.3125 * 111000.0
Δy_min = STRIDE * 0.25   * 111000.0
u_max, v_max = cfl * Δx_min / Δt, cfl * Δy_min / Δt

u = zeros(Nx+1, Ny, Nz)
for k in 1:Nz, j in 1:Ny, i in 1:Nx
    u[i,j,k] = clamp((u_cc[i,j,k]+u_cc[i==Nx ? 1 : i+1,j,k])/2, -u_max, u_max)
end
u[Nx+1,:,:] .= u[1,:,:]

v = zeros(Nx, Ny+1, Nz)
for k in 1:Nz, j in 2:Ny, i in 1:Nx
    v[i,j,k] = clamp((v_cc[i,j-1,k]+v_cc[i,j,k])/2, -v_max, v_max)
end

w = zeros(Nx, Ny, Nz+1)
for k in 2:Nz, j in 1:Ny, i in 1:Nx
    w_raw = -(omega[i,j,k-1]+omega[i,j,k])/2
    dp_min = min(mean_dp[k-1], mean_dp[k])
    w[i,j,k] = clamp(w_raw, -cfl*dp_min/Δt, cfl*dp_min/Δt)
end

velocities = (; u, v, w)

(u = [-0.10558640956878662 1.5721479058265686 … 13.093647003173828 19.091601371765137; -0.9793499112129211 0.43457523733377457 … 12.247355461120605 17.604777336120605; … ; 0.4188087582588196 2.224350690841675 … 13.601898193359375 19.965086936950684; -0.10558640956878662 1.5721479058265686 … 13.093647003173828 19.091601371765137;;; 2.0638346076011658 -0.38375914841890335 … 11.516682624816895 18.985431671142578; 1.426577389240265 -1.290219485759735 … 10.913424968719482 17.554049491882324; … ; 2.4446429014205933 0.12497882544994354 … 11.89544153213501 19.82317543029785; 2.0638346076011658 -0.38375914841890335 … 11.516682624816895 18.985431671142578;;; 2.1304547786712646 -4.0073065757751465 … 15.079345703125 20.32085132598877; 1.6254568099975586 -4.783722877502441 … 14.504448890686035 19.783745765686035; … ; 2.431405782699585 -3.555964946746826 … 15.402244567871094 20.620312690734863; 2.1304547786712646 -4.0073065757751465 … 15.079345703125 20.32085132598877;;; … ;;; -0.8081897497177124 -12.160240173339844 … 6.472219467163086 -2.144968032836914; -1.1234923005104065 -12.422808647155762 … 5.714457035064697 -2.1800742149353027; … ; -0.6176870763301849 -11.978038311004639 … 6.950274467468262 -2.1200380325317383; -0.8081897497177124 -12.160240173339844 … 6.472219467163086 -2.144968032836914;;; -0.007020831108093262 -10.414491653442383 … 5.942331790924072 -1.7783713340759277; -0.37555310130119324 -10.799869537353516 … 5.609593868255615 -1.8181405067443848; … ; 0.21413880586624146 -10.192111015319824 … 6.065232276916504 -1.7531270384788513; -0.007020831108093262 -10.414491653442383 … 5.942331790924072 -1.7783713340759277;;; 1.9195311069488525 -6.652734279632568 … 5.202366828918457 -1.4939221143722534; 1.6800243258476257 -7.1031787395477295 … 5.1494832038879395 -1.5331339836120605; … ; 2.062979817390442 -6.3842856884002686 … 5.103298664093018 -1.4689668416976929; 1.9195311069488525 -6.652734279632568 … 5.202366828918457 -1.4939221143722534], v = [0.0 19.37337875366211 … -23.073841094970703 0.0; 0.0 19.55892562866211 … -23.778919219970703 0.0; … ; 0.0 19.069664001464844 … -22.16530132293701 0.0; 0.0 19.32793426513672 … -22.929424285888672 0.0;;; 0.0 13.92320442199707 … -20.721424102783203 0.0; 0.0 14.179281234741211 … -21.147205352783203 0.0; … ; 0.0 13.56017017364502 … -20.187713623046875 0.0; 0.0 13.865142345428467 … -20.63269805908203 0.0;;; 0.0 10.864823341369629 … -9.875043869018555 0.0; 0.0 10.975961685180664 … -10.609586715698242 0.0; … ; 0.0 10.629823684692383 … -8.987610816955566 0.0; 0.0 10.831406593322754 … -9.727540016174316 0.0;;; … ;;; 0.0 6.976498603820801 … -2.672709882259369 0.0; 0.0 6.870480537414551 … -2.778300702571869 0.0; … ; 0.0 7.183637619018555 … -2.264506757259369 0.0; 0.0 7.006543159484863 … -2.626235544681549 0.0;;; 0.0 8.097809076309204 … -2.0820870101451874 0.0; 0.0 8.013580560684204 … -2.2152657210826874 0.0; … ; 0.0 8.260329246520996 … -1.8132467567920685 0.0; 0.0 8.125269174575806 … -2.0456986725330353 0.0;;; 0.0 6.32042121887207 … -1.4047837555408478 0.0; 0.0 6.48033332824707 … -1.7702622711658478 0.0; … ; 0.0 6.266854286193848 … -0.9605422914028168 0.0; 0.0 6.304211616516113 … -1.32133749127388 0.0], w = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 2.3424595383403357e-5 2.471354582667118e-5 … -1.588562736287713e-5 -4.7665374040661845e-5; 2.3424595383403357e-5 1.8201737475465052e-5 … -1.6727542970329523e-5 -4.7665374040661845e-5; … ; 2.3430894998455187e-5 3.068031037400942e-5 … -1.516436259407783e-5 -4.7672596338088624e-5; 2.3425569906976307e-5 2.6010921374108875e-5 … -1.568826769471343e-5 -4.7672571781731676e-5;;; 4.177551454631612e-5 4.281859582988545e-5 … -1.9778065961872926e-5 -7.2189919592347e-5; 4.177551454631612e-5 2.8513480174297e-5 … -2.339304774068296e-5 -7.2189919592347e-5; … ; 4.1784276618273e-5 5.5344333304674365e-5 … -1.6403771041950677e-5 -7.219900726340711e-5; 4.178954986855388e-5 4.551484016701579e-5 … -1.8872556211135816e-5 -7.218042082968168e-5;;; … ;;; 0.09335735812783241 0.049686701968312263 … -0.03265854064375162 0.04035476595163345; 0.09335735812783241 0.08353069797158241 … -0.016701661981642246 0.04035476595163345; … ; 0.09334065020084381 -0.018347961828112602 … 0.016026337631046772 0.04036410618573427; 0.0933496244251728 0.039394546300172806 … -0.02843050006777048 0.040371378883719444;;; 0.0502821309491992 0.042225491255521774 … -0.026487978640943766 0.021272031124681234; 0.0502821309491992 0.062275540083646774 … -0.018767031375318766 0.021272031124681234; … ; 0.050302449613809586 -0.009161049500107765 … 0.020783578045666218 0.021264230366796255; 0.05027688201516867 0.03446877654641867 … -0.021718944888561964 0.021272693295031786;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Step 3: Build model grid

We use a LatitudeLongitudeGrid with 72 hybrid sigma-pressure levels derived from the actual GEOS-FP pressure thickness. This ensures that Δz(k, grid) returns physically realistic pressure thicknesses at each level — crucial for CFL stability.

vc = HybridSigmaPressure(zeros(Nz+1), p_edges ./ Ps)
grid = LatitudeLongitudeGrid(CPU();
    size = (Nx, Ny, Nz),
    longitude = (lons[1], lons[end] + (lons[2]-lons[1])),
    latitude  = (lats[1], lats[end]),
    vertical  = vc)

AtmosTransport.Grids.LatitudeLongitudeGrid{Float64, AtmosTransport.Architectures.CPU, AtmosTransport.Grids.HybridSigmaPressure{Float64}, Vector{Float64}, Nothing}(AtmosTransport.Architectures.CPU(), 116, 73, 72, 3, 3, 1, [-178.4375, -175.3125, -172.1875, -169.0625, -165.9375, -162.8125, -159.6875, -156.5625, -153.4375, -150.3125  …  152.8125, 155.9375, 159.0625, 162.1875, 165.3125, 168.4375, 171.5625, 174.6875, 177.8125, 180.9375], [-180.0, -176.875, -173.75, -170.625, -167.5, -164.375, -161.25, -158.125, -155.0, -151.875  …  154.375, 157.5, 160.625, 163.75, 166.875, 170.0, 173.125, 176.25, 179.375, 182.5], [-88.76712328767124, -86.30136986301369, -83.83561643835617, -81.36986301369862, -78.9041095890411, -76.43835616438355, -73.97260273972603, -71.5068493150685, -69.04109589041096, -66.57534246575344  …  66.57534246575344, 69.04109589041096, 71.5068493150685, 73.97260273972603, 76.43835616438355, 78.9041095890411, 81.36986301369862, 83.83561643835617, 86.30136986301369, 88.76712328767124], [-90.0, -87.53424657534246, -85.06849315068493, -82.6027397260274, -80.13698630136986, -77.67123287671232, -75.20547945205479, -72.73972602739725, -70.27397260273973, -67.8082191780822  …  67.8082191780822, 70.27397260273973, 72.73972602739725, 75.20547945205479, 77.67123287671232, 80.13698630136986, 82.6027397260274, 85.06849315068493, 87.53424657534246, 90.0], 3.125, 2.4657534246575343, AtmosTransport.Grids.HybridSigmaPressure{Float64}([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0351522478708264e-5, 2.3497955162869617e-5, 3.890619885893152e-5, 5.796852588076628e-5, 8.213415258378615e-5, 0.00011355929682697324, 0.00015475007888483299, 0.00020842685084315814, 0.0002779652864824505  …  0.866814313918778, 0.8815909589132285, 0.8963678019873056, 0.9111423967089195, 0.9259151182705274, 0.9406876948401214, 0.9554598550143326, 0.9702302042542923, 0.9850007841240074, 1.0]), 6.371e6, 9.81, 101325.0, nothing, [-178.4375, -175.3125, -172.1875, -169.0625, -165.9375, -162.8125, -159.6875, -156.5625, -153.4375, -150.3125  …  152.8125, 155.9375, 159.0625, 162.1875, 165.3125, 168.4375, 171.5625, 174.6875, 177.8125, 180.9375], [-180.0, -176.875, -173.75, -170.625, -167.5, -164.375, -161.25, -158.125, -155.0, -151.875  …  154.375, 157.5, 160.625, 163.75, 166.875, 170.0, 173.125, 176.25, 179.375, 182.5], [-88.76712328767124, -86.30136986301369, -83.83561643835617, -81.36986301369862, -78.9041095890411, -76.43835616438355, -73.97260273972603, -71.5068493150685, -69.04109589041096, -66.57534246575344  …  66.57534246575344, 69.04109589041096, 71.5068493150685, 73.97260273972603, 76.43835616438355, 78.9041095890411, 81.36986301369862, 83.83561643835617, 86.30136986301369, 88.76712328767124], [-90.0, -87.53424657534246, -85.06849315068493, -82.6027397260274, -80.13698630136986, -77.67123287671232, -75.20547945205479, -72.73972602739725, -70.27397260273973, -67.8082191780822  …  67.8082191780822, 70.27397260273973, 72.73972602739725, 75.20547945205479, 77.67123287671232, 80.13698630136986, 82.6027397260274, 85.06849315068493, 87.53424657534246, 90.0])

Step 4: Initialize tracers

Start with a uniform CO₂-like background (400 ppm) plus a Gaussian perturbation centered over Europe. This makes it easy to visualize how the tracer disperses under the influence of GEOS-FP winds.

c = fill(400.0, Nx, Ny, Nz)
for k in 1:Nz, j in 1:Ny, i in 1:Nx
    Δlon = lons[i] - 10.0
    Δlat = lats[j] - 50.0
    blob = 50.0 * exp(-(Δlon^2/450 + Δlat^2/200 + (k-Nz+5)^2/8))
    c[i,j,k] += blob
end
tracers = (; c)

println("Initial tracer: min=$(minimum(c)), max=$(round(maximum(c), digits=1)), " *
        "mean=$(round(sum(c)/length(c), digits=2))")

Initial tracer: min=400.0, max=450.0, mean=400.05

Step 5: Run 24 hours of advective transport

We use the Russell-Lerner SlopesAdvection scheme (second-order, with minmod flux limiter) and symmetric Strang operator splitting:

advect_x (Δt/2)  →  advect_y (Δt/2)  →  advect_z (Δt/2)
      →  advect_z (Δt/2)  →  advect_y (Δt/2)  →  advect_x (Δt/2)

No convection or diffusion in this demo — pure advection.

scheme = SlopesAdvection(use_limiter=true)
N_hours = 24
N_steps = round(Int, N_hours * 3600 / Δt)
half = Δt / 2

mass_initial = sum(tracers.c)

for step in 1:N_steps
    advect_x!(tracers, velocities, grid, scheme, half)
    advect_y!(tracers, velocities, grid, scheme, half)
    advect_z!(tracers, velocities, grid, scheme, half)
    advect_z!(tracers, velocities, grid, scheme, half)
    advect_y!(tracers, velocities, grid, scheme, half)
    advect_x!(tracers, velocities, grid, scheme, half)

    if step % (N_steps ÷ 8) == 0
        t_hr = step * Δt / 3600
        mass_now = sum(tracers.c)
        Δm = abs(mass_now - mass_initial) / abs(mass_initial) * 100
        println("  t=$(round(t_hr, digits=1))h: " *
                "min=$(round(minimum(tracers.c), digits=1)), " *
                "max=$(round(maximum(tracers.c), digits=1)), " *
                "mass_change=$(round(Δm, sigdigits=3))%")
    end
end

  t=3.0h: min=0.0, max=5544.3, mass_change=0.0667%
  t=6.0h: min=0.0, max=13637.5, mass_change=0.203%
  t=9.0h: min=0.0, max=18784.6, mass_change=0.372%
  t=12.0h: min=0.0, max=23328.7, mass_change=0.569%
  t=15.0h: min=0.0, max=29418.3, mass_change=0.79%
  t=18.0h: min=0.0, max=36953.0, mass_change=1.03%
  t=21.0h: min=0.0, max=46453.4, mass_change=1.28%
  t=24.0h: min=0.0, max=56717.2, mass_change=1.55%

Results

The tracer blob has been transported by the GEOS-FP wind field. Key observations:

Mass conservation: The limiter introduces small mass changes (~1%), which is expected for a single-snapshot wind field that isn't divergence-free after CFL clipping. A proper simulation would use temporally interpolated, divergence-corrected winds.
Dispersion: The initial Gaussian blob spreads both horizontally (by the jet stream and surface winds) and vertically (by omega).
Stability: The SlopesAdvection scheme with CFL-limited winds runs stably for the full 24-hour period.

println("\n--- Final State ---")
println("  Steps: $N_steps")
println("  min=$(round(minimum(tracers.c), digits=2))")
println("  max=$(round(maximum(tracers.c), digits=2))")
println("  mass change: $(round(abs(sum(tracers.c)-mass_initial)/abs(mass_initial)*100, sigdigits=3))%")


--- Final State ---
  Steps: 288
  min=0.0
  max=56717.24
  mass change: 1.55%

Next Steps

From here, you can:

Add physics: Include BoundaryLayerDiffusion and TiedtkeConvection for a more complete transport model
Swap met drivers: Replace geosfp.toml with merra2.toml or era5.toml to compare reanalysis products
Run the adjoint: Use adjoint_advect_x/y/z! to compute sensitivities of an observation to upwind sources
GPU acceleration: Replace CPU() with GPU() for NVIDIA hardware (requires the CUDA extension)
Higher resolution: Reduce STRIDE for finer grids (needs proportionally smaller Δt)

This page was generated using Literate.jl.