diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 000000000..a02d45739 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,321 @@ +# ClimaOcean.jl rules for agent-assisted development + +## Project Overview + +ClimaOcean.jl is a Julia package for realistic ocean-only and coupled ocean + sea-ice simulations driven by prescribed atmospheres. It is built on top of [Oceananigans.jl](https://github.com/CliMA/Oceananigans.jl) and [ClimaSeaIce.jl](https://github.com/CliMA/ClimaSeaIce.jl). + +Key features include: +- **OceanSeaIceModel**: Core abstraction for coupling ocean, sea ice, and prescribed atmosphere +- **PrescribedAtmosphere**: Interface for driving ocean simulations with reanalysis data (JRA55) +- **Data Wrangling**: Utilities for downloading and processing bathymetry (ETOPO), ocean state (ECCO, GLORYS, EN4), and atmospheric forcing (JRA55) +- **Surface Flux Computation**: Similarity theory and coefficient-based turbulent flux calculations +- **ocean_simulation**: Convenience constructor for realistic ocean simulations with sensible defaults + +### Relationship to Oceananigans + +ClimaOcean provides the "realistic modeling" layer on top of Oceananigans: +- Oceananigans handles the fluid dynamics (solvers, grids, fields, operators) +- ClimaOcean handles coupling, forcing, and data wrangling for realistic simulations +- Users should be proficient with Oceananigans to use ClimaOcean effectively + +## Language & Environment + +- **Language**: Julia 1.10+ +- **Architectures**: CPU and GPU (via Oceananigans/KernelAbstractions.jl) +- **Key Dependencies**: + - `Oceananigans.jl` - Core fluid dynamics + - `ClimaSeaIce.jl` - Sea ice model + - `SeawaterPolynomials.jl` - TEOS-10 equation of state + - `SurfaceFluxes.jl` - Surface flux parameterizations (via Thermodynamics.jl) +- **Testing**: Test suite covering data downloading, surface fluxes, and coupled models + +## Code Style & Conventions + +### Julia Best Practices + +1. **Explicit Imports**: Use explicit imports in source code + - Import from Oceananigans modules explicitly + - Example: `using Oceananigans.Utils: launch!` + +2. **Type Stability**: Prioritize type-stable code for GPU performance + - All structs must be concretely typed + +3. **Kernel Functions**: For GPU compatibility: + - Use KernelAbstractions.jl syntax (`@kernel`, `@index`) + - Keep kernels type-stable and allocation-free + - Use `ifelse` instead of short-circuiting if-statements when possible + - No error messages inside kernels + - Models _never_ go inside kernels + +4. **Documentation**: + - Use DocStringExtensions.jl for consistent docstrings + - Include `$(SIGNATURES)` for automatic signature documentation + - Add examples in docstrings when helpful + +5. **Memory Efficiency**: + - Favor inline computations over temporary allocations + - Design solutions that work within the existing Oceananigans framework + +### Naming Conventions + +- **Files**: snake_case (e.g., `ocean_sea_ice_model.jl`, `similarity_theory_turbulent_fluxes.jl`) +- **Types**: PascalCase (e.g., `OceanSeaIceModel`, `PrescribedAtmosphere`, `SimilarityTheoryFluxes`) +- **Functions**: snake_case (e.g., `ocean_simulation`, `compute_atmosphere_ocean_fluxes!`) +- **Kernels**: May be prefixed with underscore (e.g., `_compute_flux_kernel`) +- **Variables**: Use _either_ an English long name or mathematical notation with readable unicode + +### Module Structure + +``` +src/ +├── ClimaOcean.jl # Main module, exports +├── Bathymetry.jl # Bathymetry regridding utilities +├── SeaIceSimulations.jl # Sea ice simulation setup +├── OceanSimulations/ # Ocean simulation constructors +│ ├── OceanSimulations.jl +│ ├── ocean_simulation.jl # Hydrostatic ocean simulation setup +│ └── nonhydrostatic_ocean_simulation.jl +├── OceanSeaIceModels/ # Coupled model implementation +│ ├── OceanSeaIceModels.jl # Module definition +│ ├── ocean_sea_ice_model.jl # Core model type +│ ├── PrescribedAtmospheres.jl # Atmosphere data structures +│ ├── time_step_ocean_sea_ice_model.jl +│ └── InterfaceComputations/ # Flux calculations +│ ├── interface_states.jl +│ ├── similarity_theory_turbulent_fluxes.jl +│ ├── coefficient_based_turbulent_fluxes.jl +│ ├── atmosphere_ocean_fluxes.jl +│ └── ... +├── DataWrangling/ # Data downloading and processing +│ ├── DataWrangling.jl # Module definition +│ ├── metadata.jl # Metadata types +│ ├── inpainting.jl # Gap-filling algorithms +│ ├── restoring.jl # Dataset restoring forcing +│ ├── ECCO/ # ECCO data utilities +│ ├── JRA55/ # JRA55 atmosphere data +│ ├── ETOPO/ # Bathymetry data +│ ├── Copernicus/ # GLORYS data (via extension) +│ └── EN4/ # EN4 data +├── InitialConditions/ # Initial condition utilities +│ └── diffuse_tracers.jl +└── Diagnostics/ # Diagnostic utilities + └── mixed_layer_depth.jl +``` + +## Testing Guidelines + +### Running Tests + +```julia +# All tests +Pkg.test("ClimaOcean") + +# Run specific test groups by setting TEST_GROUP environment variable +ENV["TEST_GROUP"] = "fluxes" # or "JRA55", "ecco2_monthly", "bathymetry", etc. +Pkg.test("ClimaOcean") + +# Available test groups: +# - init: Download test data +# - JRA55: JRA55 atmosphere tests +# - ecco2_monthly, ecco2_daily, ecco4_en4: Dataset tests +# - downloading, copernicus_downloading: Download tests +# - fluxes: Surface flux tests +# - bathymetry: Bathymetry regridding tests +# - ocean_sea_ice_model: Coupled model tests +# - distributed: MPI/distributed tests +# - reactant: Reactant.jl integration tests +``` + +### Writing Tests + +- Place tests in `test/` directory +- Follow the existing test group structure in `runtests.jl` +- Use `include("runtests_setup.jl")` for common setup +- Test on both CPU and GPU when possible +- Name test files `test_.jl` + +### Data Dependencies + +- Tests download data from external sources (ECCO, JRA55, ETOPO) +- Use `@get_scratch!` for caching downloaded files +- Be mindful of data download times in CI + +## Common Development Tasks + +### Adding New Data Sources + +1. Create a new subdirectory in `src/DataWrangling/` (e.g., `NewDataset/`) +2. Define dataset types and metadata structures +3. Implement required interface functions: + - `download_dataset(metadata)` + - `z_interfaces(dataset)` + - `native_grid(dataset)` +4. Add to `DataWrangling.jl` includes +5. Export from `ClimaOcean.jl` if needed +6. Add tests for downloading and usage + +### Modifying Surface Flux Calculations + +- Flux calculations are in `src/OceanSeaIceModels/InterfaceComputations/` +- `SimilarityTheoryFluxes` uses Monin-Obukhov similarity theory +- `CoefficientBasedFluxes` uses simple transfer coefficients +- Roughness lengths: `roughness_lengths.jl` +- State interpolation: `interpolate_atmospheric_state.jl` + +### Adding New Coupled Model Features + +1. Core model: `src/OceanSeaIceModels/ocean_sea_ice_model.jl` +2. Time stepping: `time_step_ocean_sea_ice_model.jl` +3. Interface fluxes: `InterfaceComputations/` +4. Update exports in `OceanSeaIceModels.jl` and `ClimaOcean.jl` + +## Documentation + +### Building Docs Locally + +```sh +julia --project=docs/ docs/make.jl +``` + +### Viewing Docs + +```julia +using LiveServer +serve(dir="docs/build") +``` + +### Documentation Style + +- Use Documenter.jl syntax for cross-references +- Include code examples in documentation pages +- Add references to papers in `climaocean.bib` +- In examples, use `using ClimaOcean` and avoid explicit imports of exported names + +### Writing Examples + +- Explain at the top of the file what a simulation is doing +- Use Literate.jl style - let code speak for itself +- Common patterns: + - Single-column simulations (e.g., `single_column_os_papa_simulation.jl`) + - Regional simulations (e.g., `mediterranean_simulation_with_ecco_restoring.jl`) + - Near-global simulations (e.g., `near_global_ocean_simulation.jl`) +- Demonstrate data wrangling, model setup, and visualization + +## Important Files to Know + +### Core Implementation + +- `src/ClimaOcean.jl` - Main module, all exports +- `src/OceanSeaIceModels/ocean_sea_ice_model.jl` - Core coupled model type +- `src/OceanSeaIceModels/PrescribedAtmospheres.jl` - Atmosphere data structures +- `src/OceanSimulations/ocean_simulation.jl` - Ocean simulation constructor +- `src/OceanSeaIceModels/InterfaceComputations/` - Surface flux calculations + +### Data Wrangling + +- `src/DataWrangling/metadata.jl` - Metadata types and interface +- `src/DataWrangling/JRA55/JRA55_prescribed_atmosphere.jl` - JRA55 atmosphere +- `src/DataWrangling/ECCO/ECCO.jl` - ECCO data utilities +- `src/Bathymetry.jl` - Bathymetry regridding + +### Configuration + +- `Project.toml` - Package dependencies and compat bounds +- `test/runtests.jl` - Test configuration + +### Examples + +- `examples/single_column_os_papa_simulation.jl` - Single column at Ocean Station Papa +- `examples/mediterranean_simulation_with_ecco_restoring.jl` - Regional with restoring +- `examples/near_global_ocean_simulation.jl` - Near-global hydrostatic +- `examples/diurnal_large_eddy_simulation.jl` - LES with diurnal forcing +- `examples/idealized_single_column_simulation.jl` - Idealized single column + +## Physics Domain Knowledge + +### Air-Sea Fluxes + +- Momentum flux (wind stress): τ = ρₐ Cᴰ |Δu| Δu +- Sensible heat flux: Qₛ = ρₐ cₚ Cᴴ |Δu| ΔT +- Latent heat flux: Qₗ = ρₐ Lᵥ Cᴱ |Δu| Δq +- Monin-Obukhov similarity theory for stability-dependent transfer coefficients +- Roughness length parameterizations (Charnock, COARE, etc.) + +### Ocean-Sea Ice Coupling + +- Freezing-limited ocean temperature +- Sea ice thermodynamics via ClimaSeaIce.jl +- Ice-ocean heat and salt fluxes +- Ice concentration and thickness tracking + +### Radiation + +- Downwelling shortwave and longwave radiation +- Albedo parameterizations (latitude-dependent, tabulated) +- Surface temperature feedback + +### Datasets + +- **JRA55**: 3-hourly atmospheric reanalysis (1958-present) +- **ECCO**: Monthly/daily ocean state estimates +- **GLORYS**: Copernicus ocean reanalysis (via extension) +- **ETOPO**: Global bathymetry +- **EN4**: Subsurface ocean temperature/salinity + +## Common Pitfalls + +1. **Type Instability**: Especially in kernel functions - check with `@code_warntype` +2. **Missing Data**: ECCO/GLORYS data needs inpainting for land points +3. **Coordinate Systems**: Be careful with longitude conventions (0-360 vs -180-180) +4. **Time Zones**: JRA55 uses UTC; be consistent with dates +5. **Grid Staggering**: Remember Oceananigans C-grid locations for fluxes +6. **Memory**: Large FieldTimeSeries can consume significant memory; use `InMemory(chunk_size)` or `OnDisk()` + +## Git Workflow + +- Follow ColPrac (Collaborative Practices for Community Packages) +- Create feature branches for new work +- Write descriptive commit messages +- Update tests and documentation with code changes +- Check CI passes before merging + +## Helpful Resources + +- ClimaOcean docs: https://clima.github.io/ClimaOceanDocumentation/stable/ +- Oceananigans docs: https://clima.github.io/OceananigansDocumentation/stable/ +- ClimaSeaIce: https://github.com/CliMA/ClimaSeaIce.jl +- Discussions: https://github.com/CliMA/Oceananigans.jl/discussions +- JRA55 documentation: https://jra.kishou.go.jp/JRA-55/ +- ECCO data portal: https://ecco-group.org/ +- ETOPO: https://www.ncei.noaa.gov/products/etopo-global-relief-model + +## When Unsure + +1. Check existing examples in `examples/` directory +2. Look at similar implementations in the codebase +3. Review tests for usage patterns +4. Ask in GitHub discussions +5. Check Oceananigans documentation for underlying functionality +6. Review the ClimaOcean preprint for high-level overview + +## AI Assistant Behavior + +- Prioritize type stability and GPU compatibility +- Follow established patterns in existing ClimaOcean code +- Add tests for new functionality +- Update exports in `ClimaOcean.jl` when adding public API +- Consider data downloading and caching implications +- Reference physics equations in comments when implementing flux calculations +- Maintain consistency with both ClimaOcean and Oceananigans coding styles + +## Current Development Focus + +Active areas of development to be aware of: + +- **Reactant.jl integration**: XLA compilation for performance (see `ext/ClimaOceanReactantExt.jl`) +- **Copernicus/GLORYS support**: Alternative to ECCO for ocean state (see `ext/ClimaOceanCopernicusMarineExt.jl`) +- **Nonhydrostatic ocean simulations**: LES with realistic forcing +- **Improved sea ice coupling**: Better ice-ocean interactions +- **Additional datasets**: EN4, multi-year JRA55, etc. +- **Surface flux improvements**: Better roughness length parameterizations +- **Distributed computing**: MPI support for large simulations diff --git a/docs/make.jl b/docs/make.jl index 99afa9e46..1a28d8616 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -18,9 +18,10 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/literated") to_be_literated = [ - "single_column_os_papa_simulation.jl", - "one_degree_simulation.jl", - "near_global_ocean_simulation.jl" + # "single_column_os_papa_simulation.jl", + # "one_degree_simulation.jl", + # "near_global_ocean_simulation.jl" + "diurnal_large_eddy_simulation.jl" ] for file in to_be_literated @@ -44,9 +45,10 @@ pages = [ "Home" => "index.md", "Examples" => [ - "Single-column ocean simulation" => "literated/single_column_os_papa_simulation.md", - "One-degree ocean--sea ice simulation" => "literated/one_degree_simulation.md", - "Near-global ocean simulation" => "literated/near_global_ocean_simulation.md", + # "Single-column ocean simulation" => "literated/single_column_os_papa_simulation.md", + # "One-degree ocean--sea ice simulation" => "literated/one_degree_simulation.md", + # "Near-global ocean simulation" => "literated/near_global_ocean_simulation.md", + "Diurnally-varying ocean large eddy simulation" => "literated/diurnal_large_eddy_simulation.md", ], "Vertical grids" => "vertical_grids.md", diff --git a/examples/diurnal_large_eddy_simulation.jl b/examples/diurnal_large_eddy_simulation.jl new file mode 100644 index 000000000..74c052a97 --- /dev/null +++ b/examples/diurnal_large_eddy_simulation.jl @@ -0,0 +1,331 @@ +# # Diurnal large eddy simulation +# +# This script assembles a large-eddy simulation that resolves the marine boundary layer +# under a clear diurnal cycle, transitioning from nighttime cooling to daytime heating. +# The setup is based on typical summer solstice conditions near San Diego (33°N latitude). +# +# The atmosphere is prescribed with time-varying temperature, humidity, wind, and +# radiation fields that follow a realistic diurnal pattern. + +using ClimaOcean +using Oceananigans +using Oceananigans.Units +using Printf +using CairoMakie + +# ## Parameters +# +# Physical and numerical parameters for the simulation. +# Solar radiation parameters: +# The solar constant (top-of-atmosphere irradiance) is ~1361 W/m². +# Clear-sky atmospheric transmittance is about 75%, so surface irradiance is ~1020 W/m². +# At 33°N on summer solstice, the solar elevation at noon is: +# +# ``θ = 90° - |ϕ - δ| = 90° - |33° - 23.5°| ≈ 80.5°`` +# +# where ``ϕ`` is latitude and ``δ`` is the solar declination (23.5° at summer solstice). +# Thus peak surface irradiance is approximately ``1361 × sin(80.5°) × 0.75 ≈ 1006`` W/m². + +day = 86400 # seconds in a day +latitude = 33 # degrees (San Diego) +solar_constant = 1361 # W/m², top-of-atmosphere +clear_sky_transmittance = 0.75 # atmospheric transmittance on clear day +surface_solar_irradiance = solar_constant * clear_sky_transmittance # ≈ 1020 W/m² + +# ## Diurnal cycle functions +# +# These functions describe the diurnal evolution of atmospheric state variables +# based on typical clear-sky summer solstice conditions for coastal Southern California. +# +# The simulation starts at midnight (t = 0). Key features: +# - Sunrise around 6 AM, sunset around 8 PM (14-hour day) +# - Maximum solar radiation ~1000 W/m² at solar noon (see solar parameters above) +# - Maximum temperature around 2-3 PM (thermal lag) +# - Minimum temperature around 6 AM (just before sunrise) + +""" + solar_elevation(t) + +Compute the solar elevation angle (radians) at time `t` for a given `latitude`. +Simplified model assuming summer solstice conditions. +""" +function solar_elevation(t) + hour_angle = 2π * (t / day - 0.5) # 0 at noon + δ = deg2rad(23.5) # summer solstice declination + ϕ = deg2rad(latitude) + sin_elevation = sin(ϕ) * sin(δ) + cos(ϕ) * cos(δ) * cos(hour_angle) + return asin(clamp(sin_elevation, -1, 1)) +end + +""" + diurnal_shortwave(t) + +Clear-sky downwelling shortwave radiation (W/m²) at time `t`. +""" +function diurnal_shortwave(t) + elevation = solar_elevation(t) + daytime_shortwave = surface_solar_irradiance * sin(elevation) + return max(0, daytime_shortwave) +end + +""" + diurnal_longwave(t, day) + +Downwelling longwave radiation (W/m²) at time `t`. +Varies with atmospheric temperature: ~300 W/m² at night, ~380 W/m² during day. +""" +function diurnal_longwave(t) + phase = 2π * t / day - π/2 - 2 * 2π / 24 + return 340 + 40 * sin(phase) +end + +""" + diurnal_temperature(t, day) + +Air temperature (Kelvin) at time `t`. +Ranges from ~17°C (290 K) before sunrise to ~24°C (297 K) in mid-afternoon. +""" +function diurnal_temperature(t) + T_min, T_max = 290, 297 + T_mean = (T_min + T_max) / 2 + T_amp = (T_max - T_min) / 2 + phase = 2π * (t / day - 6 / 24) - π/2 + return T_mean + T_amp * sin(phase) +end + +""" + diurnal_humidity(t, day) + +Specific humidity (kg/kg) at time `t`. +Higher at night (~0.012), lower during day (~0.008). +""" +function diurnal_humidity(t) + q_mean, q_amp = 0.010, 0.002 + phase = 2π * (t / day - 6 / 24) - π/2 + return q_mean - q_amp * sin(phase) +end + +""" + diurnal_wind_u(t, day) + +Eastward wind component (m/s) at time `t`. +Land-sea breeze pattern: offshore at night, onshore during day. +""" +function diurnal_wind_u(t) + phase = 2π * t / day - π + return 2 + 3 * sin(phase) +end + +diurnal_wind_v(t) = -1 # slight southward component + +# ## Ocean LES configuration +# +# 3D grid with 4m spacing: 512m × 4m × 256m (Nx × Ny × Nz = 128 × 1 × 64). +# The thin y-dimension makes this essentially a 2D slice for computational efficiency. + +arch = GPU() +Lx, Ly, Lz = 512, 512, 256 # meters +Nx, Ny, Nz = 128, 128, 64 # 4m grid spacing + +grid = RectilinearGrid(arch; size = (Nx, Ny, Nz), topology = (Periodic, Periodic, Bounded), + x = (0, Lx), y = (0, Ly), z = (-Lz, 0)) + +coriolis = FPlane(; latitude) +ocean = nonhydrostatic_ocean_simulation(grid; coriolis) + +# Initial conditions: warm mixed layer with small perturbations + +Tᵢ(x, y, z) = 22 + 0.01 * randn() # °C, typical summer SST +Sᵢ(x, y, z) = 33 + 0.001 * randn() # psu, typical coastal salinity +set!(ocean.model, T=Tᵢ, S=Sᵢ) + +# ## Prescribed diurnal atmosphere +# +# Set up the atmosphere with time-varying fields over a full diurnal cycle. +# We use hourly time snapshots for smooth interpolation. + +atmosphere_grid = RectilinearGrid(arch; size = (Nx, Ny), topology = (Periodic, Periodic, Flat), + x = (0, Lx), y = (0, Ly)) + +atmosphere_times = range(0, 24hours, length=25) +atmosphere = PrescribedAtmosphere(atmosphere_grid, collect(atmosphere_times)) + +# Use set! to fill the atmosphere with diurnal cycle functions. +# We use closures to capture the parameters. + +set!(atmosphere; + u = diurnal_wind_u, + v = diurnal_wind_v, + T = diurnal_temperature, + q = diurnal_humidity, + shortwave = diurnal_shortwave, + longwave = diurnal_longwave) + +# ## Coupled ocean–atmosphere model + +coupled_model = OceanSeaIceModel(ocean; atmosphere) +simulation = Simulation(coupled_model; Δt=1second, stop_time=24hours) + +# ## Output: snapshots and time-averages +# +# We save: +# 1. Snapshots of velocity and temperature fields for visualization +# 2. Time-averaged profiles of u and T (horizontally averaged) +# 3. Time-series of surface fluxes (both pointwise and horizontally-averaged) + +u, v, w = ocean.model.velocities +T = ocean.model.tracers.T +Q = coupled_model.interfaces.net_fluxes.ocean_surface.Q +τx = coupled_model.interfaces.net_fluxes.ocean_surface.u +τy = coupled_model.interfaces.net_fluxes.ocean_surface.v + +# Snapshot output every 10 minutes + +Nz = grid.Nz +uˢ = view(u, :, :, Nz) +vˢ = view(v, :, :, Nz) +Tˢ = view(T, :, :, Nz) +u_avg = Average(u, dims=(1, 2)) +T_avg = Average(T, dims=(1, 2)) +Q_avg = Average(Q, dims=(1, 2)) +τx_avg = Average(τx, dims=(1, 2)) + +outputs = (; uˢ, vˢ, Tˢ, Q, τx, τy, u_avg, T_avg, Q_avg, τx_avg) + +simulation.output_writers[:snapshots] = JLD2Writer(ocean.model, outputs; + filename = "diurnal_les", + schedule = TimeInterval(10minutes), + overwrite_existing = true) + +# Progress callback + +function progress(sim) + t = time(sim) + hour = t / 3600 + u, v, w = sim.model.ocean.model.velocities + T = sim.model.ocean.model.tracers.T + Tmax, Tmin = maximum(interior(T)), minimum(interior(T)) + umax, wmax = maximum(abs, u), maximum(abs, w) + Q = sim.model.interfaces.net_fluxes.ocean_surface.Q + Qmax, Qmin = maximum(Q), minimum(Q) + @info @sprintf("Hour %5.1f | Δt: %s | SST: %.2f–%.2f°C | Q: %.0f–%.0f W/m² | max|u|: %.2e, max|w|: %.2e", + hour, prettytime(sim.Δt), Tmin, Tmax, Qmin, Qmax, umax, wmax) +end + +add_callback!(simulation, progress, IterationInterval(100)) + +# ## Run the simulation + +run!(simulation) + +# ## Animation +# +# Create a multi-panel animation showing: +# - Main panels: x-z slices of u and T +# - Top panels: time-series of momentum and heat flux +# - Side panels: time-averaged profiles of u and T + +using JLD2 + +# Load data + +u_snapshots = FieldTimeSeries("diurnal_les.jld2", "uˢ") +T_snapshots = FieldTimeSeries("diurnal_les.jld2", "Tˢ") +Q_avg_ts = FieldTimeSeries("diurnal_les.jld2", "Q_avg") +τx_avg_ts = FieldTimeSeries("diurnal_les.jld2", "τx_avg") +u_profiles = FieldTimeSeries("diurnal_les.jld2", "u_avg") +T_profiles = FieldTimeSeries("diurnal_les.jld2", "T_avg") + +times = u_snapshots.times +Nt = length(times) + +# Compute time-averaged profiles + +u_mean = zeros(Nz) +T_mean = zeros(Nz) +for n in 1:Nt + u_mean .+= interior(u_profiles[n], 1, 1, :) + T_mean .+= interior(T_profiles[n], 1, 1, :) +end +u_mean ./= Nt +T_mean ./= Nt + +# Extract coordinates + +xc = xnodes(u_snapshots) +zc = znodes(u_snapshots) + +# Build time-series arrays for flux plots + +hours = times ./ 3600 +Q_avg_series = [interior(Q_avg_ts[n], 1, 1, 1) for n in 1:Nt] +τx_avg_series = [interior(τx_avg_ts[n], 1, 1, 1) for n in 1:Nt] + +# Create figure + +fig = Figure(size = (1400, 900), fontsize = 14) + +# Layout: +# Row 1: flux time-series (spans columns 2-3) +# Row 2: u profile | u heatmap | T heatmap | T profile + +ax_flux = Axis(fig[1, 2:3], xlabel = "Hour", ylabel = "Flux", + title = "Surface fluxes") + +ax_u_profile = Axis(fig[2, 1], xlabel = "⟨u⟩ (m/s)", ylabel = "z (m)", + title = "Time-avg u") + +ax_u = Axis(fig[2, 2], xlabel = "x (m)", ylabel = "z (m)", + title = "Horizontal velocity u") + +ax_T = Axis(fig[2, 3], xlabel = "x (m)", ylabel = "z (m)", + title = "Temperature T") + +ax_T_profile = Axis(fig[2, 4], xlabel = "⟨T⟩ (°C)", ylabel = "z (m)", + title = "Time-avg T") + +# Time-averaged profile plots (static) + +lines!(ax_u_profile, u_mean, zc, color = :blue, linewidth = 2) +lines!(ax_T_profile, T_mean, zc, color = :red, linewidth = 2) + +# Flux time-series (will update with vertical line marker) + +lines!(ax_flux, hours, Q_avg_series, label = "Heat flux Q (W/m²)", color = :red) +lines!(ax_flux, hours, τx_avg_series .* 1000, label = "Mom. flux τx (×10³ N/m²)", color = :blue) +axislegend(ax_flux, position = :rt) + +time_marker = Observable(0) +vlines!(ax_flux, time_marker, color = :black, linestyle = :dash, linewidth = 2) + +# Heatmap observables + +n = Observable(1) + +u_slice = @lift interior(u_snapshots[$n], :, 1, :) +T_slice = @lift interior(T_snapshots[$n], :, 1, :) + +# Determine color limits from data + +u_lim = maximum(abs, interior(u_snapshots[Nt])) +T_lim_min, T_lim_max = minimum(interior(T_snapshots[1])), maximum(interior(T_snapshots[Nt])) + +hm_u = heatmap!(ax_u, xc, zc, u_slice, colormap = :balance, colorrange = (-u_lim, u_lim)) +hm_T = heatmap!(ax_T, xc, zc, T_slice, colormap = :thermal, colorrange = (T_lim_min - 0.5, T_lim_max + 0.5)) + +Colorbar(fig[3, 2], hm_u, vertical = false, label = "u (m/s)") +Colorbar(fig[3, 3], hm_T, vertical = false, label = "T (°C)") + +# Title with time + +title = @lift "Diurnal LES: Hour " * @sprintf("%.1f", times[$n] / 3600) +Label(fig[0, :], title, fontsize = 20, font = :bold) + +# Record animation + +record(fig, "diurnal_les_animation.mp4", 1:Nt; framerate = 12) do i + n[] = i + time_marker[] = times[i] / 3600 +end + +# ![Diurnal LES Animation](diurnal_les_animation.mp4) diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl index 9f07663fc..5c59b4a82 100644 --- a/examples/single_column_os_papa_simulation.jl +++ b/examples/single_column_os_papa_simulation.jl @@ -153,7 +153,7 @@ Q = ρₒ * cₚ * JT ρτx = ρₒ * τx ρτy = ρₒ * τy N² = buoyancy_frequency(ocean.model) -κc = ocean.model.diffusivity_fields.κc +κc = ocean.model.closure_fields.κc fluxes = (; ρτx, ρτy, E, Js, Qv, Qc) auxiliary_fields = (; N², κc) diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 7ac5f01cb..84ed4b05e 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -41,6 +41,7 @@ export LinearlyTaperedPolarMask, DatasetRestoring, ocean_simulation, + nonhydrostatic_ocean_simulation, sea_ice_simulation, initialize! diff --git a/src/DataWrangling/restoring.jl b/src/DataWrangling/restoring.jl index 5355acfbb..92a75250d 100644 --- a/src/DataWrangling/restoring.jl +++ b/src/DataWrangling/restoring.jl @@ -11,7 +11,7 @@ using NCDatasets using Dates: Second import ClimaOcean: stateindex -import Oceananigans.Forcings: regularize_forcing +import Oceananigans.Forcings: materialize_forcing # Variable names for restorable data struct Temperature end @@ -230,7 +230,7 @@ function Base.show(io::IO, dsr::DatasetRestoring) "└── native_grid: ", summary(dsr.native_grid)) end -regularize_forcing(forcing::DatasetRestoring, field, field_name, model_field_names) = forcing +materialize_forcing(forcing::DatasetRestoring, field, field_name, model_field_names) = forcing ##### ##### Masks for restoring diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/OceanSeaIceModels/PrescribedAtmospheres.jl index e7a3ca18c..8d58ba756 100644 --- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl +++ b/src/OceanSeaIceModels/PrescribedAtmospheres.jl @@ -10,6 +10,7 @@ using Oceananigans.Utils: prettysummary, Time using Adapt using Thermodynamics.Parameters: AbstractThermodynamicsParameters +import Oceananigans.Fields: set! import Oceananigans.TimeSteppers: time_step!, update_state! import Thermodynamics.Parameters: @@ -445,4 +446,74 @@ Adapt.adapt_structure(to, tsdr::TwoBandDownwellingRadiation) = TwoBandDownwellingRadiation(adapt(to, tsdr.shortwave), adapt(to, tsdr.longwave)) +##### +##### set! for PrescribedAtmosphere with time-dependent functions +##### + +# Map field names to their corresponding FieldTimeSeries in the atmosphere. +# Follows the pattern: check velocities, tracers, pressure, radiation, freshwater_flux. +function get_atmosphere_field(atmos::PrescribedAtmosphere, name::Symbol) + name === :p && return atmos.pressure + name ∈ propertynames(atmos.velocities) && return getproperty(atmos.velocities, name) + name ∈ propertynames(atmos.tracers) && return getproperty(atmos.tracers, name) + name ∈ propertynames(atmos.downwelling_radiation) && return getproperty(atmos.downwelling_radiation, name) + name ∈ propertynames(atmos.freshwater_flux) && return getproperty(atmos.freshwater_flux, name) + + error("Unknown atmosphere field name: $name. " * + "Available names: :u, :v, :T, :q, :p, :shortwave, :longwave, :rain, :snow") +end + +""" + set!(atmos::PrescribedAtmosphere; kwargs...) + +Set fields of a `PrescribedAtmosphere` using functions of time. + +Each keyword argument should be a function that takes a single argument `t` (time in seconds) +and returns the field value at that time. The function is evaluated at each time in +`atmos.times` and the resulting values are used to fill the corresponding `FieldTimeSeries`. + +Available field names: +- `:u`, `:v` - velocity components (m/s) +- `:T` - temperature (K) +- `:q` - specific humidity (kg/kg) +- `:p` - pressure (Pa) +- `:shortwave` - downwelling shortwave radiation (W/m²) +- `:longwave` - downwelling longwave radiation (W/m²) +- `:rain`, `:snow` - precipitation rates (m/s) + +Example +======= + +```jldoctest +using Oceananigans +using ClimaOcean + +grid = RectilinearGrid(size=(4, 4), x=(0, 1), y=(0, 1), topology=(Periodic, Periodic, Flat)) +times = [0.0, 3600.0, 7200.0] # 0, 1, and 2 hours in seconds +atmosphere = PrescribedAtmosphere(grid, times) + +# Set temperature as a linear function of time +set!(atmosphere; T = t -> 300 + t / 3600) # increases by 1 K per hour + +atmosphere.tracers.T[1], atmosphere.tracers.T[2], atmosphere.tracers.T[3] + +# output +(300.0, 301.0, 302.0) +``` +""" +function set!(atmos::PrescribedAtmosphere; kwargs...) + times = atmos.times + + for (name, func) in pairs(kwargs) + fts = get_atmosphere_field(atmos, name) + for (n, t) in enumerate(times) + parent(fts)[n] = func(t) + end + end + + update_state!(atmos) + + return nothing +end + end # module diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 9c3bfabca..e5a56bb46 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -1,6 +1,7 @@ module OceanSimulations -export ocean_simulation +export ocean_simulation, + nonhydrostatic_ocean_simulation using Oceananigans using Oceananigans.Units @@ -42,5 +43,6 @@ default_or_override(override, alternative_default=nothing) = override include("barotropic_potential_forcing.jl") include("radiative_forcing.jl") include("ocean_simulation.jl") +include("nonhydrostatic_ocean_simulation.jl") end # module diff --git a/src/OceanSimulations/nonhydrostatic_ocean_simulation.jl b/src/OceanSimulations/nonhydrostatic_ocean_simulation.jl new file mode 100644 index 000000000..e3072715f --- /dev/null +++ b/src/OceanSimulations/nonhydrostatic_ocean_simulation.jl @@ -0,0 +1,55 @@ + +""" + nonhydrostatic_ocean_simulation(grid; kwargs...) + +Build a NonhydrostaticModel-based ocean simulation configured for high-resolution LES. +The default configuration has WENO advection, TEOS-10 equation of state, and no turbulence closure. +""" +function nonhydrostatic_ocean_simulation(grid; Δt=1, + tracers = (:T, :S), + gravitational_acceleration = Oceananigans.defaults.gravitational_acceleration, + coriolis = nothing, + closure = nothing, + forcing = NamedTuple(), + biogeochemistry = nothing, + advection = WENO(order=9), + equation_of_state = TEOS10EquationOfState(), + boundary_conditions::NamedTuple = NamedTuple(), + verbose = false) + + FT = eltype(grid) + + τx = Field{Face, Center, Nothing}(grid) + τy = Field{Center, Face, Nothing}(grid) + Jᵀ = Field{Center, Center, Nothing}(grid) + Jˢ = Field{Center, Center, Nothing}(grid) + + u_top_bc = FluxBoundaryCondition(τx) + v_top_bc = FluxBoundaryCondition(τy) + T_top_bc = FluxBoundaryCondition(Jᵀ) + S_top_bc = FluxBoundaryCondition(Jˢ) + + default_boundary_conditions = (u = FieldBoundaryConditions(top=u_top_bc), + v = FieldBoundaryConditions(top=v_top_bc), + T = FieldBoundaryConditions(top=T_top_bc), + S = FieldBoundaryConditions(top=S_top_bc)) + + boundary_conditions = merge(default_boundary_conditions, boundary_conditions) + buoyancy = SeawaterBuoyancy(; gravitational_acceleration, equation_of_state) + + ocean_model = NonhydrostaticModel(; grid, + buoyancy, + closure, + biogeochemistry, + advection, + tracers, + coriolis, + forcing, + boundary_conditions) + + ocean = Simulation(ocean_model; Δt, verbose) + conjure_time_step_wizard!(ocean, cfl=0.7) + + return ocean +end +