From c8dc98412138fff1994af28fda141ad96e8ff5c2 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Wed, 11 Dec 2024 12:55:16 -0700 Subject: [PATCH] Refactor PrescribedAtmosphere and add an idealized atmosphere simulation (#251) * Refactor PrescribedAtmosphere and add an idealized atmosphere simulation * Add arctic night case * Convert to JRA55PrescribedAtmosphere * Import thermodynamics parameters * bugfixes * Fix bug with Default coriolis * Fix syntax * Update src/OceanSeaIceModels/PrescribedAtmospheres.jl Co-authored-by: Simone Silvestri --------- Co-authored-by: Navid C. Constantinou Co-authored-by: Simone Silvestri --- README.md | 2 +- docs/src/index.md | 2 +- examples/arctic_night.jl | 118 ++++++++++++++++++ examples/generate_surface_fluxes.jl | 2 +- .../idealized_single_column_simulation.jl | 118 ++++++++++++++++++ examples/inspect_JRA55_data.jl | 1 + examples/near_global_ocean_simulation.jl | 2 +- examples/single_column_os_papa_simulation.jl | 8 +- .../lat_lon_near_global_simulation.jl | 2 +- .../prototype_omip_simulation.jl | 4 +- .../one_degree_simulation.jl | 2 +- src/ClimaOcean.jl | 7 +- src/DataWrangling/JRA55.jl | 67 +++++----- .../atmosphere_ocean_fluxes.jl | 9 +- .../latitude_dependent_albedo.jl | 1 + .../CrossRealmFluxes/radiation.jl | 7 +- src/OceanSeaIceModels/OceanSeaIceModels.jl | 1 + .../PrescribedAtmospheres.jl | 98 ++++++++++----- src/OceanSimulations/OceanSimulations.jl | 25 ++-- test/runtests.jl | 2 +- test/test_jra55.jl | 7 +- ...est_ocean_sea_ice_model_parameter_space.jl | 2 +- test/test_simulations.jl | 2 +- test/test_surface_fluxes.jl | 25 ++-- 24 files changed, 401 insertions(+), 113 deletions(-) create mode 100644 examples/arctic_night.jl create mode 100644 examples/idealized_single_column_simulation.jl diff --git a/README.md b/README.md index 3f40c8fe..0836f75a 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,7 @@ set!(ocean.model, T = ClimaOcean.ECCOMetadata(:temperature; date), S = ClimaOcean.ECCOMetadata(:salinity; date)) # Build and run an OceanSeaIceModel (with no sea ice component) forced by JRA55 reanalysis -atmosphere = ClimaOcean.JRA55_prescribed_atmosphere(arch) +atmosphere = ClimaOcean.JRA55PrescribedAtmosphere(arch) coupled_model = ClimaOcean.OceanSeaIceModel(ocean; atmosphere) simulation = Simulation(coupled_model, Δt=5minutes, stop_time=30days) run!(simulation) diff --git a/docs/src/index.md b/docs/src/index.md index 0992cc60..7ccf1f70 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -54,7 +54,7 @@ set!(ocean.model, T = ClimaOcean.ECCOMetadata(:temperature; date), S = ClimaOcean.ECCOMetadata(:salinity; date)) # Build and run an OceanSeaIceModel (with no sea ice component) forced by JRA55 reanalysis -atmosphere = ClimaOcean.JRA55_prescribed_atmosphere(arch) +atmosphere = ClimaOcean.JRA55PrescribedAtmosphere(arch) coupled_model = ClimaOcean.OceanSeaIceModel(ocean; atmosphere) simulation = Simulation(coupled_model, Δt=5minutes, stop_time=30days) run!(simulation) diff --git a/examples/arctic_night.jl b/examples/arctic_night.jl new file mode 100644 index 00000000..caeb34f6 --- /dev/null +++ b/examples/arctic_night.jl @@ -0,0 +1,118 @@ +using ClimaOcean +using Oceananigans +using Oceananigans.Units +using SeawaterPolynomials + +# Ocean state parameters +T₀ = 0 # Surface temperature, ᵒC +S₀ = 35 # Surface salinity +N² = 1e-5 # Buoyancy gradient due to temperature stratification +f = 0 # Coriolis parameter + +# Atmospheric state parameters +Tₐ = 273.15 - 5 # Kelvin +u₁₀ = 10 # wind at 10 m, m/s +qₐ = 0.0 # specific humidity +Qℓ = 0 # shortwave radiation (W m⁻², positive means heating right now) + +# Build the atmosphere +atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat)) +atmosphere_times = range(0, 1day, length=3) +atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times) + +parent(atmosphere.tracers.T) .= Tₐ # K +parent(atmosphere.velocities.u) .= u₁₀ # m/s +parent(atmosphere.tracers.q) .= qₐ # mass ratio +parent(atmosphere.downwelling_radiation.shortwave) .= Qℓ # W + +# Build ocean model at rest with initial temperature stratification +grid = RectilinearGrid(size=400, z=(-400, 0), topology=(Flat, Flat, Bounded)) +ocean = ocean_simulation(grid, coriolis=FPlane(; f)) + +eos = ocean.model.buoyancy.model.equation_of_state +g = ocean.model.buoyancy.model.gravitational_acceleration +α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos) +dTdz = N² / (α * g) +Tᵢ(z) = max(-1.8, T₀ + dTdz * z) +set!(ocean.model, T=Tᵢ, S=S₀) + +radiation = Radiation(ocean_albedo=0.1) +model = OceanSeaIceModel(ocean; atmosphere, radiation) +simulation = Simulation(model, Δt=10.0, stop_time=3days) + +Q = model.fluxes.total.ocean.heat +Ql = model.fluxes.turbulent.fields.latent_heat +Qs = model.fluxes.turbulent.fields.sensible_heat +τx = model.fluxes.turbulent.fields.x_momentum +τy = model.fluxes.turbulent.fields.y_momentum +Fv = model.fluxes.turbulent.fields.water_vapor + +fluxes = (; Q, Ql, Qs, τx, τy, Fv) +outputs = merge(ocean.model.velocities, ocean.model.tracers, fluxes) +ow = JLD2OutputWriter(ocean.model, outputs, + schedule = TimeInterval(10minutes), + filename = "idealized_atmosphere.jld2", + overwrite_existing = true) + +simulation.output_writers[:ow] = ow + +run!(simulation) + +ut = FieldTimeSeries("idealized_atmosphere.jld2", "u") +vt = FieldTimeSeries("idealized_atmosphere.jld2", "v") +Tt = FieldTimeSeries("idealized_atmosphere.jld2", "T") +St = FieldTimeSeries("idealized_atmosphere.jld2", "S") +τx = FieldTimeSeries("idealized_atmosphere.jld2", "τx") +Q = FieldTimeSeries("idealized_atmosphere.jld2", "Q") +Nt = length(St) + +using GLMakie + +fig = Figure(size=(1000, 800)) + +axτ = Axis(fig[1, 1:3], xlabel="Time (days)", ylabel="Zonal momentum \n flux (N m⁻²)", xaxisposition=:top) +axQ = Axis(fig[2, 1:3], xlabel="Time (days)", ylabel="Heat flux \n (W m⁻²)") +axT = Axis(fig[3, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)") +axS = Axis(fig[3, 2], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)") +axu = Axis(fig[3, 3], xlabel="Velocities (m s⁻¹)", ylabel="z (m)", yaxisposition=:right) + +slider = Slider(fig[4, 1:3], startvalue=1, range=1:Nt) +n = slider.value + +un = @lift ut[$n] +vn = @lift vt[$n] +Tn = @lift Tt[$n] +Sn = @lift St[$n] + +lines!(axT, Tn) +lines!(axS, Sn) +lines!(axu, un, label="u") +lines!(axu, vn, label="v") +axislegend(axu) + +t = τx.times ./ days +tn = @lift t[$n] + +lines!(axτ, t, τx[:], label="τx") +vlines!(axτ, tn) + +lines!(axQ, t, Q[:]) +vlines!(axQ, tn) + +rowsize!(fig.layout, 1, Relative(0.2)) +rowsize!(fig.layout, 2, Relative(0.2)) + +xlims!(axS, 34.9, 35.1) +xlims!(axu, -0.1, 4.0) +hideydecorations!(axS, grid=false) +hidespines!(axT, :t, :r) +hidespines!(axS, :t, :l, :r) +hidespines!(axu, :t, :l) +hidespines!(axτ, :b, :r) +hidespines!(axQ, :t, :r) + +record(fig, "idealized_atmosphere.mp4", 1:Nt, framerate=24) do nn + @info "Drawing frame $nn of $Nt..." + n[] = nn +end + diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 946a8be5..4cfd35dd 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -44,7 +44,7 @@ save("ECCO_continents.png", fig) #hide # We invoke the constructor with only the first two time indices, corresponding to # January 1st (at 00:00 AM and 03:00 AM). -atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory()) +atmosphere = JRA55PrescribedAtmosphere(1:2; backend = InMemory()) ocean = ocean_simulation(grid) # Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity diff --git a/examples/idealized_single_column_simulation.jl b/examples/idealized_single_column_simulation.jl new file mode 100644 index 00000000..3d11a552 --- /dev/null +++ b/examples/idealized_single_column_simulation.jl @@ -0,0 +1,118 @@ +using ClimaOcean +using Oceananigans +using Oceananigans.Units +using SeawaterPolynomials + +# Ocean state parameters +T₀ = 20 # Surface temperature, ᵒC +S₀ = 35 # Surface salinity +N² = 1e-4 # Buoyancy gradient due to temperature stratification +f = 0 # Coriolis parameter + +# Atmospheric state parameters +Tₐ = 273.15 + 20 # Kelvin +u₁₀ = 10 # wind at 10 m, m/s +qₐ = 0.01 # specific humidity +Qℓ = 400 # shortwave radiation (W m⁻², positive means heating right now) + +# Build the atmosphere +atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat)) +atmosphere_times = range(0, 1day, length=3) +atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times) + +parent(atmosphere.tracers.T) .= Tₐ # K +parent(atmosphere.velocities.u) .= u₁₀ # m/s +parent(atmosphere.tracers.q) .= qₐ # mass ratio +parent(atmosphere.downwelling_radiation.shortwave) .= Qℓ # W + +# Build ocean model at rest with initial temperature stratification +grid = RectilinearGrid(size=100, z=(-400, 0), topology=(Flat, Flat, Bounded)) +ocean = ocean_simulation(grid, coriolis=FPlane(; f)) + +eos = ocean.model.buoyancy.model.equation_of_state +g = ocean.model.buoyancy.model.gravitational_acceleration +α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos) +dTdz = N² / (α * g) +Tᵢ(z) = T₀ + dTdz * z +set!(ocean.model, T=Tᵢ, S=S₀) + +radiation = Radiation(ocean_albedo=0.1) +model = OceanSeaIceModel(ocean; atmosphere, radiation) +simulation = Simulation(model, Δt=5minutes, stop_time=30days) + +Q = model.fluxes.total.ocean.heat +Ql = model.fluxes.turbulent.fields.latent_heat +Qs = model.fluxes.turbulent.fields.sensible_heat +τx = model.fluxes.turbulent.fields.x_momentum +τy = model.fluxes.turbulent.fields.y_momentum +Fv = model.fluxes.turbulent.fields.water_vapor + +fluxes = (; Q, Ql, Qs, τx, τy, Fv) +outputs = merge(ocean.model.velocities, ocean.model.tracers, fluxes) +ow = JLD2OutputWriter(ocean.model, outputs, + schedule = TimeInterval(1hour), + filename = "idealized_atmosphere.jld2", + overwrite_existing = true) + +simulation.output_writers[:ow] = ow + +run!(simulation) + +ut = FieldTimeSeries("idealized_atmosphere.jld2", "u") +vt = FieldTimeSeries("idealized_atmosphere.jld2", "v") +Tt = FieldTimeSeries("idealized_atmosphere.jld2", "T") +St = FieldTimeSeries("idealized_atmosphere.jld2", "S") +τx = FieldTimeSeries("idealized_atmosphere.jld2", "τx") +Q = FieldTimeSeries("idealized_atmosphere.jld2", "Q") +Nt = length(St) + +using GLMakie + +fig = Figure(size=(1000, 800)) + +axτ = Axis(fig[1, 1:3], xlabel="Time (days)", ylabel="Zonal momentum \n flux (N m⁻²)", xaxisposition=:top) +axQ = Axis(fig[2, 1:3], xlabel="Time (days)", ylabel="Heat flux \n (W m⁻²)") +axT = Axis(fig[3, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)") +axS = Axis(fig[3, 2], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)") +axu = Axis(fig[3, 3], xlabel="Velocities (m s⁻¹)", ylabel="z (m)", yaxisposition=:right) + +slider = Slider(fig[4, 1:3], startvalue=1, range=1:Nt) +n = slider.value + +un = @lift ut[$n] +vn = @lift vt[$n] +Tn = @lift Tt[$n] +Sn = @lift St[$n] + +lines!(axT, Tn) +lines!(axS, Sn) +lines!(axu, un, label="u") +lines!(axu, vn, label="v") +axislegend(axu) + +t = τx.times ./ days +tn = @lift t[$n] + +lines!(axτ, t, τx[:], label="τx") +vlines!(axτ, tn) + +lines!(axQ, t, Q[:]) +vlines!(axQ, tn) + +rowsize!(fig.layout, 1, Relative(0.2)) +rowsize!(fig.layout, 2, Relative(0.2)) + +xlims!(axS, 34.9, 35.1) +xlims!(axu, -0.1, 4.0) +hideydecorations!(axS, grid=false) +hidespines!(axT, :t, :r) +hidespines!(axS, :t, :l, :r) +hidespines!(axu, :t, :l) +hidespines!(axτ, :b, :r) +hidespines!(axQ, :t, :r) + +record(fig, "idealized_atmosphere.mp4", 1:Nt, framerate=24) do nn + @info "Drawing frame $nn of $Nt..." + n[] = nn +end + diff --git a/examples/inspect_JRA55_data.jl b/examples/inspect_JRA55_data.jl index 23445799..708bc1ca 100644 --- a/examples/inspect_JRA55_data.jl +++ b/examples/inspect_JRA55_data.jl @@ -43,6 +43,7 @@ label = @lift string("Repeat-year JRA55 forcing on year-day ", Label(fig[0, 1:2], label, fontsize=24) sf = surface!(axsw, x, y, z, color=Qswn, colorrange=(0, 1200)) + Colorbar(fig[2, 1], sf, vertical = false, width = Relative(0.5), diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl index f72a1985..6db964bc 100644 --- a/examples/near_global_ocean_simulation.jl +++ b/examples/near_global_ocean_simulation.jl @@ -107,7 +107,7 @@ radiation = Radiation(arch) # The number of snapshots that are loaded into memory is determined by # the `backend`. Here, we load 41 snapshots at a time into memory. -atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(41)) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41)) # ## The coupled simulation diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl index 11e5b585..ed6052d4 100644 --- a/examples/single_column_os_papa_simulation.jl +++ b/examples/single_column_os_papa_simulation.jl @@ -59,10 +59,10 @@ set!(ocean.model, T=ECCOMetadata(:temperature), S=ECCOMetadata(:salinity)) simulation_days = 31 snapshots_per_day = 8 # corresponding to JRA55's 3-hour frequency last_time = simulation_days * snapshots_per_day -atmosphere = JRA55_prescribed_atmosphere(1:last_time; - longitude = λ★, - latitude = φ★, - backend = InMemory()) +atmosphere = JRA55PrescribedAtmosphere(1:last_time; + longitude = λ★, + latitude = φ★, + backend = InMemory()) # This builds a representation of the atmosphere on the small grid diff --git a/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl b/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl index 97a2144a..125b7631 100644 --- a/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl +++ b/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl @@ -116,7 +116,7 @@ ocean.output_writers[:checkpoint] = Checkpointer(ocean.model, ##### backend = JRA55NetCDFBackend(4) -atmosphere = JRA55_prescribed_atmosphere(arch; backend) +atmosphere = JRA55PrescribedAtmosphere(arch; backend) radiation = Radiation(arch, ocean_albedo=LatitudeDependentAlbedo()) coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) diff --git a/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl b/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl index 1e6cd9ac..1aa5a404 100644 --- a/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl +++ b/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl @@ -15,7 +15,7 @@ using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: Radiation, SimilarityTheory using ClimaOcean.VerticalGrids: exponential_z_faces using ClimaOcean.JRA55 using ClimaOcean.ECCO -using ClimaOcean.JRA55: JRA55NetCDFBackend, JRA55_prescribed_atmosphere +using ClimaOcean.JRA55: JRA55NetCDFBackend, JRA55PrescribedAtmosphere using ClimaOcean.ECCO: ECCO_restoring_forcing, ECCO4Monthly, ECCO2Daily, ECCOMetadata using ClimaOcean.Bathymetry using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: LatitudeDependentAlbedo @@ -119,7 +119,7 @@ initial_date = dates[1] ##### backend = JRA55NetCDFBackend(4) -atmosphere = JRA55_prescribed_atmosphere(arch; backend) +atmosphere = JRA55PrescribedAtmosphere(arch; backend) radiation = Radiation(arch) coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 08577fed..27b6818f 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -88,7 +88,7 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), ##### radiation = Radiation(arch) -atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20)) ##### ##### Coupled simulation diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index b2c4f261..14b78efb 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -15,7 +15,8 @@ export SimilarityTheoryTurbulentFluxes, SkinTemperature, BulkTemperature, - JRA55_prescribed_atmosphere, + PrescribedAtmosphere, + JRA55PrescribedAtmosphere, JRA55NetCDFBackend, regrid_bathymetry, retrieve_bathymetry, @@ -87,7 +88,9 @@ using .InitialConditions using .OceanSeaIceModels using .OceanSimulations using .DataWrangling: JRA55, ECCO -using ClimaOcean.DataWrangling.JRA55: JRA55_prescribed_atmosphere, JRA55NetCDFBackend + +using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere +using ClimaOcean.DataWrangling.JRA55: JRA55PrescribedAtmosphere, JRA55NetCDFBackend using ClimaOcean.DataWrangling.ECCO end # module diff --git a/src/DataWrangling/JRA55.jl b/src/DataWrangling/JRA55.jl index d0a7b4b1..96ee8f9e 100644 --- a/src/DataWrangling/JRA55.jl +++ b/src/DataWrangling/JRA55.jl @@ -17,6 +17,7 @@ using ClimaOcean.DataWrangling: download_progress using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere, + PrescribedAtmosphereThermodynamicsParameters, TwoBandDownwellingRadiation using CUDA: @allowscalar @@ -36,6 +37,9 @@ function __init__() global download_jra55_cache = @get_scratch!("JRA55") end +struct JRA55Data end +const JRA55PrescribedAtmosphere = PrescribedAtmosphere{<:Any, <:JRA55Data} + # A list of all variables provided in the JRA55 dataset: JRA55_variable_names = (:river_freshwater_flux, :rain_freshwater_flux, @@ -626,15 +630,15 @@ end const AA = Oceananigans.Architectures.AbstractArchitecture -JRA55_prescribed_atmosphere(time_indices=Colon(); kw...) = - JRA55_prescribed_atmosphere(CPU(), time_indices; kw...) +JRA55PrescribedAtmosphere(time_indices=Colon(); kw...) = + JRA55PrescribedAtmosphere(CPU(), time_indices; kw...) -JRA55_prescribed_atmosphere(arch::Distributed, time_indices=Colon(); kw...) = - JRA55_prescribed_atmosphere(child_architecture(arch), time_indices; kw...) +JRA55PrescribedAtmosphere(arch::Distributed, time_indices=Colon(); kw...) = + JRA55PrescribedAtmosphere(child_architecture(arch), time_indices; kw...) # TODO: allow the user to pass dates """ - JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon(); + JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon(); backend = nothing, time_indexing = Cyclical(), reference_height = 10, # meters @@ -643,7 +647,7 @@ JRA55_prescribed_atmosphere(arch::Distributed, time_indices=Colon(); kw...) = Return a `PrescribedAtmosphere` representing JRA55 reanalysis data. """ -function JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon(); +function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon(); backend = nothing, time_indexing = Cyclical(), reference_height = 10, # meters @@ -675,45 +679,44 @@ function JRA55_prescribed_atmosphere(architecture::AA, time_indices=Colon(); Ql = JRA55_field_time_series(:downwelling_longwave_radiation; kw...) Qs = JRA55_field_time_series(:downwelling_shortwave_radiation; kw...) - freshwater_flux = (rain = Fra, - snow = Fsn) - - # Remember that rivers and icebergs are on a different grid and have - # a different frequency than the rest of the JRA55 data. We use `PrescribedAtmospheres` - # "auxiliary_freshwater_flux" feature to represent them. + # In JRA55, rivers and icebergs are on a different grid and stored with + # a different frequency than the rest of the data. Here, we use the + # `PrescribedAtmosphere.auxiliary_freshwater_flux` feature to represent them. if include_rivers_and_icebergs Fri = JRA55_field_time_series(:river_freshwater_flux; kw...) Fic = JRA55_field_time_series(:iceberg_freshwater_flux; kw...) - auxiliary_freshwater_flux = (rivers = Fri, icebergs = Fic) + auxiliary_freshwater_flux = (rivers=Fri, icebergs=Fic) else auxiliary_freshwater_flux = nothing end times = ua.times - - velocities = (u = ua, - v = va) - - tracers = (T = Ta, - q = qa) - + freshwater_flux = (rain=Fra, snow=Fsn) + velocities = (u=ua, v=va) + tracers = (T=Ta, q=qa) pressure = pa - downwelling_radiation = TwoBandDownwellingRadiation(shortwave=Qs, longwave=Ql) FT = eltype(ua) + boundary_layer_height = convert(FT, 600) reference_height = convert(FT, reference_height) - - atmosphere = PrescribedAtmosphere(times, FT; - velocities, - freshwater_flux, - auxiliary_freshwater_flux, - tracers, - downwelling_radiation, - reference_height, - pressure) - - return atmosphere + thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT) + grid = ua.grid + metadata = JRA55Data() + + return PrescribedAtmosphere(grid, + metadata, + velocities, + pressure, + tracers, + freshwater_flux, + auxiliary_freshwater_flux, + downwelling_radiation, + thermodynamics_parameters, + times, + reference_height, + boundary_layer_height) end end # module + diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index b80ad76b..1d8a0f65 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -35,13 +35,11 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) atmosphere_tracers = (T = atmosphere.tracers.T.data, q = atmosphere.tracers.q.data) - atmosphere_pressure = atmosphere.pressure.data - Qs = atmosphere.downwelling_radiation.shortwave Qℓ = atmosphere.downwelling_radiation.longwave downwelling_radiation = (shortwave=Qs.data, longwave=Qℓ.data) - freshwater_flux = map(ϕ -> ϕ.data, atmosphere.freshwater_flux) + atmosphere_pressure = atmosphere.pressure.data # Extract info for time-interpolation u = atmosphere.velocities.u # for example @@ -49,13 +47,13 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) atmosphere_backend = u.backend atmosphere_time_indexing = u.time_indexing - # kernel parameters that compute fluxes in 0:Nx+1 and 0:Ny+1 Nx, Ny, Nz = size(grid) single_column_grid = Nx == 1 && Ny == 1 if single_column_grid kernel_parameters = KernelParameters(1:1, 1:1) else + # Compute fluxes into halo regions, ie from 0:Nx+1 and 0:Ny+1 kernel_parameters = KernelParameters(0:Nx+1, 0:Ny+1) end @@ -97,7 +95,7 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) # TODO: do not assume that `auxiliary_freshater_flux` is a tuple auxiliary_data = map(ϕ -> ϕ.data, auxiliary_freshwater_flux) - first_auxiliary_flux = first(auxiliary_freshwater_flux) + first_auxiliary_flux = first(auxiliary_freshwater_flux) auxiliary_grid = first_auxiliary_flux.grid auxiliary_times = first_auxiliary_flux.times auxiliary_backend = first_auxiliary_flux.backend @@ -469,3 +467,4 @@ end # Note: positive implies _upward_ heat flux, and therefore cooling. return ϵ * σ * Tₒ^4 end + diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl b/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl index 26cee4fe..25b0edba 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl @@ -51,3 +51,4 @@ Base.show(io::IO, α::LatitudeDependentAlbedo) = print(io, summary(α)) α₁ = α.direct return α₀ - α₁ * hack_cosd(2φ) end + diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl index 54706edc..19ecf9d0 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl @@ -9,10 +9,9 @@ struct Radiation{FT, E, R} stefan_boltzmann_constant :: FT end -Adapt.adapt_structure(to, r :: Radiation) = - Radiation(Adapt.adapt(to, r.emission), - Adapt.adapt(to, r.reflection), - Adapt.adapt(to, r.stefan_boltzmann_constant)) +Adapt.adapt_structure(to, r :: Radiation) = Radiation(Adapt.adapt(to, r.emission), + Adapt.adapt(to, r.reflection), + Adapt.adapt(to, r.stefan_boltzmann_constant)) """ Radiation([arch = CPU(), FT=Float64]; diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 85e69d94..779a1f49 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -43,6 +43,7 @@ include("PrescribedAtmospheres.jl") using .PrescribedAtmospheres: PrescribedAtmosphere, + PrescribedAtmosphereThermodynamicsParameters, TwoBandDownwellingRadiation include("CrossRealmFluxes/CrossRealmFluxes.jl") diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/OceanSeaIceModels/PrescribedAtmospheres.jl index 63dcc85f..79a8fafd 100644 --- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl +++ b/src/OceanSeaIceModels/PrescribedAtmospheres.jl @@ -2,7 +2,8 @@ module PrescribedAtmospheres using Oceananigans.Grids: grid_name using Oceananigans.Utils: prettysummary -using Oceananigans.OutputReaders: update_field_time_series!, extract_field_time_series +using Oceananigans.Fields: Center +using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series using Adapt using Thermodynamics.Parameters: AbstractThermodynamicsParameters @@ -95,12 +96,12 @@ end const CP{FT} = ConstitutiveParameters{FT} where FT -@inline gas_constant(p::CP) = p.gas_constant -@inline molmass_dryair(p::CP) = p.dry_air_molar_mass -@inline molmass_water(p::CP) = p.water_molar_mass -@inline molmass_ratio(p::CP) = molmass_dryair(p) / molmass_water(p) -@inline R_v(p::CP) = gas_constant(p) / molmass_water(p) -@inline R_d(p::CP) = gas_constant(p) / molmass_dryair(p) +@inline gas_constant(p::CP) = p.gas_constant +@inline molmass_dryair(p::CP) = p.dry_air_molar_mass +@inline molmass_water(p::CP) = p.water_molar_mass +@inline molmass_ratio(p::CP) = molmass_dryair(p) / molmass_water(p) +@inline R_v(p::CP) = gas_constant(p) / molmass_water(p) +@inline R_d(p::CP) = gas_constant(p) / molmass_dryair(p) struct HeatCapacityParameters{FT} <: AbstractThermodynamicsParameters{FT} dry_air_adiabatic_exponent :: FT @@ -286,8 +287,9 @@ const PATP = PrescribedAtmosphereThermodynamicsParameters ##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic) ##### -struct PrescribedAtmosphere{FT, G, U, P, C, F, I, R, TP, TI} +struct PrescribedAtmosphere{FT, M, G, U, P, C, F, I, R, TP, TI} grid :: G + metadata :: M velocities :: U pressure :: P tracers :: C @@ -314,36 +316,72 @@ function Base.show(io::IO, pa::PrescribedAtmosphere) print(io, "└── boundary_layer_height: ", prettysummary(pa.boundary_layer_height)) end +function default_atmosphere_velocities(grid, times) + ua = FieldTimeSeries{Center, Center, Nothing}(grid, times) + va = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return (u=ua, v=va) +end + +function default_atmosphere_tracers(grid, times) + Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times) + qa = FieldTimeSeries{Center, Center, Nothing}(grid, times) + parent(Ta) .= 273.15 + 20 + return (T=Ta, q=qa) +end + +function default_downwelling_radiation(grid, times) + Qℓ = FieldTimeSeries{Center, Center, Nothing}(grid, times) + Qs = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return TwoBandDownwellingRadiation(shortwave=Qs, longwave=Qℓ) +end + +function default_freshwater_flux(grid, times) + rain = FieldTimeSeries{Center, Center, Nothing}(grid, times) + snow = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return (; rain, snow) +end + +function default_atmosphere_pressure(grid, times) + pa = FieldTimeSeries{Center, Center, Nothing}(grid, times) + parent(pa) .= 101325 + return pa +end + """ - PrescribedAtmosphere(times; - reference_height, - velocities = nothing, - pressure = nothing, - freshwater_flux = nothing, - downwelling_radiation = nothing, - tracers = nothing) + PrescribedAtmosphere(grid, times; + metadata = nothing, + reference_height = 10, # meters + boundary_layer_height = 600 # meters, + thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT), + auxiliary_freshwater_flux = nothing, + velocities = default_atmosphere_velocities(grid, times), + tracers = default_atmosphere_tracers(grid, times), + pressure = default_atmosphere_pressure(grid, times), + freshwater_flux = default_freshwater_flux(grid, times), + downwelling_radiation = default_downwelling_radiation(grid, times)) Return a representation of a prescribed time-evolving atmospheric state with data given at `times`. """ -function PrescribedAtmosphere(times, FT=Float64; - reference_height, - velocities = nothing, - boundary_layer_height = convert(FT, 600), - pressure = nothing, - freshwater_flux = nothing, +function PrescribedAtmosphere(grid, times; + metadata = nothing, + reference_height = convert(eltype(grid), 10), + boundary_layer_height = convert(eltype(grid), 600), + thermodynamics_parameters = nothing, auxiliary_freshwater_flux = nothing, - downwelling_radiation = nothing, - thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT), - grid = nothing, - tracers = nothing) - - if isnothing(grid) # try to find it - u = first(velocities) - grid = u.grid + velocities = default_atmosphere_velocities(grid, times), + tracers = default_atmosphere_tracers(grid, times), + pressure = default_atmosphere_pressure(grid, times), + freshwater_flux = default_freshwater_flux(grid, times), + downwelling_radiation = default_downwelling_radiation(grid, times)) + + FT = eltype(grid) + if isnothing(thermodynamics_parameters) + thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT) end return PrescribedAtmosphere(grid, + metadata, velocities, pressure, tracers, @@ -384,7 +422,7 @@ TwoBandDownwellingRadiation(; shortwave=nothing, longwave=nothing) = Adapt.adapt_structure(to, tsdr::TwoBandDownwellingRadiation) = TwoBandDownwellingRadiation(adapt(to, tsdr.shortwave), - adapt(to, tsdr.longwave)) + adapt(to, tsdr.longwave)) end # module diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 8c95965b..8040ed07 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -6,6 +6,7 @@ using Oceananigans using Oceananigans.Units using Oceananigans.Utils: with_tracers using Oceananigans.Advection: FluxFormAdvection +using Oceananigans.BoundaryConditions: DefaultBoundaryCondition using Oceananigans.Coriolis: ActiveCellEnstrophyConserving using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node using OrthogonalSphericalShellGrids @@ -36,7 +37,7 @@ For example, the default bottom drag should be 0 for a single column model, but We therefore need a way to specify both the "normal" default 0.003 as well as the "alternative default" 0, all while respecting user input and changing this to a new value if specified. """ -default_or_override(default::Default, possibly_alternative_default=default.value) = possibly_alternative_default +default_or_override(default::Default, possibly_alternative_default=default.value) = possibly_alternative_default default_or_override(override, alternative_default=nothing) = override # Some defaults @@ -84,7 +85,7 @@ function ocean_simulation(grid; gravitational_acceleration = g_Earth, bottom_drag_coefficient = Default(0.003), forcing = NamedTuple(), - coriolis = HydrostaticSphericalCoriolis(; rotation_rate), + coriolis = Default(HydrostaticSphericalCoriolis(; rotation_rate)), momentum_advection = default_momentum_advection(), equation_of_state = TEOS10EquationOfState(; reference_density), tracer_advection = default_tracer_advection(), @@ -92,6 +93,12 @@ function ocean_simulation(grid; FT = eltype(grid) + if grid isa RectilinearGrid # turn off Coriolis unless user-supplied + coriolis = default_or_override(coriolis, nothing) + else + coriolis = default_or_override(coriolis) + end + # Detect whether we are on a single column grid Nx, Ny, _ = size(grid) single_column_simulation = Nx == 1 && Ny == 1 @@ -105,8 +112,8 @@ function ocean_simulation(grid; momentum_advection = nothing # No immersed boundaries in a single column grid - u_immersed_bc = nothing - v_immersed_bc = nothing + u_immersed_bc = DefaultBoundaryCondition() + v_immersed_bc = DefaultBoundaryCondition() else if !(grid isa ImmersedBoundaryGrid) msg = """Are you totally, 100% sure that you want to build a simulation on @@ -143,11 +150,11 @@ function ocean_simulation(grid; u_bot_bc = FluxBoundaryCondition(u_quadratic_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) v_bot_bc = FluxBoundaryCondition(v_quadratic_bottom_drag, discrete_form=true, parameters=bottom_drag_coefficient) - - ocean_boundary_conditions = (u = FieldBoundaryConditions(top = u_top_bc, bottom = u_bot_bc, immersed = u_immersed_bc), - v = FieldBoundaryConditions(top = v_top_bc, bottom = v_bot_bc, immersed = v_immersed_bc), - T = FieldBoundaryConditions(top = T_top_bc), - S = FieldBoundaryConditions(top = S_top_bc)) + + ocean_boundary_conditions = (u = FieldBoundaryConditions(top=u_top_bc, bottom=u_bot_bc, immersed=u_immersed_bc), + v = FieldBoundaryConditions(top=v_top_bc, bottom=v_bot_bc, immersed=v_immersed_bc), + T = FieldBoundaryConditions(top=T_top_bc), + S = FieldBoundaryConditions(top=S_top_bc)) # Use the TEOS10 equation of state teos10 = TEOS10EquationOfState(; reference_density) diff --git a/test/runtests.jl b/test/runtests.jl index 0c71ec17..932f48ef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,7 +23,7 @@ if test_group == :init || test_group == :all #### Download JRA55 data #### - atmosphere = JRA55_prescribed_atmosphere() + atmosphere = JRA55PrescribedAtmosphere() #### #### Download ECCO data diff --git a/test/test_jra55.jl b/test/test_jra55.jl index 47f0b2e7..3fc3e8b6 100644 --- a/test/test_jra55.jl +++ b/test/test_jra55.jl @@ -123,13 +123,12 @@ using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere ##### JRA55 prescribed atmosphere ##### - Nb = 4 - backend = JRA55NetCDFBackend(Nb) - atmosphere = JRA55_prescribed_atmosphere(arch; backend, include_rivers_and_icebergs=false) + backend = JRA55NetCDFBackend(2) + atmosphere = JRA55PrescribedAtmosphere(arch; backend, include_rivers_and_icebergs=false) @test atmosphere isa PrescribedAtmosphere @test isnothing(atmosphere.auxiliary_freshwater_flux) - atmosphere = JRA55_prescribed_atmosphere(arch; backend, include_rivers_and_icebergs=true) + atmosphere = JRA55PrescribedAtmosphere(arch; backend, include_rivers_and_icebergs=true) @test haskey(atmosphere.auxiliary_freshwater_flux, :rivers) @test haskey(atmosphere.auxiliary_freshwater_flux, :icebergs) end diff --git a/test/test_ocean_sea_ice_model_parameter_space.jl b/test/test_ocean_sea_ice_model_parameter_space.jl index ec4c2ccf..a1c9d9fa 100644 --- a/test/test_ocean_sea_ice_model_parameter_space.jl +++ b/test/test_ocean_sea_ice_model_parameter_space.jl @@ -30,7 +30,7 @@ model = ocean.model start_time = time_ns() backend = JRA55NetCDFBackend(4) -atmosphere = JRA55_prescribed_atmosphere(arch; backend) +atmosphere = JRA55PrescribedAtmosphere(arch; backend) radiation = Radiation(arch) elapsed = 1e-9 * (time_ns() - start_time) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 123199de..b089d489 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -25,7 +25,7 @@ using OrthogonalSphericalShellGrids ocean = ocean_simulation(grid; free_surface) backend = JRA55NetCDFBackend(4) - atmosphere = JRA55_prescribed_atmosphere(arch; backend) + atmosphere = JRA55PrescribedAtmosphere(arch; backend) radiation = Radiation(arch) # Fluxes are computed when the model is constructed, so we just test that this works. diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index 6aba08c7..ecbc337a 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -32,18 +32,19 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, @testset "Test surface fluxes" begin for arch in test_architectures grid = LatitudeLongitudeGrid(arch; - size = 1, - latitude = 0, - longitude = 0, - z = (-1, 0), - topology = (Flat, Flat, Bounded)) + size = 1, + latitude = 0, + longitude = 0, + z = (-1, 0), + topology = (Flat, Flat, Bounded)) - ocean = ocean_simulation(grid; momentum_advection = nothing, - tracer_advection = nothing, - closure = nothing, - bottom_drag_coefficient = 0.0) + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + bottom_drag_coefficient = 0.0) - atmosphere = JRA55_prescribed_atmosphere(1:2; grid, architecture = arch, backend = InMemory()) + atmosphere = JRA55PrescribedAtmosphere(1:2; grid, architecture = arch, backend = InMemory()) CUDA.@allowscalar begin h = atmosphere.reference_height @@ -177,7 +178,7 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, closure = nothing, bottom_drag_coefficient = 0.0) - atmosphere = JRA55_prescribed_atmosphere(1:2; grid, architecture = arch, backend = InMemory()) + atmosphere = JRA55PrescribedAtmosphere(1:2; grid, architecture = arch, backend = InMemory()) fill!(ocean.model.tracers.T, -2.0) @@ -235,7 +236,7 @@ end set!(ocean.model; T=T_metadata, S=S_metadata) - atmosphere = JRA55_prescribed_atmosphere(1:10; grid, architecture = arch, backend = InMemory()) + atmosphere = JRA55PrescribedAtmosphere(1:10; grid, architecture = arch, backend = InMemory()) radiation = Radiation(ocean_albedo=0.1, ocean_emissivity=1.0) sea_ice = nothing