diff --git a/.gitignore b/.gitignore index 4acc5052a..c553b3bb0 100644 --- a/.gitignore +++ b/.gitignore @@ -67,3 +67,4 @@ CondaPkg.toml # claude .claude +test/Manifest.toml diff --git a/docs/src/developers/slab_ocean.jl b/docs/src/developers/slab_ocean.jl index e05a1674d..9f29dd87d 100644 --- a/docs/src/developers/slab_ocean.jl +++ b/docs/src/developers/slab_ocean.jl @@ -141,7 +141,7 @@ set!(sea_ice.model, h=Metadatum(:sea_ice_thickness, dataset=ECCO4Monthly()), ℵ=Metadatum(:sea_ice_concentration, dataset=ECCO4Monthly())) interfaces = ComponentInterfaces(atmosphere, slab_ocean, sea_ice; exchange_grid=grid) -coupled_model = NumericalEarth.EarthSystemModel(atmosphere, slab_ocean, sea_ice; interfaces) +coupled_model = NumericalEarth.EarthSystemModel(; atmosphere, sea_ice, ocean=slab_ocean, interfaces) simulation = Simulation(coupled_model, Δt=60minutes, stop_time=120days) run!(simulation) diff --git a/docs/src/earth_system_model.md b/docs/src/earth_system_model.md index f50e40868..d7f7a278c 100644 --- a/docs/src/earth_system_model.md +++ b/docs/src/earth_system_model.md @@ -23,6 +23,7 @@ model = OceanOnlyModel(ocean) # output EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) +├── radiation: Nothing ├── atmosphere: Nothing ├── land: Nothing ├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} @@ -51,6 +52,7 @@ model # output EarthSystemModel{CPU}(time = 1 hour, iteration = 3) +├── radiation: Nothing ├── atmosphere: Nothing ├── land: Nothing ├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} @@ -87,10 +89,11 @@ a sea ice component as positional arguments: ```jldoctest esm ocean = ocean_simulation(grid, timestepper = :QuasiAdamsBashforth2) sea_ice = FreezingLimitedOceanTemperature() -model = OceanSeaIceModel(ocean, sea_ice) +model = OceanSeaIceModel(sea_ice, ocean) # output EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) +├── radiation: Nothing ├── atmosphere: Nothing ├── land: Nothing ├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} @@ -108,9 +111,12 @@ Oceananigans `Simulation`, would also work. EarthSystemModel ``` -The full constructor takes positional arguments `(atmosphere, ocean, sea_ice)` and -gives access to every knob: radiation parameters, reference densities, heat capacities, -and -- most importantly -- the `interfaces` keyword, which controls how fluxes are computed. +The full constructor takes positional arguments +`(radiation, atmosphere, land, sea_ice, ocean)` -- the components in struct +order, top to bottom -- and gives access to every knob: reference densities, +heat capacities, and -- most importantly -- the `interfaces` keyword, which +controls how fluxes are computed. Pass `nothing` for components that are +absent. ## Customizing flux formulations @@ -130,6 +136,7 @@ model = OceanOnlyModel(ocean; interfaces) # output EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) +├── radiation: Nothing ├── atmosphere: Nothing ├── land: Nothing ├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} @@ -148,6 +155,7 @@ model = OceanOnlyModel(ocean; interfaces) # output EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) +├── radiation: Nothing ├── atmosphere: Nothing ├── land: Nothing ├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 3043bef0f..1e09eab1a 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -62,11 +62,12 @@ S_metadata = ECCOMetadatum(:salinity; date=DateTime(1993, 1, 1)) set!(ocean.model; T=T_metadata, S=S_metadata) # Finally, we construct a coupled model, which will compute fluxes during construction. -# We omit `sea_ice` so the model is ocean-only, and use the default `Radiation()` that -# uses the two-band shortwave (visible and UV) + longwave (mid and far infrared) -# decomposition of the radiation spectrum. +# We omit `sea_ice` so the model is ocean-only, and pair the JRA55 atmosphere with a +# matching `JRA55PrescribedRadiation` that supplies downwelling shortwave and +# longwave radiation as well as ocean / sea-ice surface properties. -coupled_model = OceanOnlyModel(ocean; atmosphere, radiation=Radiation(eltype)) +radiation = JRA55PrescribedRadiation(; backend = JRA55NetCDFBackend(2)) +coupled_model = OceanOnlyModel(ocean; atmosphere, radiation) # Now that the surface fluxes are computed, we can extract and visualize them. # The turbulent fluxes are stored in `coupled_model.interfaces.atmosphere_ocean_interface.fluxes`. diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl index 6fbac8d07..1eac5f6e3 100644 --- a/examples/global_climate_simulation.jl +++ b/examples/global_climate_simulation.jl @@ -98,10 +98,12 @@ nothing #hide nothing #hide # We build the complete coupled `earth_model` and the coupled simulation. -# Since radiation is idealized in this example, we set the emissivities to zero. +# NumericalEarth still computes turbulent (sensible heat, water vapor) fluxes +# here using its own bulk-formula machinery and writes them back to Speedy. +# Top-level radiation, however, is not yet wired up against SpeedyWeather's +# upwelling-LW path, so we leave `radiation = nothing` for now. -radiation = Radiation(ocean_emissivity=0.0, sea_ice_emissivity=0.0) -earth_model = EarthSystemModel(atmosphere, ocean, sea_ice; radiation) +earth_model = EarthSystemModel(; atmosphere, sea_ice, ocean) # ## Building and running the simulation # diff --git a/examples/idealized_single_column_simulation.jl b/examples/idealized_single_column_simulation.jl index c6e4a498c..ec76bece1 100644 --- a/examples/idealized_single_column_simulation.jl +++ b/examples/idealized_single_column_simulation.jl @@ -18,15 +18,20 @@ qᵃᵗ = 0.01 # specific humidity ℐꜜˢʷ = 400 # shortwave radiation (W m⁻², positive means heating right now) # Build the atmosphere -radiation = Radiation(ocean_albedo=0.1) atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat)) atmosphere_times = range(0, 1days, length=3) atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times) +# Build the radiation component (lives on the same grid + times as the atmosphere +# in this single-column setup) and prescribe a constant shortwave forcing. +# Override the default ocean albedo (0.05) but keep the default emissivity (0.97). +radiation = PrescribedRadiation(atmosphere_grid, atmosphere_times; + ocean_surface = SurfaceRadiationProperties(albedo=0.1)) + 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) .= ℐꜜˢʷ # W +parent(radiation.downwelling_shortwave) .= ℐꜜˢʷ # W # Build ocean model at rest with initial temperature stratification grid = RectilinearGrid(size=20, z=(-100, 0), topology=(Flat, Flat, Bounded)) @@ -40,8 +45,8 @@ Tᵢ(z) = T₀ + dTdz * z set!(ocean.model, T=Tᵢ, S=S₀) atmosphere_ocean_fluxes = SimilarityTheoryFluxes(stability_functions=nothing) -interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes) -model = OceanOnlyModel(ocean; atmosphere, interfaces) +interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes, radiation) +model = OceanOnlyModel(ocean; atmosphere, radiation, interfaces) 𝒬ᵛ = model.interfaces.atmosphere_ocean_interface.fluxes.latent_heat 𝒬ᵀ = model.interfaces.atmosphere_ocean_interface.fluxes.sensible_heat diff --git a/examples/meridional_heat_transport_ecco.jl b/examples/meridional_heat_transport_ecco.jl index 58b3e7849..fd6ad9766 100755 --- a/examples/meridional_heat_transport_ecco.jl +++ b/examples/meridional_heat_transport_ecco.jl @@ -42,10 +42,10 @@ ecco_sea_ice_concentration = Metadatum(:sea_ice_concentration; date, dataset) set!(ocean.model, T=ecco_temperature, S=ecco_salinity) set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) -radiation = Radiation(arch) -atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80), - include_rivers_and_icebergs = false) -esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +jra55_backend = JRA55NetCDFBackend(80) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend) +radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend) +esm = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation) simulation = Simulation(esm; Δt=20minutes, stop_time=5*365days) diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl index f57a5e8e0..1ec9c3ce0 100644 --- a/examples/near_global_ocean_simulation.jl +++ b/examples/near_global_ocean_simulation.jl @@ -88,24 +88,14 @@ set!(ocean.model, T=Metadatum(:temperature, dataset=ECCO4Monthly()), # ### Prescribed atmosphere and radiation # -# Next we build a prescribed atmosphere state and radiation model, -# which will drive the ocean simulation. We use the default `Radiation` model, - -# The radiation model specifies an ocean albedo emissivity to compute the net radiative -# fluxes. The default ocean albedo is based on Payne (1982) and depends on cloud cover -# (calculated from the ratio of maximum possible incident solar radiation to actual -# incident solar radiation) and latitude. The ocean emissivity is set to 0.97. - -radiation = Radiation(arch) - -# The atmospheric data is prescribed using the JRA55 dataset. -# The JRA55 dataset provides atmospheric data such as temperature, humidity, and winds -# to calculate turbulent fluxes using bulk formulae, see [`InterfaceComputations`](@ref NumericalEarth.EarthSystemModels.InterfaceComputations). -# The number of snapshots that are loaded into memory is determined by -# the `backend`. Here, we load 41 snapshots at a time into memory. +# Next we build a prescribed atmosphere state and radiation component, +# which together drive the ocean simulation. The atmospheric data and +# downwelling shortwave / longwave radiation are both prescribed using JRA55. +# The number of snapshots loaded into memory is set by the backend. jra55_backend = JRA55NetCDFBackend(41) atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend) +radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend) land = JRA55PrescribedLand(arch; backend=jra55_backend) # ## The coupled simulation diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl index a08bba029..53e09aaf9 100644 --- a/examples/one_degree_simulation.jl +++ b/examples/one_degree_simulation.jl @@ -95,9 +95,13 @@ set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) # ### Atmospheric forcing # We force the simulation with a JRA55-do atmospheric reanalysis. -radiation = Radiation(arch) jra55_backend = JRA55NetCDFBackend(80) atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend) +# Use a latitude-dependent ocean albedo (Large & Yeager 2009); keep the +# default ocean emissivity (0.97) and sea-ice surface (albedo 0.7, +# emissivity 1.0). +radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend, + ocean_surface = SurfaceRadiationProperties(albedo = LatitudeDependentAlbedo())) land = JRA55PrescribedLand(arch; backend=jra55_backend) # ### Coupled simulation @@ -107,7 +111,7 @@ land = JRA55PrescribedLand(arch; backend=jra55_backend) # With Runge-Kutta 3rd order time-stepping we can safely use a timestep of 20 minutes. -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, land, radiation) +coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, land, radiation) simulation = Simulation(coupled_model; Δt=20minutes, stop_time=365days) # ### A progress messenger diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl index d0a38baf6..1420687fd 100644 --- a/examples/single_column_os_papa_simulation.jl +++ b/examples/single_column_os_papa_simulation.jl @@ -68,6 +68,11 @@ atmosphere = JRA55PrescribedAtmosphere(longitude = λ★, end_date = DateTime(1990, 1, 31), # Last day of the simulation backend = InMemory()) +radiation = JRA55PrescribedRadiation(longitude = λ★, + latitude = φ★, + end_date = DateTime(1990, 1, 31), + backend = InMemory()) + # This builds a representation of the atmosphere on the small grid atmosphere.grid @@ -101,7 +106,6 @@ lines!(axq, t_days, qa) current_figure() # We continue constructing a simulation. -radiation = Radiation() coupled_model = OceanOnlyModel(ocean; atmosphere, radiation) simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) @@ -205,8 +209,8 @@ ua = atmosphere.velocities.u va = atmosphere.velocities.v Ta = atmosphere.tracers.T qa = atmosphere.tracers.q -ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave -ℐꜜˢʷ = atmosphere.downwelling_radiation.shortwave +ℐꜜˡʷ = radiation.downwelling_longwave +ℐꜜˢʷ = radiation.downwelling_shortwave Pr = atmosphere.freshwater_flux.rain Ps = atmosphere.freshwater_flux.snow diff --git a/examples/veros_ocean_forced_simulation.jl b/examples/veros_ocean_forced_simulation.jl index bb95573f5..598582683 100644 --- a/examples/veros_ocean_forced_simulation.jl +++ b/examples/veros_ocean_forced_simulation.jl @@ -69,12 +69,11 @@ ocean.set_forcing = set_forcing_tke_only # radiation, as well as freshwater fluxes. atmos = JRA55PrescribedAtmosphere(; backend = JRA55NetCDFBackend(10)) +radiation = JRA55PrescribedRadiation(; backend = JRA55NetCDFBackend(10)) -# The coupled ocean--atmosphere model. -# We use the default radiation model and we do not couple an ice model for simplicity. +# The coupled ocean--atmosphere model. We do not couple an ice model for simplicity. -radiation = Radiation() -coupled_model = OceanSeaIceModel(ocean, nothing; atmosphere=atmos, radiation) +coupled_model = OceanSeaIceModel(nothing, ocean; atmosphere=atmos, radiation) simulation = Simulation(coupled_model; Δt = 1800, stop_time = 60days) # We set up a progress callback that will print the current time, iteration, and maximum velocities diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index afabdf826..eee1fdb48 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -90,13 +90,13 @@ set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset), ##### atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40)) -radiation = Radiation() +radiation = JRA55PrescribedRadiation(arch; backend=JRA55NetCDFBackend(40)) ##### ##### Arctic coupled model ##### -arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +arctic = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation) arctic = Simulation(arctic, Δt=5minutes, stop_time=365days) # Sea-ice variables diff --git a/experiments/coupled_simulation/earth_system_coupled_simulation.jl b/experiments/coupled_simulation/earth_system_coupled_simulation.jl index 334ebce91..73a81d5d5 100644 --- a/experiments/coupled_simulation/earth_system_coupled_simulation.jl +++ b/experiments/coupled_simulation/earth_system_coupled_simulation.jl @@ -79,14 +79,17 @@ salinity = ECCOMetadata(:salinity; dir="./") ice_thickness = ECCOMetadata(:sea_ice_thickness; dir="./") ice_concentration = ECCOMetadata(:sea_ice_concentration; dir="./") -atmosphere = JRA55PrescribedAtmosphere(arch, backend=JRA55NetCDFBackend(20)) -radiation = Radiation(ocean_albedo = LatitudeDependentAlbedo(), sea_ice_albedo=0.6) +jra55_backend = JRA55NetCDFBackend(20) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend) +radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend, + ocean_surface = SurfaceRadiationProperties(LatitudeDependentAlbedo(), 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.6, 1.0)) set!(ocean.model, T=temperature, S=salinity) set!(sea_ice.model.ice_thickness, ice_thickness, inpainting=NearestNeighborInpainting(1)) set!(sea_ice.model.ice_concentration, ice_concentration, inpainting=NearestNeighborInpainting(1)) -earth_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +earth_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation) earth = Simulation(earth_model; Δt=30minutes, stop_iteration=10, stop_time=30days) u, v, _ = ocean.model.velocities diff --git a/experiments/flux_climatology/flux_climatology.jl b/experiments/flux_climatology/flux_climatology.jl index 20ad23d60..dfda99ab8 100644 --- a/experiments/flux_climatology/flux_climatology.jl +++ b/experiments/flux_climatology/flux_climatology.jl @@ -197,7 +197,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(1000)) ##### A prescribed earth... ##### -earth_model = OceanOnlyModel(ocean; atmosphere, radiation = Radiation(ocean_emissivity=0, ocean_albedo=1)) +earth_model = OceanOnlyModel(ocean; atmosphere) # radiation off (radiation=nothing default) earth = Simulation(earth_model, Δt=3hours, stop_time=365days) wall_time = Ref(time_ns()) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index e48bc4d56..215beafd6 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -45,8 +45,9 @@ ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surfac set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), S=ECCOMetadata(:salinity; dates=first(dates))) -radiation = Radiation(arch) -atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41)) +jra55_backend = JRA55NetCDFBackend(41) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=jra55_backend) +radiation = JRA55PrescribedRadiation(arch; backend=jra55_backend) coupled_model = OceanOnlyModel(ocean; atmosphere, radiation) simulation = Simulation(coupled_model; Δt=10minutes, stop_iteration=100) diff --git a/ext/NumericalEarthReactantExt.jl b/ext/NumericalEarthReactantExt.jl index 3a5ccab54..93c3a6f06 100644 --- a/ext/NumericalEarthReactantExt.jl +++ b/ext/NumericalEarthReactantExt.jl @@ -13,9 +13,9 @@ const OceananigansReactantExt = Base.get_extension( Oceananigans, :OceananigansReactantExt ) -const ReactantOSIM{I, A, L, O, F, C} = Union{ - EarthSystemModel{I, A, L, O, F, C, <:ReactantState}, - EarthSystemModel{I, A, L, O, F, C, <:Distributed{ReactantState}}, +const ReactantOSIM{R, A, L, I, O, F, C} = Union{ + EarthSystemModel{R, A, L, I, O, F, C, <:ReactantState}, + EarthSystemModel{R, A, L, I, O, F, C, <:Distributed{ReactantState}}, } reconcile_state!(model::ReactantOSIM) = nothing diff --git a/src/Atmospheres/interpolate_atmospheric_state.jl b/src/Atmospheres/interpolate_atmospheric_state.jl index 56cf48727..2f0b365ec 100644 --- a/src/Atmospheres/interpolate_atmospheric_state.jl +++ b/src/Atmospheres/interpolate_atmospheric_state.jl @@ -27,9 +27,6 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c atmosphere_tracers = (T = atmosphere.tracers.T.data, q = atmosphere.tracers.q.data) - ℐꜜˢʷ = atmosphere.downwelling_radiation.shortwave - ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave - downwelling_radiation = (shortwave=ℐꜜˢʷ.data, longwave=ℐꜜˡʷ.data) freshwater_flux = map(ϕ -> ϕ.data, atmosphere.freshwater_flux) atmosphere_pressure = atmosphere.pressure.data @@ -44,14 +41,12 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c # Simplify NamedTuple to reduce parameter space consumption. # See https://github.com/CliMA/NumericalEarth.jl/issues/116. - atmosphere_data = (u = atmosphere_fields.u.data, - v = atmosphere_fields.v.data, - T = atmosphere_fields.T.data, - p = atmosphere_fields.p.data, - q = atmosphere_fields.q.data, - ℐꜜˢʷ = atmosphere_fields.ℐꜜˢʷ.data, - ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data, - Jᶜ = atmosphere_fields.Jᶜ.data) + atmosphere_data = (u = atmosphere_fields.u.data, + v = atmosphere_fields.v.data, + T = atmosphere_fields.T.data, + p = atmosphere_fields.p.data, + q = atmosphere_fields.q.data, + Jᶜ = atmosphere_fields.Jᶜ.data) kernel_parameters = interface_kernel_parameters(grid) @@ -72,7 +67,6 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c atmosphere_velocities, atmosphere_tracers, atmosphere_pressure, - downwelling_radiation, freshwater_flux, atmosphere_backend, atmosphere_time_indexing) @@ -99,7 +93,6 @@ end atmos_velocities, atmos_tracers, atmos_pressure, - downwelling_radiation, prescribed_freshwater_flux, atmos_backend, atmos_time_indexing) @@ -121,9 +114,6 @@ end qᵃᵗ = interp_atmos_time_series(atmos_tracers.q, atmos_args...) pᵃᵗ = interp_atmos_time_series(atmos_pressure, atmos_args...) - ℐꜜˢʷ = interp_atmos_time_series(downwelling_radiation.shortwave, atmos_args...) - ℐꜜˡʷ = interp_atmos_time_series(downwelling_radiation.longwave, atmos_args...) - # Usually precipitation Mh = interp_atmos_time_series(prescribed_freshwater_flux, atmos_args...) @@ -138,8 +128,6 @@ end surface_atmos_state.T[i, j, 1] = Tᵃᵗ surface_atmos_state.p[i, j, 1] = pᵃᵗ surface_atmos_state.q[i, j, 1] = qᵃᵗ - surface_atmos_state.ℐꜜˢʷ[i, j, 1] = ℐꜜˢʷ - surface_atmos_state.ℐꜜˡʷ[i, j, 1] = ℐꜜˡʷ surface_atmos_state.Jᶜ[i, j, 1] = Mh end end diff --git a/src/Atmospheres/prescribed_atmosphere.jl b/src/Atmospheres/prescribed_atmosphere.jl index b753a506e..bbb5158ca 100644 --- a/src/Atmospheres/prescribed_atmosphere.jl +++ b/src/Atmospheres/prescribed_atmosphere.jl @@ -2,14 +2,13 @@ ##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic) ##### -mutable struct PrescribedAtmosphere{FT, G, T, U, P, C, F, R, TP, TI} +mutable struct PrescribedAtmosphere{FT, G, T, U, P, C, F, TP, TI} grid :: G clock :: Clock{T} velocities :: U pressure :: P tracers :: C freshwater_flux :: F - downwelling_radiation :: R thermodynamics_parameters :: TP times :: TI surface_layer_height :: FT @@ -43,12 +42,6 @@ function default_atmosphere_tracers(grid, times) return (T=Ta, q=qa) end -function default_downwelling_radiation(grid, times) - ℐꜜˡʷ = FieldTimeSeries{Center, Center, Nothing}(grid, times) - ℐꜜˢʷ = FieldTimeSeries{Center, Center, Nothing}(grid, times) - return TwoBandDownwellingRadiation(shortwave=ℐꜜˢʷ, longwave=ℐꜜˡʷ) -end - function default_freshwater_flux(grid, times) rain = FieldTimeSeries{Center, Center, Nothing}(grid, times) snow = FieldTimeSeries{Center, Center, Nothing}(grid, times) @@ -93,27 +86,29 @@ update_net_fluxes!(coupled_model, ::PrescribedAtmosphere) = nothing PrescribedAtmosphere(grid, times=[zero(grid)]; clock = Clock{Float64}(time = 0), surface_layer_height = 10, # meters - boundary_layer_height = 512 # meters, + boundary_layer_height = 512, # meters thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(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)) + 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)) Return a representation of a prescribed time-evolving atmospheric state with data given at `times`. + +Note: downwelling shortwave / longwave radiation is now part of the +top-level `radiation` component (see `PrescribedRadiation`, +`JRA55PrescribedRadiation`). """ function PrescribedAtmosphere(grid, times=[zero(grid)]; clock = Clock{Float64}(time = 0), surface_layer_height = 10, boundary_layer_height = 512, thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(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)) + 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)) FT = eltype(grid) if isnothing(thermodynamics_parameters) @@ -126,7 +121,6 @@ function PrescribedAtmosphere(grid, times=[zero(grid)]; pressure, tracers, freshwater_flux, - downwelling_radiation, thermodynamics_parameters, times, convert(FT, surface_layer_height), @@ -136,36 +130,17 @@ function PrescribedAtmosphere(grid, times=[zero(grid)]; return atmosphere end -struct TwoBandDownwellingRadiation{SW, LW} - shortwave :: SW - longwave :: LW -end - -""" - TwoBandDownwellingRadiation(shortwave=nothing, longwave=nothing) - -Return a two-band model for downwelling radiation (split into a shortwave band -and a longwave band) that passes through the atmosphere and arrives at the surface of ocean -or sea ice. -""" -TwoBandDownwellingRadiation(; shortwave=nothing, longwave=nothing) = - TwoBandDownwellingRadiation(shortwave, longwave) - -Adapt.adapt_structure(to, tsdr::TwoBandDownwellingRadiation) = - TwoBandDownwellingRadiation(adapt(to, tsdr.shortwave), - adapt(to, tsdr.longwave)) - ##### ##### Chekpointing ##### import Oceananigans: prognostic_state, restore_prognostic_state! -function prognostic_state(atmos::PrescribedAtmosphere) +function prognostic_state(atmos::PrescribedAtmosphere) return (; clock = prognostic_state(atmos.clock)) end -function restore_prognostic_state!(atmos::PrescribedAtmosphere, state) +function restore_prognostic_state!(atmos::PrescribedAtmosphere, state) restore_prognostic_state!(atmos.clock, state.clock) update_state!(atmos) return atmos diff --git a/src/Atmospheres/prescribed_atmosphere_regridder.jl b/src/Atmospheres/prescribed_atmosphere_regridder.jl index da93d1766..9ed5ebe1f 100644 --- a/src/Atmospheres/prescribed_atmosphere_regridder.jl +++ b/src/Atmospheres/prescribed_atmosphere_regridder.jl @@ -2,14 +2,12 @@ function ComponentExchanger(atmosphere::PrescribedAtmosphere, grid) regridder = atmosphere_regridder(atmosphere, grid) - state = (; u = Field{Center, Center, Nothing}(grid), - v = Field{Center, Center, Nothing}(grid), - T = Field{Center, Center, Nothing}(grid), - p = Field{Center, Center, Nothing}(grid), - q = Field{Center, Center, Nothing}(grid), - ℐꜜˢʷ = Field{Center, Center, Nothing}(grid), - ℐꜜˡʷ = Field{Center, Center, Nothing}(grid), - Jᶜ = Field{Center, Center, Nothing}(grid)) + state = (; u = Field{Center, Center, Nothing}(grid), + v = Field{Center, Center, Nothing}(grid), + T = Field{Center, Center, Nothing}(grid), + p = Field{Center, Center, Nothing}(grid), + q = Field{Center, Center, Nothing}(grid), + Jᶜ = Field{Center, Center, Nothing}(grid)) return ComponentExchanger(state, regridder) end diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl index 576e81014..b85ecbf7a 100644 --- a/src/DataWrangling/ECCO/ECCO.jl +++ b/src/DataWrangling/ECCO/ECCO.jl @@ -2,7 +2,7 @@ module ECCO export ECCOMetadatum, ECCO_immersed_grid, adjusted_ECCO_tracers, initialize! export ECCO2Monthly, ECCO4Monthly, ECCO2Daily -export ECCOPrescribedAtmosphere +export ECCOPrescribedAtmosphere, ECCOPrescribedRadiation export ECCO2DarwinMonthly, ECCO4DarwinMonthly export retrieve_data @@ -370,6 +370,7 @@ end inpainted_metadata_path(metadata::ECCOMetadatum) = joinpath(metadata.dir, inpainted_metadata_filename(metadata)) include("ECCO_atmosphere.jl") +include("ECCO_radiation.jl") ##### ##### Column Field for ECCO datasets (which always download globally) diff --git a/src/DataWrangling/ECCO/ECCO_atmosphere.jl b/src/DataWrangling/ECCO/ECCO_atmosphere.jl index 3481c4e03..6f54ea1d4 100644 --- a/src/DataWrangling/ECCO/ECCO_atmosphere.jl +++ b/src/DataWrangling/ECCO/ECCO_atmosphere.jl @@ -1,6 +1,6 @@ using NumericalEarth.DataWrangling: DatasetBackend using Oceananigans.OutputReaders -using NumericalEarth.Atmospheres: PrescribedAtmosphere, TwoBandDownwellingRadiation +using NumericalEarth.Atmospheres: PrescribedAtmosphere """ ECCOPrescribedAtmosphere([architecture = CPU(), FT = Float32]; @@ -19,7 +19,10 @@ The atmospheric data will be held in `FieldTimeSeries` objects containing - air temperature and humidity: T, q - surface pressure: p - freshwater flux: rain -- downwelling radiation: ℐꜜˢʷ, ℐꜜˡʷ + +Note: downwelling shortwave / longwave radiation is now part of the +top-level `radiation` component. Use [`ECCOPrescribedRadiation`](@ref) to +load ECCO SW/LW into a `PrescribedRadiation`. """ function ECCOPrescribedAtmosphere(architecture = CPU(), FT = Float32; dataset = ECCO4Monthly(), @@ -36,8 +39,6 @@ function ECCOPrescribedAtmosphere(architecture = CPU(), FT = Float32; Ta_meta = Metadata(:air_temperature; dataset, start_date, end_date, dir) qa_meta = Metadata(:air_specific_humidity; dataset, start_date, end_date, dir) pa_meta = Metadata(:sea_level_pressure; dataset, start_date, end_date, dir) - ℐꜜˡʷ_meta = Metadata(:downwelling_longwave; dataset, start_date, end_date, dir) - ℐꜜˢʷ_meta = Metadata(:downwelling_shortwave; dataset, start_date, end_date, dir) Fr_meta = Metadata(:rain_freshwater_flux; dataset, start_date, end_date, dir) kw = (; time_indices_in_memory, time_indexing) @@ -48,10 +49,8 @@ function ECCOPrescribedAtmosphere(architecture = CPU(), FT = Float32; Ta = FieldTimeSeries(Ta_meta, architecture; kw...) qa = FieldTimeSeries(qa_meta, architecture; kw...) pa = FieldTimeSeries(pa_meta, architecture; kw...) - ℐꜜˡʷ = FieldTimeSeries(ℐꜜˡʷ_meta, architecture; kw...) - ℐꜜˢʷ = FieldTimeSeries(ℐꜜˢʷ_meta, architecture; kw...) Fr = FieldTimeSeries(Fr_meta, architecture; kw...) - + freshwater_flux = (; rain = Fr) times = ua.times @@ -61,8 +60,6 @@ function ECCOPrescribedAtmosphere(architecture = CPU(), FT = Float32; tracers = (T = Ta, q = qa) pressure = pa - downwelling_radiation = TwoBandDownwellingRadiation(shortwave=ℐꜜˢʷ, longwave=ℐꜜˡʷ) - FT = eltype(ua) surface_layer_height = convert(FT, surface_layer_height) @@ -70,7 +67,6 @@ function ECCOPrescribedAtmosphere(architecture = CPU(), FT = Float32; velocities, freshwater_flux, tracers, - downwelling_radiation, surface_layer_height, pressure) diff --git a/src/DataWrangling/ECCO/ECCO_radiation.jl b/src/DataWrangling/ECCO/ECCO_radiation.jl new file mode 100644 index 000000000..01a1eebcf --- /dev/null +++ b/src/DataWrangling/ECCO/ECCO_radiation.jl @@ -0,0 +1,44 @@ +using NumericalEarth.Radiations: PrescribedRadiation, SurfaceRadiationProperties, default_stefan_boltzmann_constant + +""" + ECCOPrescribedRadiation([architecture = CPU(), FT = Float32]; + dataset = ECCO4Monthly(), + start_date = first_date(dataset, :downwelling_shortwave), + end_date = last_date(dataset, :downwelling_shortwave), + dir = default_download_directory(dataset), + time_indices_in_memory = 10, + time_indexing = Cyclical(), + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0), + stefan_boltzmann_constant = default_stefan_boltzmann_constant, + other_kw...) + +Return a [`PrescribedRadiation`](@ref) backed by ECCO downwelling shortwave +and longwave fields. +""" +function ECCOPrescribedRadiation(architecture = CPU(), FT = Float32; + dataset = ECCO4Monthly(), + start_date = first_date(dataset, :downwelling_shortwave), + end_date = last_date(dataset, :downwelling_shortwave), + dir = default_download_directory(dataset), + time_indexing = Cyclical(), + time_indices_in_memory = 10, + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0), + stefan_boltzmann_constant = default_stefan_boltzmann_constant, + other_kw...) + + ℐꜜˢʷ_meta = Metadata(:downwelling_shortwave; dataset, start_date, end_date, dir) + ℐꜜˡʷ_meta = Metadata(:downwelling_longwave; dataset, start_date, end_date, dir) + + kw = (; time_indices_in_memory, time_indexing) + kw = merge(kw, other_kw) + + ℐꜜˢʷ = FieldTimeSeries(ℐꜜˢʷ_meta, architecture; kw...) + ℐꜜˡʷ = FieldTimeSeries(ℐꜜˡʷ_meta, architecture; kw...) + + return PrescribedRadiation(ℐꜜˢʷ, ℐꜜˡʷ; + ocean_surface, + sea_ice_surface, + stefan_boltzmann_constant) +end diff --git a/src/DataWrangling/JRA55/JRA55.jl b/src/DataWrangling/JRA55/JRA55.jl index 5751b9dcb..57cdcd6b4 100644 --- a/src/DataWrangling/JRA55/JRA55.jl +++ b/src/DataWrangling/JRA55/JRA55.jl @@ -1,6 +1,11 @@ module JRA55 -export JRA55FieldTimeSeries, JRA55PrescribedAtmosphere, JRA55PrescribedLand, RepeatYearJRA55, MultiYearJRA55 +export JRA55FieldTimeSeries, + JRA55PrescribedAtmosphere, + JRA55PrescribedLand, + JRA55PrescribedRadiation, + RepeatYearJRA55, + MultiYearJRA55 using Oceananigans using Oceananigans.Units @@ -14,9 +19,8 @@ using Oceananigans.OutputReaders: Cyclical, TotallyInMemory, AbstractInMemoryBac using NumericalEarth -using NumericalEarth.Atmospheres: - PrescribedAtmosphere, - TwoBandDownwellingRadiation +using NumericalEarth.Atmospheres: PrescribedAtmosphere +using NumericalEarth.Radiations: PrescribedRadiation, SurfaceRadiationProperties, default_stefan_boltzmann_constant using GPUArraysCore: @allowscalar @@ -40,5 +44,6 @@ include("JRA55_metadata.jl") include("JRA55_field_time_series.jl") include("JRA55_prescribed_atmosphere.jl") include("JRA55_prescribed_land.jl") +include("JRA55_prescribed_radiation.jl") end # module diff --git a/src/DataWrangling/JRA55/JRA55_prescribed_atmosphere.jl b/src/DataWrangling/JRA55/JRA55_prescribed_atmosphere.jl index ae395238d..373566c54 100644 --- a/src/DataWrangling/JRA55/JRA55_prescribed_atmosphere.jl +++ b/src/DataWrangling/JRA55/JRA55_prescribed_atmosphere.jl @@ -15,8 +15,11 @@ JRA55PrescribedAtmosphere(arch::Distributed, FT = Float32; kw...) = other_kw...) Return a [`PrescribedAtmosphere`](@ref) representing JRA55 reanalysis data. -The atmospheric data will be held in `JRA55FieldTimeSeries` objects containing. -For a detailed description of the keyword arguments, see the [`JRA55FieldTimeSeries`](@ref) constructor. +The atmospheric data will be held in `JRA55FieldTimeSeries` objects. + +Note: downwelling shortwave / longwave radiation is now part of the +top-level `radiation` component. Use [`JRA55PrescribedRadiation`](@ref) to +load JRA55 SW/LW into a `PrescribedRadiation`. """ function JRA55PrescribedAtmosphere(architecture = CPU(), FT = Float32; dataset = RepeatYearJRA55(), @@ -30,15 +33,13 @@ function JRA55PrescribedAtmosphere(architecture = CPU(), FT = Float32; kw = (; time_indexing, backend, start_date, end_date, dataset) kw = merge(kw, other_kw) - ua = JRA55FieldTimeSeries(:eastward_velocity, architecture, FT; kw...) - va = JRA55FieldTimeSeries(:northward_velocity, architecture, FT; kw...) - Ta = JRA55FieldTimeSeries(:temperature, architecture, FT; kw...) - qa = JRA55FieldTimeSeries(:specific_humidity, architecture, FT; kw...) - pa = JRA55FieldTimeSeries(:sea_level_pressure, architecture, FT; kw...) - Fra = JRA55FieldTimeSeries(:rain_freshwater_flux, architecture, FT; kw...) - Fsn = JRA55FieldTimeSeries(:snow_freshwater_flux, architecture, FT; kw...) - ℐꜜˡʷ = JRA55FieldTimeSeries(:downwelling_longwave_radiation, architecture, FT; kw...) - ℐꜜˢʷ = JRA55FieldTimeSeries(:downwelling_shortwave_radiation, architecture, FT; kw...) + ua = JRA55FieldTimeSeries(:eastward_velocity, architecture, FT; kw...) + va = JRA55FieldTimeSeries(:northward_velocity, architecture, FT; kw...) + Ta = JRA55FieldTimeSeries(:temperature, architecture, FT; kw...) + qa = JRA55FieldTimeSeries(:specific_humidity, architecture, FT; kw...) + pa = JRA55FieldTimeSeries(:sea_level_pressure, architecture, FT; kw...) + Fra = JRA55FieldTimeSeries(:rain_freshwater_flux, architecture, FT; kw...) + Fsn = JRA55FieldTimeSeries(:snow_freshwater_flux, architecture, FT; kw...) freshwater_flux = (rain = Fra, snow = Fsn) @@ -54,8 +55,6 @@ function JRA55PrescribedAtmosphere(architecture = CPU(), FT = Float32; pressure = pa - downwelling_radiation = TwoBandDownwellingRadiation(shortwave=ℐꜜˢʷ, longwave=ℐꜜˡʷ) - FT = eltype(ua) surface_layer_height = convert(FT, surface_layer_height) @@ -63,7 +62,6 @@ function JRA55PrescribedAtmosphere(architecture = CPU(), FT = Float32; velocities, freshwater_flux, tracers, - downwelling_radiation, surface_layer_height, pressure) diff --git a/src/DataWrangling/JRA55/JRA55_prescribed_radiation.jl b/src/DataWrangling/JRA55/JRA55_prescribed_radiation.jl new file mode 100644 index 000000000..da28be94a --- /dev/null +++ b/src/DataWrangling/JRA55/JRA55_prescribed_radiation.jl @@ -0,0 +1,42 @@ +JRA55PrescribedRadiation(arch::Distributed, FT = Float32; kw...) = + JRA55PrescribedRadiation(child_architecture(arch); kw...) + +""" + JRA55PrescribedRadiation([architecture = CPU(), FT = Float32]; + dataset = RepeatYearJRA55(), + start_date = first_date(dataset, :downwelling_shortwave_radiation), + end_date = last_date(dataset, :downwelling_shortwave_radiation), + backend = JRA55NetCDFBackend(10), + time_indexing = Cyclical(), + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0), + stefan_boltzmann_constant = default_stefan_boltzmann_constant, + other_kw...) + +Return a [`PrescribedRadiation`](@ref) backed by JRA55 downwelling shortwave +and longwave `JRA55FieldTimeSeries`. Surface radiative properties (albedo, +emissivity) for ocean and sea-ice surfaces default to standard values; pass +`ocean_surface = nothing` (or `sea_ice_surface = nothing`) to omit a surface. +""" +function JRA55PrescribedRadiation(architecture = CPU(), FT = Float32; + dataset = RepeatYearJRA55(), + start_date = first_date(dataset, :downwelling_shortwave_radiation), + end_date = last_date(dataset, :downwelling_shortwave_radiation), + backend = JRA55NetCDFBackend(10), + time_indexing = Cyclical(), + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0), + stefan_boltzmann_constant = default_stefan_boltzmann_constant, + other_kw...) + + kw = (; time_indexing, backend, start_date, end_date, dataset) + kw = merge(kw, other_kw) + + ℐꜜˢʷ = JRA55FieldTimeSeries(:downwelling_shortwave_radiation, architecture, FT; kw...) + ℐꜜˡʷ = JRA55FieldTimeSeries(:downwelling_longwave_radiation, architecture, FT; kw...) + + return PrescribedRadiation(ℐꜜˢʷ, ℐꜜˡʷ; + ocean_surface, + sea_ice_surface, + stefan_boltzmann_constant) +end diff --git a/src/DataWrangling/OSPapa/OSPapa.jl b/src/DataWrangling/OSPapa/OSPapa.jl index 0499ac0d5..056145019 100644 --- a/src/DataWrangling/OSPapa/OSPapa.jl +++ b/src/DataWrangling/OSPapa/OSPapa.jl @@ -1,6 +1,7 @@ module OSPapa export OSPapaPrescribedAtmosphere +export OSPapaPrescribedRadiation export os_papa_prescribed_fluxes export os_papa_prescribed_flux_boundary_conditions export OSPapaHourly @@ -14,7 +15,7 @@ using Downloads using Thermodynamics: q_vap_from_RH, Liquid using NumericalEarth.DataWrangling: download_progress -using NumericalEarth.Atmospheres: PrescribedAtmosphere, TwoBandDownwellingRadiation, AtmosphereThermodynamicsParameters +using NumericalEarth.Atmospheres: PrescribedAtmosphere, AtmosphereThermodynamicsParameters using NumericalEarth.Oceans: reference_density, heat_capacity using NumericalEarth.DataWrangling: @@ -80,6 +81,7 @@ end include("OSPapa_ocean_observations.jl") include("OSPapa_flux_observations.jl") include("OSPapa_prescribed_atmosphere.jl") +include("OSPapa_prescribed_radiation.jl") include("OSPapa_prescribed_fluxes.jl") end # module diff --git a/src/DataWrangling/OSPapa/OSPapa_prescribed_atmosphere.jl b/src/DataWrangling/OSPapa/OSPapa_prescribed_atmosphere.jl index c06d3f747..9990eee11 100644 --- a/src/DataWrangling/OSPapa/OSPapa_prescribed_atmosphere.jl +++ b/src/DataWrangling/OSPapa/OSPapa_prescribed_atmosphere.jl @@ -69,8 +69,6 @@ function OSPapaPrescribedAtmosphere(architecture = CPU(), FT = Float32; va = ospapa_fts(:northward_wind) Ta = ospapa_fts(:air_temperature) # K (Celsius conversion) Pa = ospapa_fts(:sea_level_pressure) # Pa (Millibar conversion) - swa = ospapa_fts(:shortwave_radiation) - lwa = ospapa_fts(:longwave_radiation) rain = ospapa_fts(:rain) # kg/m²/s (MillimetersPerHour conversion) thermo_params = AtmosphereThermodynamicsParameters(FT) @@ -82,7 +80,6 @@ function OSPapaPrescribedAtmosphere(architecture = CPU(), FT = Float32; tracers = (T=Ta, q=qa), pressure = Pa, freshwater_flux = (; rain), - downwelling_radiation = TwoBandDownwellingRadiation(shortwave=swa, longwave=lwa), thermodynamics_parameters = thermo_params, surface_layer_height = convert(FT, surface_layer_height)) end diff --git a/src/DataWrangling/OSPapa/OSPapa_prescribed_radiation.jl b/src/DataWrangling/OSPapa/OSPapa_prescribed_radiation.jl new file mode 100644 index 000000000..04ea5b798 --- /dev/null +++ b/src/DataWrangling/OSPapa/OSPapa_prescribed_radiation.jl @@ -0,0 +1,43 @@ +using NumericalEarth.Radiations: PrescribedRadiation, SurfaceRadiationProperties, default_stefan_boltzmann_constant + +""" + OSPapaPrescribedRadiation(architecture = CPU(), FT = Float32; + start_date = first_date(OSPapaHourly(), :shortwave_radiation), + end_date = last_date(OSPapaHourly(), :shortwave_radiation), + dir = download_OSPapa_cache, + max_gap_hours = 72, + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = nothing, + stefan_boltzmann_constant = default_stefan_boltzmann_constant) + +Construct a `PrescribedRadiation` from Ocean Station Papa buoy SW/LW +observations on a single-column grid. +""" +function OSPapaPrescribedRadiation(architecture = CPU(), FT = Float32; + start_date = first_date(OSPapaHourly(), :shortwave_radiation), + end_date = last_date(OSPapaHourly(), :shortwave_radiation), + dir = download_OSPapa_cache, + max_gap_hours = 72, + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = nothing, + stefan_boltzmann_constant = default_stefan_boltzmann_constant) + + mdkw = (; dataset = OSPapaHourly(), start_date, end_date, dir) + surface_grid = RectilinearGrid(architecture, FT; size=(), topology=(Flat, Flat, Flat)) + + function ospapa_fts(name) + md = Metadata(name; mdkw...) + download_dataset(md) + fts = FieldTimeSeries(md, surface_grid; time_indices_in_memory = length(md)) + fill_gaps!(fts; max_gap = max_gap_hours) + return fts + end + + swa = ospapa_fts(:shortwave_radiation) + lwa = ospapa_fts(:longwave_radiation) + + return PrescribedRadiation(swa, lwa; + ocean_surface, + sea_ice_surface, + stefan_boltzmann_constant) +end diff --git a/src/Diagnostics/meridional_heat_transport.jl b/src/Diagnostics/meridional_heat_transport.jl index 6bc6338db..9ef034fd2 100644 --- a/src/Diagnostics/meridional_heat_transport.jl +++ b/src/Diagnostics/meridional_heat_transport.jl @@ -61,7 +61,7 @@ sea_ice = sea_ice_simulation(grid, ocean) atmosphere = PrescribedAtmosphere(grid, [0.0]) -esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation = Radiation()) +esm = OceanSeaIceModel(sea_ice, ocean; atmosphere) mht = meridional_heat_transport(esm) diff --git a/src/EarthSystemModels/EarthSystemModels.jl b/src/EarthSystemModels/EarthSystemModels.jl index 5f55ea644..80de85734 100644 --- a/src/EarthSystemModels/EarthSystemModels.jl +++ b/src/EarthSystemModels/EarthSystemModels.jl @@ -8,8 +8,6 @@ export SimilarityTheoryFluxes, CoefficientBasedFluxes, FreezingLimitedOceanTemperature, - Radiation, - LatitudeDependentAlbedo, SkinTemperature, BulkTemperature, compute_atmosphere_ocean_fluxes!, @@ -80,10 +78,10 @@ const NoOceanInterface = ComponentInterfaces{<:Nothing, <:AtmosphereInterface, const NoAtmosInterface = ComponentInterfaces{<:Nothing, <:Nothing, <:SeaIceOceanInterface} const NoInterface = ComponentInterfaces{<:Nothing, <:Nothing, <:Nothing} -const NoSeaIceInterfaceModel = EarthSystemModel{I, A, L, O, <:NoSeaIceInterface} where {I, A, L, O} -const NoAtmosInterfaceModel = EarthSystemModel{I, A, L, O, <:NoAtmosInterface} where {I, A, L, O} -const NoOceanInterfaceModel = EarthSystemModel{I, A, L, O, <:NoOceanInterface} where {I, A, L, O} -const NoInterfaceModel = EarthSystemModel{I, A, L, O, <:NoInterface} where {I, A, L, O} +const NoSeaIceInterfaceModel = EarthSystemModel{R, A, L, I, O, <:NoSeaIceInterface} where {R, A, L, I, O} +const NoAtmosInterfaceModel = EarthSystemModel{R, A, L, I, O, <:NoAtmosInterface} where {R, A, L, I, O} +const NoOceanInterfaceModel = EarthSystemModel{R, A, L, I, O, <:NoOceanInterface} where {R, A, L, I, O} +const NoInterfaceModel = EarthSystemModel{R, A, L, I, O, <:NoInterface} where {R, A, L, I, O} InterfaceComputations.compute_atmosphere_sea_ice_fluxes!(::NoSeaIceInterfaceModel) = nothing InterfaceComputations.compute_sea_ice_ocean_fluxes!(::NoSeaIceInterfaceModel) = nothing diff --git a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl index ed93f8c9c..62fb8f15d 100644 --- a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl +++ b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl @@ -6,9 +6,7 @@ using Oceananigans.Utils: KernelParameters using Adapt export - Radiation, ComponentInterfaces, - LatitudeDependentAlbedo, SimilarityTheoryFluxes, MomentumRoughnessLength, ScalarRoughnessLength, @@ -41,6 +39,26 @@ import Oceananigans.Simulations: initialize! net_fluxes(::Nothing) = nothing +##### +##### Radiation hooks: declared here so the turbulent flux kernels can +##### resolve them at parse time. The `Radiations` module extends them +##### with concrete methods for `PrescribedRadiation`. +##### + +# `nothing` fallback (radiation is off). Concrete methods for +# `PrescribedRadiation` (and future radiation types) are added in `Radiations`. +@inline kernel_radiation_properties(::Nothing) = nothing + +@inline function air_sea_interface_radiation_state(::Nothing, ::Nothing, i, j, k, grid, time) + z = zero(eltype(grid)) + return (σ = z, α = z, ϵ = z, ℐꜜˢʷ = z, ℐꜜˡʷ = z) +end + +@inline function air_sea_ice_interface_radiation_state(::Nothing, ::Nothing, i, j, k, grid, time) + z = zero(eltype(grid)) + return (σ = z, α = z, ϵ = z, ℐꜜˢʷ = z, ℐꜜˡʷ = z) +end + ##### ##### Utilities ##### @@ -63,11 +81,6 @@ function interface_kernel_parameters(grid) return kernel_parameters end -# Radiation -include("radiation.jl") -include("latitude_dependent_albedo.jl") -include("tabulated_albedo.jl") - # Turbulent fluxes include("roughness_lengths.jl") include("interface_states.jl") diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_ocean_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_ocean_fluxes.jl index a32f36c01..f0ea3af2b 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_ocean_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_ocean_fluxes.jl @@ -11,7 +11,7 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) # Simplify NamedTuple to reduce parameter space consumption. # See https://github.com/CliMA/NumericalEarth.jl/issues/116. - atmosphere_data = merge(atmosphere_fields, + atmosphere_data = merge(atmosphere_fields, (; h_bℓ = boundary_layer_height(coupled_model.atmosphere))) flux_formulation = coupled_model.interfaces.atmosphere_ocean_interface.flux_formulation @@ -23,6 +23,15 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) surface_layer_height = surface_layer_height(coupled_model.atmosphere), gravitational_acceleration = coupled_model.interfaces.properties.gravitational_acceleration) + # Radiation state for the interface solve (used by SkinTemperature). + # When `radiation === nothing` these are `nothing`s and the getter + # returns zero-valued radiative state, so SkinTemperature degrades to + # a turbulent-only flux balance. + radiation = coupled_model.radiation + radiation_kernel_props = kernel_radiation_properties(radiation) + radiation_exchanger = exchanger.radiation + radiation_state = isnothing(radiation_exchanger) ? nothing : radiation_exchanger.state + kernel_parameters = interface_kernel_parameters(grid) launch!(arch, grid, kernel_parameters, @@ -36,12 +45,14 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) atmosphere_data, interface_properties, atmosphere_properties, - ocean_properties) + ocean_properties, + radiation_kernel_props, + radiation_state) return nothing end -""" Compute turbulent fluxes between an atmosphere and a interface state using similarity theory """ +""" Compute turbulent fluxes between an atmosphere and an interface state using similarity theory """ @kernel function _compute_atmosphere_ocean_interface_state!(interface_fluxes, interface_temperature, grid, @@ -51,7 +62,9 @@ end atmosphere_state, interface_properties, atmosphere_properties, - ocean_properties) + ocean_properties, + radiation_kernel_props, + radiation_exchanger_state) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) # index of the top ocean cell @@ -63,11 +76,8 @@ end Tᵃᵗ = atmosphere_state.T[i, j, 1] pᵃᵗ = atmosphere_state.p[i, j, 1] qᵃᵗ = atmosphere_state.q[i, j, 1] - ℐꜜˢʷ = atmosphere_state.ℐꜜˢʷ[i, j, 1] - ℐꜜˡʷ = atmosphere_state.ℐꜜˡʷ[i, j, 1] - # Extract state variables at cell centers - # Ocean state + # Ocean state at cell centers uᵒᶜ = ℑxᶜᵃᵃ(i, j, kᴺ, grid, interior_state.u) vᵒᶜ = ℑyᵃᶜᵃ(i, j, kᴺ, grid, interior_state.v) Tᵒᶜ = interior_state.T[i, j, kᴺ] @@ -76,8 +86,6 @@ end end # Build thermodynamic and dynamic states in the atmosphere and interface. - # Notation: - # ⋅ 𝒰 ≡ "dynamic" state vector (thermodynamics + reference height + velocity) ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters zᵃᵗ = atmosphere_properties.surface_layer_height # elevation of atmos variables relative to interface @@ -90,7 +98,12 @@ end h_bℓ = atmosphere_state.h_bℓ) local_interior_state = (u=uᵒᶜ, v=vᵒᶜ, T=Tᵒᶜ, S=Sᵒᶜ) - downwelling_radiation = (; ℐꜜˢʷ, ℐꜜˡʷ) + + # Local radiative state at this cell. Returns zero-valued state when + # radiation is off. + radiation_state = air_sea_interface_radiation_state(radiation_kernel_props, + radiation_exchanger_state, + i, j, kᴺ, grid, time) # Estimate initial interface state FT = typeof(Tᵒᶜ) @@ -106,16 +119,6 @@ end needs_to_converge = stop_criteria isa ConvergenceStopCriteria not_water = inactive_node(i, j, kᴺ, grid, Center(), Center(), Center()) - # Compute local radiative properties and rebuild the interface properties - α = stateindex(interface_properties.radiation.α, i, j, kᴺ, grid, time, (Center, Center, Center), ℐꜜˢʷ) - ϵ = stateindex(interface_properties.radiation.ϵ, i, j, kᴺ, grid, time, (Center, Center, Center)) - σ = interface_properties.radiation.σ - - interface_properties = InterfaceProperties((; α, ϵ, σ), - interface_properties.specific_humidity_formulation, - interface_properties.temperature_formulation, - interface_properties.velocity_formulation) - if needs_to_converge && not_water interface_state = zero_interface_state(FT) else @@ -123,7 +126,7 @@ end initial_interface_state, local_atmosphere_state, local_interior_state, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, ocean_properties) diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index 1292d6311..420b6ca82 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -14,9 +14,7 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model) atmosphere_fields = exchanger.atmosphere.state - # Simplify NamedTuple to reduce parameter space consumption. - # See https://github.com/CliMA/NumericalEarth.jl/issues/116. - atmosphere_data = merge(atmosphere_fields, + atmosphere_data = merge(atmosphere_fields, (; h_bℓ = boundary_layer_height(coupled_model.atmosphere))) flux_formulation = coupled_model.interfaces.atmosphere_sea_ice_interface.flux_formulation @@ -30,6 +28,11 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model) surface_layer_height = surface_layer_height(coupled_model.atmosphere), gravitational_acceleration = coupled_model.interfaces.properties.gravitational_acceleration) + radiation = coupled_model.radiation + radiation_kernel_props = kernel_radiation_properties(radiation) + radiation_exchanger = exchanger.radiation + radiation_state = isnothing(radiation_exchanger) ? nothing : radiation_exchanger.state + kernel_parameters = interface_kernel_parameters(grid) launch!(arch, grid, kernel_parameters, @@ -44,12 +47,14 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model) interface_properties, atmosphere_properties, sea_ice_properties, - ocean_properties) + ocean_properties, + radiation_kernel_props, + radiation_state) return nothing end -""" Compute turbulent fluxes between an atmosphere and a interface state using similarity theory """ +""" Compute turbulent fluxes between an atmosphere and an interface state using similarity theory """ @kernel function _compute_atmosphere_sea_ice_interface_state!(interface_fluxes, interface_temperature, grid, @@ -60,11 +65,14 @@ end interface_properties, atmosphere_properties, sea_ice_properties, - ocean_properties) + ocean_properties, + radiation_kernel_props, + radiation_exchanger_state) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) # index of the top ocean cell FT = eltype(grid) + time = Time(clock.time) @inbounds begin uᵃᵗ = atmosphere_state.u[i, j, 1] @@ -72,10 +80,7 @@ end Tᵃᵗ = atmosphere_state.T[i, j, 1] pᵃᵗ = atmosphere_state.p[i, j, 1] qᵃᵗ = atmosphere_state.q[i, j, 1] - ℐꜜˢʷ = atmosphere_state.ℐꜜˢʷ[i, j, 1] - ℐꜜˡʷ = atmosphere_state.ℐꜜˡʷ[i, j, 1] - # Extract state variables at cell centers # Ocean properties below sea ice Tᵒᶜ = interior_state.Tᵒᶜ[i, j, kᴺ] Tᵒᶜ = convert_to_kelvin(ocean_properties.temperature_units, Tᵒᶜ) @@ -91,9 +96,6 @@ end Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ) end - # Build thermodynamic and dynamic states in the atmosphere and interface. - # Notation: - # ⋅ 𝒰 ≡ "dynamic" state vector (thermodynamics + reference height + velocity) ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters zᵃᵗ = atmosphere_properties.surface_layer_height # elevation of atmos variables relative to interface @@ -105,9 +107,14 @@ end q = qᵃᵗ, h_bℓ = atmosphere_state.h_bℓ) - downwelling_radiation = (; ℐꜜˢʷ, ℐꜜˡʷ) local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, h=hˢⁱ, hc=hc) - + + # Local radiative state at this cell. Returns zero-valued state when + # radiation is off. + radiation_state = air_sea_ice_interface_radiation_state(radiation_kernel_props, + radiation_exchanger_state, + i, j, kᴺ, grid, time) + # Estimate initial interface state (FP32 compatible) u★ = convert(FT, 1f-4) @@ -131,7 +138,7 @@ end initial_interface_state, local_atmosphere_state, local_interior_state, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, sea_ice_properties) diff --git a/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl b/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl index 65d6adbe8..f97404b09 100644 --- a/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl +++ b/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl @@ -70,22 +70,18 @@ struct AtmosphereOceanFluxes{F} friction_velocity :: F temperature_scale :: F water_vapor_scale :: F - upwelling_longwave :: F - downwelling_longwave :: F - downwelling_shortwave :: F end function AtmosphereOceanFluxes(grid) F = Field{Center, Center, Nothing} return AtmosphereOceanFluxes(F(grid), F(grid), F(grid), - F(grid), F(grid), F(grid), F(grid), F(grid), F(grid), F(grid), F(grid)) end -AtmosphereOceanFluxes(::Nothing) = AtmosphereOceanFluxes(ntuple(_ -> ZeroField(), 11)...) +AtmosphereOceanFluxes(::Nothing) = AtmosphereOceanFluxes(ntuple(_ -> ZeroField(), 8)...) -Adapt.adapt_structure(to, fluxes::AtmosphereOceanFluxes) = +Adapt.adapt_structure(to, fluxes::AtmosphereOceanFluxes) = AtmosphereOceanFluxes(Adapt.adapt(to, fluxes.latent_heat), Adapt.adapt(to, fluxes.sensible_heat), Adapt.adapt(to, fluxes.water_vapor), @@ -93,12 +89,9 @@ Adapt.adapt_structure(to, fluxes::AtmosphereOceanFluxes) = Adapt.adapt(to, fluxes.y_momentum), Adapt.adapt(to, fluxes.friction_velocity), Adapt.adapt(to, fluxes.temperature_scale), - Adapt.adapt(to, fluxes.water_vapor_scale), - Adapt.adapt(to, fluxes.upwelling_longwave), - Adapt.adapt(to, fluxes.downwelling_longwave), - Adapt.adapt(to, fluxes.downwelling_shortwave)) + Adapt.adapt(to, fluxes.water_vapor_scale)) -on_architecture(arch, fluxes::AtmosphereOceanFluxes) = +on_architecture(arch, fluxes::AtmosphereOceanFluxes) = AtmosphereOceanFluxes(on_architecture(arch, fluxes.latent_heat), on_architecture(arch, fluxes.sensible_heat), on_architecture(arch, fluxes.water_vapor), @@ -106,10 +99,7 @@ on_architecture(arch, fluxes::AtmosphereOceanFluxes) = on_architecture(arch, fluxes.y_momentum), on_architecture(arch, fluxes.friction_velocity), on_architecture(arch, fluxes.temperature_scale), - on_architecture(arch, fluxes.water_vapor_scale), - on_architecture(arch, fluxes.upwelling_longwave), - on_architecture(arch, fluxes.downwelling_longwave), - on_architecture(arch, fluxes.downwelling_shortwave)) + on_architecture(arch, fluxes.water_vapor_scale)) struct AtmosphereSeaIceFluxes{F} latent_heat :: F @@ -174,7 +164,8 @@ on_architecture(arch, fluxes::SeaIceOceanFluxes) = # ZeroFluxes is returned by computed_fluxes(::Nothing) for absent interfaces. # It contains the union of all flux field names across interface types. struct ZeroFluxes{Z} - # Atmosphere-ocean and atmosphere-sea-ice flux fields + # Atmosphere-ocean and atmosphere-sea-ice flux fields (turbulent only; + # radiative diagnostic fields live on the radiation component) latent_heat :: Z sensible_heat :: Z water_vapor :: Z @@ -183,16 +174,13 @@ struct ZeroFluxes{Z} friction_velocity :: Z temperature_scale :: Z water_vapor_scale :: Z - upwelling_longwave :: Z - downwelling_longwave :: Z - downwelling_shortwave :: Z # Sea ice-ocean flux fields interface_heat :: Z frazil_heat :: Z salt :: Z end -ZeroFluxes() = ZeroFluxes(ntuple(_ -> ZeroField(), 14)...) +ZeroFluxes() = ZeroFluxes(ntuple(_ -> ZeroField(), 11)...) @inline computed_fluxes(::Nothing) = ZeroFluxes() @@ -224,10 +212,9 @@ atmosphere_ocean_interface(grid, ::Nothing, ocean, args...) = nothing atmosphere_ocean_interface(grid, ::Nothing, ::Nothing, args...) = nothing atmosphere_ocean_interface(grid, atmosphere, ::Nothing, args...) = nothing -function atmosphere_ocean_interface(grid, +function atmosphere_ocean_interface(grid, atmosphere, ocean, - radiation, ao_flux_formulation, temperature_formulation, velocity_formulation, @@ -235,13 +222,7 @@ function atmosphere_ocean_interface(grid, ao_fluxes = AtmosphereOceanFluxes(grid) - σ = radiation.stefan_boltzmann_constant - αₐₒ = radiation.reflection.ocean - ϵₐₒ = radiation.emission.ocean - radiation = (σ=σ, α=αₐₒ, ϵ=ϵₐₒ) - - ao_properties = InterfaceProperties(radiation, - specific_humidity_formulation, + ao_properties = InterfaceProperties(specific_humidity_formulation, temperature_formulation, velocity_formulation) @@ -258,26 +239,19 @@ atmosphere_sea_ice_interface(grid, atmos, ::Nothing, args...) = nothing atmosphere_sea_ice_interface(grid, ::Nothing, sea_ice, args...) = nothing atmosphere_sea_ice_interface(grid, ::Nothing, ::Nothing, args...) = nothing -function atmosphere_sea_ice_interface(grid, +function atmosphere_sea_ice_interface(grid, atmosphere, sea_ice, - radiation, ai_flux_formulation, temperature_formulation, velocity_formulation) fluxes = AtmosphereSeaIceFluxes(grid) - σ = radiation.stefan_boltzmann_constant - αₐᵢ = radiation.reflection.sea_ice - ϵₐᵢ = radiation.emission.sea_ice - radiation = (σ=σ, α=αₐᵢ, ϵ=ϵₐᵢ) - phase = AtmosphericThermodynamics.Ice() specific_humidity_formulation = ImpureSaturationSpecificHumidity(phase) - properties = InterfaceProperties(radiation, - specific_humidity_formulation, + properties = InterfaceProperties(specific_humidity_formulation, temperature_formulation, velocity_formulation) @@ -362,6 +336,7 @@ Keyword Arguments - `ThreeEquationHeatFlux()`: coupled heat/salt/freezing point system (default) - `radiation`: radiation component. Default: `Radiation()`. ++ `radiation`: radiation component. Default: `nothing`. - `freshwater_density`: reference density of freshwater. Default: `default_freshwater_density`. - `atmosphere_ocean_fluxes`: flux formulation for atmosphere-ocean interface. Default: `SimilarityTheoryFluxes()`. - `atmosphere_sea_ice_fluxes`: flux formulation for atmosphere-sea ice interface. Default: `SimilarityTheoryFluxes()`. @@ -377,9 +352,9 @@ Keyword Arguments - `gravitational_acceleration`: gravitational acceleration. Default: `default_gravitational_acceleration`. """ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; + radiation = nothing, land = nothing, exchange_grid = exchange_grid(atmosphere, ocean, sea_ice), - radiation = Radiation(), freshwater_density = default_freshwater_density, atmosphere_ocean_fluxes = SimilarityTheoryFluxes(eltype(exchange_grid)), atmosphere_sea_ice_fluxes = atmosphere_sea_ice_similarity_theory(eltype(exchange_grid)), @@ -429,7 +404,6 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; ao_interface = atmosphere_ocean_interface(exchange_grid, atmosphere, ocean, - radiation, atmosphere_ocean_fluxes, atmosphere_ocean_interface_temperature, atmosphere_ocean_velocity_difference, @@ -437,10 +411,9 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; io_interface = sea_ice_ocean_interface(exchange_grid, sea_ice, ocean, sea_ice_ocean_heat_flux) - ai_interface = atmosphere_sea_ice_interface(exchange_grid, + ai_interface = atmosphere_sea_ice_interface(exchange_grid, atmosphere, sea_ice, - radiation, atmosphere_sea_ice_fluxes, atmosphere_sea_ice_interface_temperature, atmosphere_sea_ice_velocity_difference) @@ -449,7 +422,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; sea_ice = net_fluxes(sea_ice), atmosphere = net_fluxes(atmosphere)) - exchanger = StateExchanger(exchange_grid, atmosphere, land, ocean, sea_ice) + exchanger = StateExchanger(exchange_grid, radiation, atmosphere, land, ocean, sea_ice) properties = (; gravitational_acceleration) diff --git a/src/EarthSystemModels/InterfaceComputations/compute_interface_state.jl b/src/EarthSystemModels/InterfaceComputations/compute_interface_state.jl index e52b5a53a..b8703a79d 100644 --- a/src/EarthSystemModels/InterfaceComputations/compute_interface_state.jl +++ b/src/EarthSystemModels/InterfaceComputations/compute_interface_state.jl @@ -32,7 +32,7 @@ end initial_interface_state, atmosphere_state, interior_state, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, interior_properties) @@ -47,7 +47,7 @@ end Ψₛ⁻ = Ψₛⁿ Ψₛⁿ = iterate_interface_state(flux_formulation, Ψₛ⁻, Ψₐ, Ψᵢ, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, interior_properties) @@ -70,7 +70,7 @@ and interior properties `ℙₛ`, `ℙₐ`, and `ℙᵢ`. approximate_interface_state, atmosphere_state, interior_state, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, interior_properties) @@ -79,7 +79,7 @@ and interior properties `ℙₛ`, `ℙₐ`, and `ℙᵢ`. approximate_interface_state, atmosphere_state, interior_state, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, interior_properties) diff --git a/src/EarthSystemModels/InterfaceComputations/interface_states.jl b/src/EarthSystemModels/InterfaceComputations/interface_states.jl index 4f1ae6c2d..255f00bbb 100644 --- a/src/EarthSystemModels/InterfaceComputations/interface_states.jl +++ b/src/EarthSystemModels/InterfaceComputations/interface_states.jl @@ -9,8 +9,7 @@ using Thermodynamics: Liquid, Ice ##### Interface properties ##### -struct InterfaceProperties{R, Q, T, V} - radiation :: R +struct InterfaceProperties{Q, T, V} specific_humidity_formulation :: Q temperature_formulation :: T velocity_formulation :: V @@ -313,7 +312,7 @@ end interface_state, atmosphere_state, interior_state, - downwelling_radiation, + radiation_state, interface_properties, atmosphere_properties, interior_properties) @@ -329,16 +328,18 @@ end #ℰv = 0 #AtmosphericThermodynamics.latent_heat_vapor(ℂᵃᵗ, Tᵃᵗ) ℒⁱ = AtmosphericThermodynamics.latent_heat_sublim(ℂᵃᵗ, Tᵃᵗ) - # upwelling radiation is calculated explicitly + # upwelling radiation is calculated explicitly. radiation_state is + # produced by `air_sea_interface_radiation_state` (or its sea-ice + # variant) and contains zero-valued σ/α/ϵ/SW/LW when radiation is off. Tₛ⁻ = interface_state.T # approximate interface temperature from previous iteration - σ = interface_properties.radiation.σ - ϵ = interface_properties.radiation.ϵ - α = interface_properties.radiation.α - - ℐꜜˢʷ = downwelling_radiation.ℐꜜˢʷ - ℐꜜˡʷ = downwelling_radiation.ℐꜜˡʷ - ℐꜛˡʷ = emitted_longwave_radiation(Tₛ⁻, σ, ϵ) - Qd = net_absorbed_interface_radiation(ℐꜜˢʷ, ℐꜜˡʷ, α, ϵ) + σ = radiation_state.σ + ϵ = radiation_state.ϵ + α = radiation_state.α + + ℐꜜˢʷ = radiation_state.ℐꜜˢʷ + ℐꜜˡʷ = radiation_state.ℐꜜˡʷ + ℐꜛˡʷ = σ * ϵ * Tₛ⁻^4 + Qd = - (1 - α) * ℐꜜˢʷ - ϵ * ℐꜜˡʷ u★ = interface_state.u★ θ★ = interface_state.θ★ diff --git a/src/EarthSystemModels/InterfaceComputations/radiation.jl b/src/EarthSystemModels/InterfaceComputations/radiation.jl deleted file mode 100644 index 0e60d56af..000000000 --- a/src/EarthSystemModels/InterfaceComputations/radiation.jl +++ /dev/null @@ -1,114 +0,0 @@ -using Oceananigans.Grids: φnode - -@inline hack_cosd(φ) = cos(π * φ / 180) -@inline hack_sind(φ) = sin(π * φ / 180) - -struct Radiation{FT, E, R} - emission :: E - reflection :: 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)) - -""" - Radiation([arch = CPU(), FT=Float64]; - ocean_emissivity = 0.97, - sea_ice_emissivity = 1.0, - ocean_albedo = 0.05, - sea_ice_albedo = 0.7, - stefan_boltzmann_constant = 5.67e-8) - -Constructs a `Radiation` object that represents the radiation properties of the ocean and sea ice. - -Arguments -========= - -- `arch`: The architecture of the system. Default: `CPU()`. -- `FT`: The floating-point type to use. Default: `Float64`. - -Keyword Arguments -================= - -- `ocean_emissivity`: The emissivity of the ocean surface. Default: `0.97`. -- `sea_ice_emissivity`: The emissivity of the sea ice surface. Default: `1.0`. -- `ocean_albedo`: The albedo of the ocean surface. Default: `0.05`. -- `sea_ice_albedo`: The albedo of the sea ice surface. Default: `0.7`. -- `stefan_boltzmann_constant`: The Stefan-Boltzmann constant. Default: `5.67e-8`. -""" -function Radiation(arch = CPU(), FT=Oceananigans.defaults.FloatType; - ocean_emissivity = 0.97, - sea_ice_emissivity = 1.0, - ocean_albedo = 0.05, - sea_ice_albedo = 0.7, - stefan_boltzmann_constant = 5.67e-8) - - ocean_emissivity isa Number && (ocean_emissivity = convert(FT, ocean_emissivity)) - sea_ice_emissivity isa Number && (sea_ice_emissivity = convert(FT, sea_ice_emissivity)) - ocean_albedo isa Number && (ocean_albedo = convert(FT, ocean_albedo)) - sea_ice_albedo isa Number && (sea_ice_albedo = convert(FT, sea_ice_albedo)) - - emission = SurfaceProperties(ocean_emissivity, sea_ice_emissivity) - reflection = SurfaceProperties(ocean_albedo, sea_ice_albedo) - - return Radiation(emission, - reflection, - convert(FT, stefan_boltzmann_constant)) -end - -Base.summary(r::Radiation{FT}) where FT = "Radiation{$FT}" - -function Base.show(io::IO, r::Radiation) - σ = r.stefan_boltzmann_constant - - print(io, summary(r), ":", '\n') - print(io, "├── stefan_boltzmann_constant: ", prettysummary(σ), '\n') - print(io, "├── emission: ", summary(r.emission), '\n') - print(io, "│ ├── ocean: ", prettysummary(r.emission.ocean), '\n') - print(io, "│ └── sea_ice: ", prettysummary(r.emission.ocean), '\n') - print(io, "└── reflection: ", summary(r.reflection), '\n') - print(io, " ├── ocean: ", prettysummary(r.reflection.ocean), '\n') - print(io, " └── sea_ice: ", prettysummary(r.reflection.sea_ice)) -end - -struct SurfaceProperties{O, I} - ocean :: O - sea_ice :: I -end - -Adapt.adapt_structure(to, s :: SurfaceProperties) = - SurfaceProperties(Adapt.adapt(to, s.ocean), - Adapt.adapt(to, s.sea_ice)) - -Base.summary(properties::SurfaceProperties) = "SurfaceProperties" - -function Base.show(io::IO, properties::SurfaceProperties) - print(io, "SurfaceProperties:", '\n') - print(io, "├── ocean: ", summary(properties.ocean), '\n') - print(io, "└── sea_ice: ", summary(properties.sea_ice)) -end - -const CCC = (Center, Center, Center) - -@inline function emitted_longwave_radiation(i, j, k, grid, time, T, σ, ϵ) - ϵi = stateindex(ϵ, i, j, k, grid, time, CCC) - return σ * ϵi * T^4 -end - -# Split the individual bands -@inline function absorbed_longwave_radiation(i, j, k, grid, time, ϵ, ℐꜜˡʷ) - ϵi = stateindex(ϵ, i, j, k, grid, time, CCC) - return - ϵi * ℐꜜˡʷ -end - -@inline function transmitted_shortwave_radiation(i, j, k, grid, time, α, ℐꜜˢʷ) - αi = stateindex(α, i, j, k, grid, time, CCC, ℐꜜˢʷ) - return - (1 - αi) * ℐꜜˢʷ -end - -# Inside the solver we lose both spatial and temporal information, but the -# radiative properties have already been computed correctly -@inline net_absorbed_interface_radiation(ℐꜜˢʷ, ℐꜜˡʷ, α, ϵ) = - (1 - α) * ℐꜜˢʷ - ϵ * ℐꜜˡʷ -@inline emitted_longwave_radiation(T, σ, ϵ) = σ * ϵ * T^4 diff --git a/src/EarthSystemModels/InterfaceComputations/state_exchanger.jl b/src/EarthSystemModels/InterfaceComputations/state_exchanger.jl index fadf03585..85830cb4d 100644 --- a/src/EarthSystemModels/InterfaceComputations/state_exchanger.jl +++ b/src/EarthSystemModels/InterfaceComputations/state_exchanger.jl @@ -2,8 +2,9 @@ ComponentExchanger(component, exchange_grid) Holds a regridder and a buffer of `state` fields used to bring data from a -component (atmosphere, land, ocean, sea ice) onto a shared `exchange_grid`, -where atmosphere--ocean and atmosphere--sea-ice fluxes are computed. +component (radiation, atmosphere, land, ocean, sea ice) onto a shared +`exchange_grid`, where atmosphere--ocean and atmosphere--sea-ice fluxes are +computed. """ struct ComponentExchanger{S, EX} state :: S @@ -11,35 +12,39 @@ struct ComponentExchanger{S, EX} end """ - StateExchanger(grid, atmosphere, land, ocean, sea_ice) + StateExchanger(grid, radiation, atmosphere, land, ocean, sea_ice) Container for one `ComponentExchanger` per component. The `grid` is the shared exchange grid onto which each component's state is regridded each time step. """ -struct StateExchanger{G, A, L, O, S} +struct StateExchanger{G, R, A, L, O, S} grid :: G + radiation :: R atmosphere :: A land :: L ocean :: O sea_ice :: S - function StateExchanger(grid, atmosphere, land, ocean, sea_ice) + function StateExchanger(grid, radiation, atmosphere, land, ocean, sea_ice) + radiation_exchanger = ComponentExchanger(radiation, grid) atmosphere_exchanger = ComponentExchanger(atmosphere, grid) land_exchanger = ComponentExchanger(land, grid) ocean_exchanger = ComponentExchanger(ocean, grid) sea_ice_exchanger = ComponentExchanger(sea_ice, grid) G = typeof(grid) + R = typeof(radiation_exchanger) A = typeof(atmosphere_exchanger) L = typeof(land_exchanger) O = typeof(ocean_exchanger) S = typeof(sea_ice_exchanger) - return new{G, A, L, O, S}(grid, - atmosphere_exchanger, - land_exchanger, - ocean_exchanger, - sea_ice_exchanger) + return new{G, R, A, L, O, S}(grid, + radiation_exchanger, + atmosphere_exchanger, + land_exchanger, + ocean_exchanger, + sea_ice_exchanger) end end @@ -47,6 +52,7 @@ end ComponentExchanger(::Nothing, grid) = nothing function initialize!(exchanger::StateExchanger, model) + initialize!(exchanger.radiation, exchanger.grid, model.radiation) initialize!(exchanger.atmosphere, exchanger.grid, model.atmosphere) initialize!(exchanger.land, exchanger.grid, model.land) initialize!(exchanger.ocean, exchanger.grid, model.ocean) diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl index 113cce6b5..1e7f59c45 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -4,9 +4,10 @@ using Oceananigans: SeawaterBuoyancy using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using KernelAbstractions: @kernel, @index -mutable struct EarthSystemModel{I, A, L, O, F, C, Arch} <: AbstractModel{Nothing, Arch} +mutable struct EarthSystemModel{R, A, L, I, O, F, C, Arch} <: AbstractModel{Nothing, Arch} architecture :: Arch clock :: C + radiation :: R atmosphere :: A land :: L sea_ice :: I @@ -38,6 +39,7 @@ function Base.show(io::IO, cm::ESM) ocean_summary = summary(cm.ocean) end + print(io, "├── radiation: ", summary(cm.radiation), "\n") print(io, "├── atmosphere: ", summary(cm.atmosphere), "\n") print(io, "├── land: ", summary(cm.land), "\n") print(io, "├── sea_ice: ", sea_ice_summary, "\n") @@ -80,9 +82,14 @@ reference_density(unsupported) = heat_capacity(unsupported) = throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))")) +# Hook called after `interfaces` is constructed and the exchange grid is known. +# Concrete radiation types (e.g. `PrescribedRadiation`) overload this to +# allocate `interface_fluxes` per-surface on the exchange grid. +allocate_interface_fluxes!(::Any, exchange_grid, surfaces) = nothing +allocate_interface_fluxes!(::Nothing, exchange_grid, surfaces) = nothing + """ - EarthSystemModel(atmosphere, ocean, sea_ice; - radiation = Radiation(), + EarthSystemModel(radiation, atmosphere, land, sea_ice, ocean; clock = Clock{Float64}(time=0), ocean_reference_density = reference_density(ocean), ocean_heat_capacity = heat_capacity(ocean), @@ -90,26 +97,32 @@ heat_capacity(unsupported) = sea_ice_heat_capacity = heat_capacity(sea_ice), interfaces = nothing) -Construct a coupled earth system model with an atmosphere, ocean, and sea ice component. -For simpler configurations, see [`OceanOnlyModel`](@ref) and [`OceanSeaIceModel`](@ref). +Construct a coupled earth system model. Components are passed in struct order +(top to bottom): radiation, atmosphere, land, sea_ice, ocean. Pass `nothing` +for components that are absent. For simpler configurations, see +[`OceanOnlyModel`](@ref), [`OceanSeaIceModel`](@ref), and +[`AtmosphereOceanModel`](@ref). Arguments ========== -- `atmosphere`: A representation of a possibly time-dependent atmospheric state. +- `radiation`: Radiation component, or `nothing` for a radiatively decoupled surface. + Pass a `PrescribedRadiation` (e.g. `JRA55PrescribedRadiation(...)`) to + enable radiative forcing. +- `atmosphere`: A representation of a possibly time-dependent atmospheric state, or `nothing`. For a prognostic atmosphere, use `atmosphere_simulation`. For prescribed atmospheric forcing, use `JRA55PrescribedAtmosphere` or `PrescribedAtmosphere`. -- `ocean`: A representation of a possibly time-dependent ocean state. Currently, only `Oceananigans.Simulation`s - of `Oceananigans.HydrostaticFreeSurfaceModel` are tested. -- `sea_ice`: A representation of a possibly time-dependent sea ice state. +- `land`: Land component, or `nothing`. +- `sea_ice`: A representation of a possibly time-dependent sea ice state, or `nothing`. For example, the minimalist `FreezingLimitedOceanTemperature` represents oceanic latent heating during freezing only, but does not evolve sea ice variables. For prognostic sea ice use an `Oceananigans.Simulation` of `ClimaSeaIce.SeaIceModel`. +- `ocean`: A representation of a possibly time-dependent ocean state. Currently, only `Oceananigans.Simulation`s + of `Oceananigans.HydrostaticFreeSurfaceModel` are tested. Keyword Arguments ================== -- `radiation`: Radiation component used to compute surface fluxes at the bottom of the atmosphere. - `clock`: Keeps track of time. - `ocean_reference_density`: Reference density for the ocean. Defaults to value from ocean model. - `ocean_heat_capacity`: Heat capacity for the ocean. Defaults to value from ocean model. @@ -117,59 +130,30 @@ Keyword Arguments - `sea_ice_heat_capacity`: Heat capacity for sea ice. Defaults to value from sea ice model. - `interfaces`: Component interfaces for coupling. Defaults to `nothing` and will be constructed automatically. To customize the sea ice-ocean heat flux formulation, create interfaces manually using `ComponentInterfaces`. - -Stability Functions -==================== - -The model uses similarity theory for turbulent fluxes between components. You can customize the stability functions -by creating a new [`SimilarityTheoryFluxes`](@ref) object with your desired stability functions. For example: - -```jldoctest earth_system_model -using NumericalEarth -using Oceananigans - -grid = RectilinearGrid(size=10, z=(-100, 0), topology=(Flat, Flat, Bounded)) -ocean = ocean_simulation(grid, timestepper = :QuasiAdamsBashforth2) - -# Three choices for stability function: -# "No stability function", which also apply to neutral boundary layers -stability_functions = nothing - -# Atmosphere-sea ice specific stability functions -stability_functions = NumericalEarth.EarthSystemModels.atmosphere_sea_ice_stability_functions(Float64) - -# Edson et al. (2013) stability functions -stability_functions = NumericalEarth.EarthSystemModels.atmosphere_ocean_stability_functions(Float64) - -atmosphere_ocean_fluxes = SimilarityTheoryFluxes(; stability_functions) -interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(nothing, ocean; atmosphere_ocean_fluxes) -model = OceanOnlyModel(ocean; interfaces) - -# output -EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) -├── atmosphere: Nothing -├── land: Nothing -├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} -├── ocean: HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -└── interfaces: ComponentInterfaces -``` - -The available stability function options include: -- `atmosphere_ocean_stability_functions`: Based on [Edson et al. (2013)](@cite edson2013exchange) -- `atmosphere_sea_ice_stability_functions`: Specifically designed for atmosphere-sea ice interactions -- `nothing`: No stability functions will be used -- Custom stability functions can be created by defining functions of the "stability parameter" - (the flux Richardson number), `ζ`. """ -function EarthSystemModel(atmosphere, ocean, sea_ice; - land = nothing, - radiation = Radiation(), +function EarthSystemModel(radiation, atmosphere, land, sea_ice, ocean; clock = Clock{Float64}(time=0), ocean_reference_density = reference_density(ocean), ocean_heat_capacity = heat_capacity(ocean), sea_ice_reference_density = reference_density(sea_ice), sea_ice_heat_capacity = heat_capacity(sea_ice), interfaces = nothing) + if isnothing(radiation) && atmosphere isa PrescribedAtmosphere + @warn """ + `EarthSystemModel` was constructed with a `PrescribedAtmosphere` but \ + `radiation = nothing`. This means no upwelling longwave (ϵσT⁴), no \ + absorbed downwelling longwave, and no shortwave absorption — \ + results will be physically inconsistent. + + If you previously relied on `Radiation()` defaults: pass \ + `radiation = JRA55PrescribedRadiation(arch; backend, ...)` (or \ + `ECCOPrescribedRadiation` / `OSPapaPrescribedRadiation`) to restore \ + radiative forcing. Pass `radiation = PrescribedRadiation(grid)` for \ + emission-only mode. To suppress this warning, build the model \ + without a `PrescribedAtmosphere` (radiatively decoupled is the \ + intended `nothing` semantics). + """ maxlog=1 + end if ocean isa Simulation if !isnothing(ocean.callbacks) @@ -193,23 +177,28 @@ function EarthSystemModel(atmosphere, ocean, sea_ice; # Contains information about flux contributions: bulk formula, prescribed fluxes, etc. if isnothing(interfaces) && !(isnothing(atmosphere) && isnothing(sea_ice)) interfaces = ComponentInterfaces(atmosphere, ocean, sea_ice; + radiation, land, ocean_reference_density, ocean_heat_capacity, sea_ice_reference_density, - sea_ice_heat_capacity, - radiation) + sea_ice_heat_capacity) end arch = architecture(interfaces.exchanger.grid) + # Allocate per-surface InterfaceRadiationFlux on the exchange grid. + surfaces = present_surfaces(ocean, sea_ice) + allocate_interface_fluxes!(radiation, interfaces.exchanger.grid, surfaces) + earth_system_model = EarthSystemModel(arch, - clock, - atmosphere, - land, - sea_ice, - ocean, - interfaces) + clock, + radiation, + atmosphere, + land, + sea_ice, + ocean, + interfaces) # Make sure the initial temperature of the ocean # is not below freezing and above melting near the surface @@ -220,7 +209,42 @@ function EarthSystemModel(atmosphere, ocean, sea_ice; end """ - OceanOnlyModel(ocean; atmosphere=nothing, radiation=Radiation(), kw...) + EarthSystemModel(; radiation = nothing, + atmosphere = nothing, + land = nothing, + sea_ice = nothing, + ocean = nothing, + kw...) + +Keyword-only constructor for `EarthSystemModel`. Equivalent to the positional +form, but lets you pass only the components you actually have: + +```julia +EarthSystemModel(; atmosphere, ocean) # ocean + atmosphere +EarthSystemModel(; atmosphere, sea_ice, ocean, radiation) # full coupled +``` + +All keyword arguments accepted by the positional constructor are forwarded. +""" +EarthSystemModel(; radiation = nothing, + atmosphere = nothing, + land = nothing, + sea_ice = nothing, + ocean = nothing, + kw...) = + EarthSystemModel(radiation, atmosphere, land, sea_ice, ocean; kw...) + +# Determine which surfaces are present in the model — used to allocate +# per-surface diagnostic radiation flux buffers. +function present_surfaces(ocean, sea_ice) + surfaces = Symbol[] + isnothing(ocean) || push!(surfaces, :ocean) + isnothing(sea_ice) || push!(surfaces, :sea_ice) + return Tuple(surfaces) +end + +""" + OceanOnlyModel(ocean; atmosphere=nothing, radiation=nothing, land=nothing, kw...) Construct an ocean-only model without a sea ice component. This is a convenience constructor for [`EarthSystemModel`](@ref) that sets `sea_ice` @@ -229,58 +253,20 @@ to `FreezingLimitedOceanTemperature` (a simple freezing limiter that does not ev The `atmosphere` keyword can be used to specify a prescribed atmospheric forcing (e.g., `JRA55PrescribedAtmosphere`). All other keyword arguments are forwarded to `EarthSystemModel`. - -```jldoctest -using NumericalEarth -using Oceananigans - -grid = RectilinearGrid(size=10, z=(-100, 0), topology=(Flat, Flat, Bounded)) -ocean = ocean_simulation(grid, timestepper = :QuasiAdamsBashforth2) -model = OceanOnlyModel(ocean) - -# output -EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) -├── atmosphere: Nothing -├── land: Nothing -├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} -├── ocean: HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -└── interfaces: ComponentInterfaces -``` """ -OceanOnlyModel(ocean; atmosphere=nothing, land=nothing, kw...) = - EarthSystemModel(atmosphere, ocean, default_sea_ice(); land, kw...) +OceanOnlyModel(ocean; atmosphere=nothing, land=nothing, radiation=nothing, kw...) = + EarthSystemModel(radiation, atmosphere, land, default_sea_ice(), ocean; kw...) """ - OceanSeaIceModel(ocean, sea_ice; atmosphere=nothing, radiation=Radiation(), kw...) + OceanSeaIceModel(sea_ice, ocean; atmosphere=nothing, radiation=nothing, land=nothing, kw...) Construct a coupled ocean--sea ice model. This is a convenience constructor for [`EarthSystemModel`](@ref) with an explicit sea ice component -and an optional prescribed atmosphere. - -The `atmosphere` keyword can be used to specify a prescribed atmospheric forcing -(e.g., `JRA55PrescribedAtmosphere`). All other keyword arguments are forwarded -to `EarthSystemModel`. - -```jldoctest -using NumericalEarth -using Oceananigans - -grid = RectilinearGrid(size=10, z=(-100, 0), topology=(Flat, Flat, Bounded)) -ocean = ocean_simulation(grid, timestepper = :QuasiAdamsBashforth2) -sea_ice = FreezingLimitedOceanTemperature() -model = OceanSeaIceModel(ocean, sea_ice) - -# output -EarthSystemModel{CPU}(time = 0 seconds, iteration = 0) -├── atmosphere: Nothing -├── land: Nothing -├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}} -├── ocean: HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -└── interfaces: ComponentInterfaces -``` +and an optional prescribed atmosphere. Positional arguments follow the +struct convention (top→bottom): `sea_ice` then `ocean`. """ -OceanSeaIceModel(ocean, sea_ice; atmosphere=nothing, land=nothing, kw...) = - EarthSystemModel(atmosphere, ocean, sea_ice; land, kw...) +OceanSeaIceModel(sea_ice, ocean; atmosphere=nothing, land=nothing, radiation=nothing, kw...) = + EarthSystemModel(radiation, atmosphere, land, sea_ice, ocean; kw...) """ AtmosphereOceanModel(atmosphere, ocean; kw...) @@ -289,8 +275,8 @@ Construct a coupled atmosphere--ocean model. Convenience constructor for [`EarthSystemModel`](@ref) with an atmosphere and ocean but no sea ice. All keyword arguments are forwarded to `EarthSystemModel`. """ -AtmosphereOceanModel(atmosphere, ocean; land=nothing, kw...) = - EarthSystemModel(atmosphere, ocean, nothing; land, kw...) +AtmosphereOceanModel(atmosphere, ocean; land=nothing, radiation=nothing, kw...) = + EarthSystemModel(radiation, atmosphere, land, nothing, ocean; kw...) time(coupled_model::EarthSystemModel) = coupled_model.clock.time @@ -334,6 +320,7 @@ above_freezing_ocean_temperature!(ocean, grid, ::Nothing) = nothing function prognostic_state(osm::EarthSystemModel) return (clock = prognostic_state(osm.clock), + radiation = prognostic_state(osm.radiation), ocean = prognostic_state(osm.ocean), atmosphere = prognostic_state(osm.atmosphere), land = prognostic_state(osm.land), @@ -343,6 +330,10 @@ end function restore_prognostic_state!(osm::EarthSystemModel, state) restore_prognostic_state!(osm.clock, state.clock) + # Backwards-compatible: older checkpoints may not have a `radiation` entry + if hasproperty(state, :radiation) + restore_prognostic_state!(osm.radiation, state.radiation) + end restore_prognostic_state!(osm.ocean, state.ocean) restore_prognostic_state!(osm.atmosphere, state.atmosphere) restore_prognostic_state!(osm.land, state.land) diff --git a/src/EarthSystemModels/time_step_earth_system_model.jl b/src/EarthSystemModels/time_step_earth_system_model.jl index bf041418e..acd174ae8 100644 --- a/src/EarthSystemModels/time_step_earth_system_model.jl +++ b/src/EarthSystemModels/time_step_earth_system_model.jl @@ -7,26 +7,26 @@ using ClimaSeaIce: SeaIceModel, SeaIceThermodynamics using Oceananigans.Grids: φnode using Printf +# Hooks called from `update_state!` to apply radiative contributions on top of +# turbulent fluxes. Concrete radiation types overload these (no-op when +# `coupled_model.radiation === nothing`). +apply_air_sea_radiative_fluxes!(::Any) = nothing +apply_air_sea_ice_radiative_fluxes!(::Any) = nothing + function time_step!(coupled_model::EarthSystemModel, Δt; callbacks=[]) maybe_prepare_first_time_step!(coupled_model, callbacks) - - ocean = coupled_model.ocean - sea_ice = coupled_model.sea_ice - atmosphere = coupled_model.atmosphere - land = coupled_model.land - - # Eventually, split out into OceanOnlyModel - !isnothing(sea_ice) && time_step!(sea_ice, Δt) - # TODO after ice time-step: - # - Adjust ocean heat flux if the ice completely melts? - !isnothing(ocean) && time_step!(ocean, Δt) + radiation = coupled_model.radiation + atmosphere = coupled_model.atmosphere + land = coupled_model.land + sea_ice = coupled_model.sea_ice + ocean = coupled_model.ocean - # Time step the atmosphere + !isnothing(radiation) && time_step!(radiation, Δt) !isnothing(atmosphere) && time_step!(atmosphere, Δt) - - # Time step land - !isnothing(land) && time_step!(land, Δt) + !isnothing(land) && time_step!(land, Δt) + !isnothing(sea_ice) && time_step!(sea_ice, Δt) + !isnothing(ocean) && time_step!(ocean, Δt) # TODO: # - Store fractional ice-free / ice-covered _time_ for more @@ -39,31 +39,37 @@ end function update_state!(coupled_model::EarthSystemModel, callbacks=[]) - # The four components - ocean = coupled_model.ocean - sea_ice = coupled_model.sea_ice + radiation = coupled_model.radiation atmosphere = coupled_model.atmosphere land = coupled_model.land + sea_ice = coupled_model.sea_ice + ocean = coupled_model.ocean exchanger = coupled_model.interfaces.exchanger grid = exchanger.grid - # This function needs to be specialized to allow different component models + # Phase 1: bring all component states onto the exchange grid + interpolate_state!(exchanger.radiation, grid, radiation, coupled_model) interpolate_state!(exchanger.atmosphere, grid, atmosphere, coupled_model) interpolate_state!(exchanger.land, grid, land, coupled_model) - interpolate_state!(exchanger.ocean, grid, ocean, coupled_model) interpolate_state!(exchanger.sea_ice, grid, sea_ice, coupled_model) + interpolate_state!(exchanger.ocean, grid, ocean, coupled_model) - # Compute interface states + # Phase 2: compute interface turbulent fluxes compute_atmosphere_ocean_fluxes!(coupled_model) compute_atmosphere_sea_ice_fluxes!(coupled_model) compute_sea_ice_ocean_fluxes!(coupled_model) - # This function needs to be specialized to allow different component models + # Phase 3: assemble net component fluxes (turbulent only) + update_net_fluxes!(coupled_model, radiation) update_net_fluxes!(coupled_model, atmosphere) update_net_fluxes!(coupled_model, land) - update_net_fluxes!(coupled_model, ocean) update_net_fluxes!(coupled_model, sea_ice) + update_net_fluxes!(coupled_model, ocean) + + # Phase 4: add radiative contributions on top + apply_air_sea_radiative_fluxes!(coupled_model) + apply_air_sea_ice_radiative_fluxes!(coupled_model) return nothing end diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index fc5731492..0e8c5873e 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -15,8 +15,14 @@ export SlabOcean, default_sea_ice, FreezingLimitedOceanTemperature, - Radiation, + PrescribedRadiation, + SurfaceRadiationProperties, + InterfaceRadiationFlux, LatitudeDependentAlbedo, + TabulatedAlbedo, + JRA55PrescribedRadiation, + ECCOPrescribedRadiation, + OSPapaPrescribedRadiation, SimilarityTheoryFluxes, CoefficientBasedFluxes, MomentumRoughnessLength, @@ -111,6 +117,7 @@ include("EarthSystemModels/EarthSystemModels.jl") include("Oceans/Oceans.jl") include("Atmospheres/Atmospheres.jl") include("Lands/Lands.jl") +include("Radiations/Radiations.jl") include("SeaIces/SeaIces.jl") include("InitialConditions/InitialConditions.jl") include("DataWrangling/DataWrangling.jl") @@ -124,6 +131,7 @@ using .InitialConditions using .EarthSystemModels using .Atmospheres using .Lands +using .Radiations using .Oceans using .SeaIces using .Diagnostics diff --git a/src/Oceans/assemble_net_ocean_fluxes.jl b/src/Oceans/assemble_net_ocean_fluxes.jl index 5fe6f81fe..6c364a434 100644 --- a/src/Oceans/assemble_net_ocean_fluxes.jl +++ b/src/Oceans/assemble_net_ocean_fluxes.jl @@ -6,16 +6,13 @@ using NumericalEarth.EarthSystemModels: EarthSystemModel, NoOceanInterfaceModel, using NumericalEarth.EarthSystemModels.InterfaceComputations: interface_kernel_parameters, computed_fluxes, - sea_ice_concentration, - convert_to_kelvin, - emitted_longwave_radiation, - absorbed_longwave_radiation, - transmitted_shortwave_radiation + sea_ice_concentration @inline τᶜᶜᶜ(i, j, k, grid, ρᵒᶜ⁻¹, ℵ, ρτᶜᶜᶜ) = @inbounds ρᵒᶜ⁻¹ * (1 - ℵ[i, j, k]) * ρτᶜᶜᶜ[i, j, k] ##### -##### Generic flux assembler +##### Generic flux assembler — turbulent + sea-ice contributions only. +##### Radiative contributions are added later by `apply_air_sea_radiative_fluxes!`. ##### # Fallback for an ocean-only model (it has no interfaces!) @@ -24,7 +21,6 @@ update_net_fluxes!(coupled_model::Union{NoOceanInterfaceModel, NoInterfaceModel} update_net_fluxes!(coupled_model, ocean::Simulation{<:HydrostaticFreeSurfaceModel}) = update_net_ocean_fluxes!(coupled_model, ocean, ocean.model.grid) -# A generic ocean flux assembler for a coupled model with both an atmosphere and sea ice function update_net_ocean_fluxes!(coupled_model, ocean_model, grid) sea_ice = coupled_model.sea_ice arch = architecture(grid) @@ -34,13 +30,7 @@ function update_net_ocean_fluxes!(coupled_model, ocean_model, grid) atmos_ocean_fluxes = computed_fluxes(coupled_model.interfaces.atmosphere_ocean_interface) sea_ice_ocean_fluxes = computed_fluxes(coupled_model.interfaces.sea_ice_ocean_interface) - # Simplify NamedTuple to reduce parameter space consumption. - # See https://github.com/CliMA/NumericalEarth.jl/issues/116. atmosphere_fields = coupled_model.interfaces.exchanger.atmosphere.state - - downwelling_radiation = (ℐꜜˢʷ = atmosphere_fields.ℐꜜˢʷ.data, - ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data) - freshwater_flux = atmosphere_fields.Jᶜ.data # Extract land freshwater flux if land component is present @@ -49,27 +39,19 @@ function update_net_ocean_fluxes!(coupled_model, ocean_model, grid) ice_concentration = sea_ice_concentration(sea_ice) ocean_surface_salinity = EarthSystemModels.ocean_surface_salinity(ocean_model) - atmos_ocean_properties = coupled_model.interfaces.atmosphere_ocean_interface.properties ocean_properties = coupled_model.interfaces.ocean_properties - ocean_surface_temperature = coupled_model.interfaces.atmosphere_ocean_interface.temperature - penetrating_radiation = get_radiative_forcing(ocean_model) - launch!(arch, grid, :xy, _assemble_net_ocean_fluxes!, net_ocean_fluxes, - penetrating_radiation, grid, clock, atmos_ocean_fluxes, sea_ice_ocean_fluxes, ocean_surface_salinity, - ocean_surface_temperature, ice_concentration, - downwelling_radiation, freshwater_flux, land_freshwater_flux, - atmos_ocean_properties, ocean_properties) return nothing @@ -79,23 +61,18 @@ end Base.@propagate_inbounds get_land_freshwater_flux(i, j, flux) = flux[i, j, 1] @kernel function _assemble_net_ocean_fluxes!(net_ocean_fluxes, - penetrating_radiation, grid, clock, atmos_ocean_fluxes, sea_ice_ocean_fluxes, ocean_surface_salinity, - ocean_surface_temperature, sea_ice_concentration, - downwelling_radiation, freshwater_flux, land_freshwater_flux, - atmos_ocean_properties, ocean_properties) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) - time = Time(clock.time) ρτˣᵃᵒ = atmos_ocean_fluxes.x_momentum # atmosphere - ocean zonal momentum flux ρτʸᵃᵒ = atmos_ocean_fluxes.y_momentum # atmosphere - ocean meridional momentum flux ρτˣⁱᵒ = sea_ice_ocean_fluxes.x_momentum # sea_ice - ocean zonal momentum flux @@ -104,54 +81,24 @@ Base.@propagate_inbounds get_land_freshwater_flux(i, j, flux) = flux[i, j, 1] @inbounds begin ℵᵢ = sea_ice_concentration[i, j, 1] Sᵒᶜ = ocean_surface_salinity[i, j, 1] - Tₛ = ocean_surface_temperature[i, j, 1] - Tₛ = convert_to_kelvin(ocean_properties.temperature_units, Tₛ) - - Jᶜ = freshwater_flux[i, j, 1] + get_land_freshwater_flux(i, j, land_freshwater_flux) # Prescribed freshwater flux (atmos + land) - ℐꜜˢʷ = downwelling_radiation.ℐꜜˢʷ[i, j, 1] # Downwelling shortwave radiation - ℐꜜˡʷ = downwelling_radiation.ℐꜜˡʷ[i, j, 1] # Downwelling longwave radiation - 𝒬ᵀ = atmos_ocean_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux - 𝒬ᵛ = atmos_ocean_fluxes.latent_heat[i, j, 1] # latent heat flux - Jᵛ = atmos_ocean_fluxes.water_vapor[i, j, 1] # mass flux of water vapor - end - - # Compute radiation fluxes (radiation is multiplied by the fraction of ocean, 1 - sea ice concentration) - σ = atmos_ocean_properties.radiation.σ - α = atmos_ocean_properties.radiation.α - ϵ = atmos_ocean_properties.radiation.ϵ - ℐꜛˡʷ = emitted_longwave_radiation(i, j, kᴺ, grid, time, Tₛ, σ, ϵ) - ℐₐˡʷ = absorbed_longwave_radiation(i, j, kᴺ, grid, time, ϵ, ℐꜜˡʷ) - - # Compute the interior + surface absorbed shortwave radiation - ℐₜˢʷ = transmitted_shortwave_radiation(i, j, kᴺ, grid, time, α, ℐꜜˢʷ) - ℐₐˡʷ *= (1 - ℵᵢ) - ℐₜˢʷ *= (1 - ℵᵢ) - - Qss = shortwave_radiative_forcing(i, j, grid, penetrating_radiation, ℐₜˢʷ, ocean_properties) - - # Compute the total heat flux - ΣQao = (ℐꜛˡʷ + 𝒬ᵀ + 𝒬ᵛ) * (1 - ℵᵢ) + ℐₐˡʷ + Qss - - @inbounds begin - # Write radiative components of the heat flux for diagnostic purposes - atmos_ocean_fluxes.upwelling_longwave[i, j, 1] = ℐꜛˡʷ - atmos_ocean_fluxes.downwelling_longwave[i, j, 1] = - ℐₐˡʷ - atmos_ocean_fluxes.downwelling_shortwave[i, j, 1] = - ℐₜˢʷ + Jᶜ = freshwater_flux[i, j, 1] + get_land_freshwater_flux(i, j, land_freshwater_flux) + 𝒬ᵀ = atmos_ocean_fluxes.sensible_heat[i, j, 1] + 𝒬ᵛ = atmos_ocean_fluxes.latent_heat[i, j, 1] + Jᵛ = atmos_ocean_fluxes.water_vapor[i, j, 1] end - # Convert from a mass flux to a volume flux (aka velocity) - # by dividing with the ocean reference density. - # Also switch the sign, for some reason we are given freshwater flux as positive down. + # Turbulent contributions to surface heat flux (radiation added later) + ΣQao = (𝒬ᵀ + 𝒬ᵛ) * (1 - ℵᵢ) + + # Convert mass flux to volume flux; sign-flip (prescribed flux is positive down) ρᵒᶜ⁻¹ = 1 / ocean_properties.reference_density ΣFao = - Jᶜ * ρᵒᶜ⁻¹ - # Add the contribution from the turbulent water vapor flux, which has - # a different sign convention as the prescribed water mass fluxes (positive upwards) + # Add turbulent water vapor flux (positive upward sign convention) Jᵛᵒᶜ = Jᵛ * ρᵒᶜ⁻¹ ΣFao += Jᵛᵒᶜ - # Compute fluxes for u, v, T, and S from momentum, heat, and freshwater fluxes τˣ = net_ocean_fluxes.u τʸ = net_ocean_fluxes.v Jᵀ = net_ocean_fluxes.T @@ -166,7 +113,7 @@ Base.@propagate_inbounds get_land_freshwater_flux(i, j, flux) = flux[i, j, 1] Jᵀao = ΣQao * ρᵒᶜ⁻¹ * cᵒᶜ⁻¹ Jᵀio = 𝒬ⁱⁿᵗ * ρᵒᶜ⁻¹ * cᵒᶜ⁻¹ - # salinity flux > 0 extracts salinity from the ocean --- the opposite of a water vapor flux + # salinity flux > 0 extracts salinity (opposite of water vapor flux sign) Jˢao = - Sᵒᶜ * ΣFao τˣᵃᵒ = ℑxᶠᵃᵃ(i, j, 1, grid, τᶜᶜᶜ, ρᵒᶜ⁻¹, ℵ, ρτˣᵃᵒ) @@ -174,12 +121,11 @@ Base.@propagate_inbounds get_land_freshwater_flux(i, j, flux) = flux[i, j, 1] τˣⁱᵒ = ρτˣⁱᵒ[i, j, 1] * ρᵒᶜ⁻¹ * ℑxᶠᵃᵃ(i, j, 1, grid, ℵ) τʸⁱᵒ = ρτʸⁱᵒ[i, j, 1] * ρᵒᶜ⁻¹ * ℑyᵃᶠᵃ(i, j, 1, grid, ℵ) - # Stresses τˣ[i, j, 1] = ifelse(inactive, zero(grid), τˣᵃᵒ + τˣⁱᵒ) τʸ[i, j, 1] = ifelse(inactive, zero(grid), τʸᵃᵒ + τʸⁱᵒ) - # Tracer fluxes - Jᵀ[i, j, 1] = ifelse(inactive, zero(grid), Jᵀao + Jᵀio) # Jᵀao is already multiplied by the sea ice concentration + # Tracer fluxes — radiative contributions added later by apply_air_sea_radiative_fluxes! + Jᵀ[i, j, 1] = ifelse(inactive, zero(grid), Jᵀao + Jᵀio) Jˢ[i, j, 1] = ifelse(inactive, zero(grid), (1 - ℵᵢ) * Jˢao + Jˢio) end end diff --git a/src/Radiations/Radiations.jl b/src/Radiations/Radiations.jl new file mode 100644 index 000000000..fd31f94c4 --- /dev/null +++ b/src/Radiations/Radiations.jl @@ -0,0 +1,55 @@ +module Radiations + +export PrescribedRadiation, + SurfaceRadiationProperties, + InterfaceRadiationFlux, + LatitudeDependentAlbedo, + TabulatedAlbedo, + default_stefan_boltzmann_constant + +# CODATA 2018 value of the Stefan–Boltzmann constant, in W m⁻² K⁻⁴. +const default_stefan_boltzmann_constant = 5.670374419e-8 + +using Oceananigans +using Oceananigans.Architectures: architecture +using Oceananigans.Fields: Center, Face, Field, ZeroField, FractionalIndices, instantiate, interpolator +using Oceananigans.Grids: grid_name, prettysummary, ηnode, _node, topology, Flat, on_architecture +using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series, cpu_interpolating_time_indices +using Oceananigans.TimeSteppers: Clock, tick! +using Oceananigans.Units: Time +using Oceananigans.Utils: launch! + +using Adapt +using KernelAbstractions: @kernel, @index + +import NumericalEarth: stateindex +using NumericalEarth.EarthSystemModels.InterfaceComputations: interface_kernel_parameters + +import Oceananigans.TimeSteppers: time_step!, update_state! +import Oceananigans.Architectures: on_architecture + +import NumericalEarth.EarthSystemModels: interpolate_state!, + update_net_fluxes!, + apply_air_sea_radiative_fluxes!, + apply_air_sea_ice_radiative_fluxes!, + allocate_interface_fluxes! + +import NumericalEarth.EarthSystemModels.InterfaceComputations: ComponentExchanger, + initialize!, + kernel_radiation_properties, + air_sea_interface_radiation_state, + air_sea_ice_interface_radiation_state + +include("surface_radiation_properties.jl") +include("interface_radiation_flux.jl") +include("radiation_kernels.jl") +include("latitude_dependent_albedo.jl") +include("tabulated_albedo.jl") +include("prescribed_radiation.jl") +include("prescribed_radiation_regridder.jl") +include("interpolate_radiation_state.jl") +include("air_sea_interface_radiation_state.jl") +include("apply_air_sea_radiative_fluxes.jl") +include("apply_air_sea_ice_radiative_fluxes.jl") + +end # module Radiations diff --git a/src/Radiations/air_sea_interface_radiation_state.jl b/src/Radiations/air_sea_interface_radiation_state.jl new file mode 100644 index 000000000..c7748fa8d --- /dev/null +++ b/src/Radiations/air_sea_interface_radiation_state.jl @@ -0,0 +1,26 @@ +# PrescribedRadiation-aware methods for the radiation getter functions +# declared (with `nothing` fallbacks) in InterfaceComputations. + +@inline kernel_radiation_properties(r::PrescribedRadiation) = + (σ = r.stefan_boltzmann_constant, + surface_properties = r.surface_properties) + +@inline function air_sea_interface_radiation_state(rk, exchanger_state, i, j, k, grid, time) + σ = rk.σ + @inbounds ℐꜜˢʷ = exchanger_state.ℐꜜˢʷ[i, j, 1] + @inbounds ℐꜜˡʷ = exchanger_state.ℐꜜˡʷ[i, j, 1] + s = rk.surface_properties.ocean + α = stateindex(s.albedo, i, j, k, grid, time, (Center, Center, Center), ℐꜜˢʷ) + ϵ = stateindex(s.emissivity, i, j, k, grid, time, (Center, Center, Center)) + return (; σ, α, ϵ, ℐꜜˢʷ, ℐꜜˡʷ) +end + +@inline function air_sea_ice_interface_radiation_state(rk, exchanger_state, i, j, k, grid, time) + σ = rk.σ + @inbounds ℐꜜˢʷ = exchanger_state.ℐꜜˢʷ[i, j, 1] + @inbounds ℐꜜˡʷ = exchanger_state.ℐꜜˡʷ[i, j, 1] + s = rk.surface_properties.sea_ice + α = stateindex(s.albedo, i, j, k, grid, time, (Center, Center, Center), ℐꜜˢʷ) + ϵ = stateindex(s.emissivity, i, j, k, grid, time, (Center, Center, Center)) + return (; σ, α, ϵ, ℐꜜˢʷ, ℐꜜˡʷ) +end diff --git a/src/Radiations/apply_air_sea_ice_radiative_fluxes.jl b/src/Radiations/apply_air_sea_ice_radiative_fluxes.jl new file mode 100644 index 000000000..eb2d9249d --- /dev/null +++ b/src/Radiations/apply_air_sea_ice_radiative_fluxes.jl @@ -0,0 +1,93 @@ +using ClimaSeaIce: SeaIceModel + +""" + apply_air_sea_ice_radiative_fluxes!(coupled_model) + +Add the radiative contribution to the net sea-ice top heat flux and write +the diagnostic radiative fluxes (upwelling LW, absorbed LW, transmitted SW) +into `coupled_model.radiation.interface_fluxes.sea_ice`. + +When `coupled_model.radiation === nothing`, this is a no-op. +Also a no-op when sea-ice is not a `Simulation{<:SeaIceModel}`. +""" +apply_air_sea_ice_radiative_fluxes!(::EarthSystemModel{<:Nothing}) = nothing + +apply_air_sea_ice_radiative_fluxes!(coupled_model::EarthSystemModel) = + _apply_air_sea_ice_radiative_fluxes_dispatch!(coupled_model, coupled_model.sea_ice) + +# No sea-ice or non-prognostic sea-ice: nothing to do. +_apply_air_sea_ice_radiative_fluxes_dispatch!(coupled_model, ::Any) = nothing + +function _apply_air_sea_ice_radiative_fluxes_dispatch!(coupled_model::EarthSystemModel, + sea_ice::Simulation{<:SeaIceModel}) + radiation = coupled_model.radiation + interface_fluxes = radiation.interface_fluxes + isnothing(interface_fluxes) && return nothing + haskey(interface_fluxes, :sea_ice) || return nothing + + grid = sea_ice.model.grid + arch = architecture(grid) + clock = coupled_model.clock + + top_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice.top.heat + radiation_state = coupled_model.interfaces.exchanger.radiation.state + rk = kernel_radiation_properties(radiation) + + ice_concentration = sea_ice_concentration(sea_ice) + surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature + sea_ice_properties = coupled_model.interfaces.sea_ice_properties + + launch!(arch, grid, :xy, + _apply_air_sea_ice_radiative_fluxes!, + top_heat_flux, + interface_fluxes.sea_ice, + grid, + clock, + rk, + radiation_state, + ice_concentration, + surface_temperature, + sea_ice_properties) + + return nothing +end + +@kernel function _apply_air_sea_ice_radiative_fluxes!(top_heat_flux, + interface_radiative_flux, + grid, + clock, + rk, + radiation_state, + ice_concentration, + surface_temperature, + sea_ice_properties) + + i, j = @index(Global, NTuple) + kᴺ = size(grid, 3) + time = Time(clock.time) + + @inbounds begin + ℵᵢ = ice_concentration[i, j, 1] + Ts = surface_temperature[i, j, kᴺ] + end + Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts) + + rs = air_sea_ice_interface_radiation_state(rk, radiation_state, i, j, kᴺ, grid, time) + + ℐꜛˡʷ = rs.σ * rs.ϵ * Ts^4 + ℐₐˡʷ = - rs.ϵ * rs.ℐꜜˡʷ + ℐₜˢʷ = - (1 - rs.α) * rs.ℐꜜˢʷ + + # Sea-ice radiation contributes only where ice exists. + ice_present = ℵᵢ > 0 + ΣQ_rad_ice = (ℐꜛˡʷ + ℐₐˡʷ + ℐₜˢʷ) * ice_present + + inactive = inactive_node(i, j, kᴺ, grid, Center(), Center(), Center()) + + @inbounds begin + top_heat_flux[i, j, 1] += ifelse(inactive, zero(grid), ΣQ_rad_ice) + interface_radiative_flux.upwelling_longwave[i, j, 1] = ℐꜛˡʷ + interface_radiative_flux.downwelling_longwave[i, j, 1] = - ℐₐˡʷ + interface_radiative_flux.downwelling_shortwave[i, j, 1] = - ℐₜˢʷ + end +end diff --git a/src/Radiations/apply_air_sea_radiative_fluxes.jl b/src/Radiations/apply_air_sea_radiative_fluxes.jl new file mode 100644 index 000000000..8441d2b63 --- /dev/null +++ b/src/Radiations/apply_air_sea_radiative_fluxes.jl @@ -0,0 +1,110 @@ +using Oceananigans.Grids: inactive_node +using NumericalEarth.EarthSystemModels: EarthSystemModel +using NumericalEarth.EarthSystemModels.InterfaceComputations: convert_to_kelvin, sea_ice_concentration +using NumericalEarth.Oceans: shortwave_radiative_forcing, get_radiative_forcing + +""" + apply_air_sea_radiative_fluxes!(coupled_model) + +Add the radiative contribution to the net ocean heat flux `Jᵀ` and write +the diagnostic radiative fluxes (upwelling LW, absorbed LW, transmitted SW) +into `coupled_model.radiation.interface_fluxes.ocean`. + +When `coupled_model.radiation === nothing`, this is a no-op. +""" +apply_air_sea_radiative_fluxes!(::EarthSystemModel{<:Nothing}) = nothing + +function apply_air_sea_radiative_fluxes!(coupled_model::EarthSystemModel) + ocean = coupled_model.ocean + isnothing(ocean) && return nothing + + # No atmosphere--ocean interface (no atmosphere or no ocean): nothing to do. + ao_interface = coupled_model.interfaces.atmosphere_ocean_interface + isnothing(ao_interface) && return nothing + + radiation = coupled_model.radiation + interface_fluxes = radiation.interface_fluxes + isnothing(interface_fluxes) && return nothing + haskey(interface_fluxes, :ocean) || return nothing + + grid = ocean.model.grid + arch = architecture(grid) + clock = coupled_model.clock + + net_ocean_fluxes = coupled_model.interfaces.net_fluxes.ocean + radiation_state = coupled_model.interfaces.exchanger.radiation.state + rk = kernel_radiation_properties(radiation) + + sea_ice = coupled_model.sea_ice + ice_concentration = sea_ice_concentration(sea_ice) + + interface_temperature = ao_interface.temperature + ocean_properties = coupled_model.interfaces.ocean_properties + penetrating_radiation = get_radiative_forcing(ocean.model) + + launch!(arch, grid, :xy, + _apply_air_sea_radiative_fluxes!, + net_ocean_fluxes, + interface_fluxes.ocean, + penetrating_radiation, + grid, + clock, + rk, + radiation_state, + ice_concentration, + interface_temperature, + ocean_properties) + + return nothing +end + +@kernel function _apply_air_sea_radiative_fluxes!(net_ocean_fluxes, + interface_radiative_flux, + penetrating_radiation, + grid, + clock, + rk, + radiation_state, + ice_concentration, + ocean_surface_temperature, + ocean_properties) + + i, j = @index(Global, NTuple) + kᴺ = size(grid, 3) + time = Time(clock.time) + + @inbounds begin + ℵᵢ = ice_concentration[i, j, 1] + Tₛ = ocean_surface_temperature[i, j, 1] + end + Tₛ = convert_to_kelvin(ocean_properties.temperature_units, Tₛ) + + rs = air_sea_interface_radiation_state(rk, radiation_state, i, j, kᴺ, grid, time) + + ℐꜛˡʷ = rs.σ * rs.ϵ * Tₛ^4 + ℐₐˡʷ = - rs.ϵ * rs.ℐꜜˡʷ + ℐₜˢʷ = - (1 - rs.α) * rs.ℐꜜˢʷ + + # Multiply by ocean fraction (only the parts not blocked by ice) + ℐₐˡʷ *= (1 - ℵᵢ) + ℐₜˢʷ *= (1 - ℵᵢ) + ℐꜛˡʷ_ocean = ℐꜛˡʷ * (1 - ℵᵢ) + + Qss = shortwave_radiative_forcing(i, j, grid, penetrating_radiation, ℐₜˢʷ, ocean_properties) + + # Total radiative contribution to surface heat flux + ΣQ_rad = ℐꜛˡʷ_ocean + ℐₐˡʷ + Qss + + ρᵒᶜ⁻¹ = 1 / ocean_properties.reference_density + cᵒᶜ⁻¹ = 1 / ocean_properties.heat_capacity + Jᵀ_rad = ΣQ_rad * ρᵒᶜ⁻¹ * cᵒᶜ⁻¹ + + inactive = inactive_node(i, j, kᴺ, grid, Center(), Center(), Center()) + + @inbounds begin + net_ocean_fluxes.T[i, j, 1] += ifelse(inactive, zero(grid), Jᵀ_rad) + interface_radiative_flux.upwelling_longwave[i, j, 1] = ℐꜛˡʷ + interface_radiative_flux.downwelling_longwave[i, j, 1] = - ℐₐˡʷ + interface_radiative_flux.downwelling_shortwave[i, j, 1] = - ℐₜˢʷ + end +end diff --git a/src/Radiations/interface_radiation_flux.jl b/src/Radiations/interface_radiation_flux.jl new file mode 100644 index 000000000..85bf1d376 --- /dev/null +++ b/src/Radiations/interface_radiation_flux.jl @@ -0,0 +1,34 @@ +""" + InterfaceRadiationFlux{F} + +Container for the diagnostic radiative fluxes at an air–surface interface. +The same struct type is instantiated per surface (ocean, sea ice, snow, ...). + +Fields +====== +- `upwelling_longwave :: F` ϵσT⁴ +- `downwelling_longwave :: F` ϵℐꜜˡʷ (absorbed by the surface) +- `downwelling_shortwave :: F` (1−α)ℐꜜˢʷ (transmitted into the surface) +""" +struct InterfaceRadiationFlux{F} + upwelling_longwave :: F + downwelling_longwave :: F + downwelling_shortwave :: F +end + +function InterfaceRadiationFlux(grid) + F = Field{Center, Center, Nothing} + return InterfaceRadiationFlux(F(grid), F(grid), F(grid)) +end + +InterfaceRadiationFlux(::Nothing) = InterfaceRadiationFlux(ntuple(_ -> ZeroField(), 3)...) + +Adapt.adapt_structure(to, fluxes::InterfaceRadiationFlux) = + InterfaceRadiationFlux(Adapt.adapt(to, fluxes.upwelling_longwave), + Adapt.adapt(to, fluxes.downwelling_longwave), + Adapt.adapt(to, fluxes.downwelling_shortwave)) + +on_architecture(arch, fluxes::InterfaceRadiationFlux) = + InterfaceRadiationFlux(on_architecture(arch, fluxes.upwelling_longwave), + on_architecture(arch, fluxes.downwelling_longwave), + on_architecture(arch, fluxes.downwelling_shortwave)) diff --git a/src/Radiations/interpolate_radiation_state.jl b/src/Radiations/interpolate_radiation_state.jl new file mode 100644 index 000000000..277149bb0 --- /dev/null +++ b/src/Radiations/interpolate_radiation_state.jl @@ -0,0 +1,69 @@ +using NumericalEarth.Atmospheres: interp_atmos_time_series + +"""Interpolate the prescribed downwelling radiation onto the exchange grid.""" +function interpolate_state!(exchanger, grid, radiation::PrescribedRadiation, coupled_model) + arch = architecture(grid) + clock = coupled_model.clock + + ℐꜜˢʷ = radiation.downwelling_shortwave + ℐꜜˡʷ = radiation.downwelling_longwave + downwelling = (shortwave = ℐꜜˢʷ.data, longwave = ℐꜜˡʷ.data) + + radiation_times = ℐꜜˢʷ.times + radiation_backend = ℐꜜˢʷ.backend + radiation_time_indexing = ℐꜜˢʷ.time_indexing + + space_fractional_indices = exchanger.regridder + state = exchanger.state + state_data = (ℐꜜˢʷ = state.ℐꜜˢʷ.data, + ℐꜜˡʷ = state.ℐꜜˡʷ.data) + + t = clock.time + time_interpolator = cpu_interpolating_time_indices(arch, radiation_times, + radiation_time_indexing, t) + + kernel_parameters = interface_kernel_parameters(grid) + + launch!(arch, grid, kernel_parameters, + _interpolate_radiation_state!, + state_data, + space_fractional_indices, + time_interpolator, + grid, + downwelling, + radiation_backend, + radiation_time_indexing) + + return nothing +end + +@inline get_fractional_index(i, j, ::Nothing) = nothing +@inline get_fractional_index(i, j, frac) = @inbounds frac[i, j, 1] + +@kernel function _interpolate_radiation_state!(state, + space_fractional_indices, + time_interpolator, + exchange_grid, + downwelling, + rad_backend, + rad_time_indexing) + + i, j = @index(Global, NTuple) + + ii = space_fractional_indices.i + jj = space_fractional_indices.j + fi = get_fractional_index(i, j, ii) + fj = get_fractional_index(i, j, jj) + + x_itp = FractionalIndices(fi, fj, nothing) + t_itp = time_interpolator + args = (x_itp, t_itp, rad_backend, rad_time_indexing) + + ℐꜜˢʷ = interp_atmos_time_series(downwelling.shortwave, args...) + ℐꜜˡʷ = interp_atmos_time_series(downwelling.longwave, args...) + + @inbounds begin + state.ℐꜜˢʷ[i, j, 1] = ℐꜜˢʷ + state.ℐꜜˡʷ[i, j, 1] = ℐꜜˡʷ + end +end diff --git a/src/EarthSystemModels/InterfaceComputations/latitude_dependent_albedo.jl b/src/Radiations/latitude_dependent_albedo.jl similarity index 100% rename from src/EarthSystemModels/InterfaceComputations/latitude_dependent_albedo.jl rename to src/Radiations/latitude_dependent_albedo.jl diff --git a/src/Radiations/prescribed_radiation.jl b/src/Radiations/prescribed_radiation.jl new file mode 100644 index 000000000..0c1a3332a --- /dev/null +++ b/src/Radiations/prescribed_radiation.jl @@ -0,0 +1,160 @@ +""" + PrescribedRadiation{G, T, FT, SW, LW, S, IF, TI} + +Top-level radiation component holding prescribed downwelling shortwave and +longwave radiation as `FieldTimeSeries`, plus per-surface radiative properties +(albedo, emissivity) and the Stefan–Boltzmann constant. Diagnostic radiative +fluxes (one `InterfaceRadiationFlux` per surface) are populated by the +`apply_air_sea_*_radiative_fluxes!` kernels at every step; `interface_fluxes` +is `nothing` until the radiation is paired with an `EarthSystemModel` (which +allocates the per-surface buffers on the exchange grid). +""" +mutable struct PrescribedRadiation{G, T, FT, SW, LW, S, TI} + grid :: G + clock :: Clock{T} + downwelling_shortwave :: SW + downwelling_longwave :: LW + surface_properties :: S + stefan_boltzmann_constant :: FT + # NamedTuple of `InterfaceRadiationFlux`, allocated at `EarthSystemModel` + # construction time once the exchange grid is known. Untyped so the field + # can be reassigned from `nothing` to a populated NamedTuple. + interface_fluxes + times :: TI +end + +function Base.summary(r::PrescribedRadiation) + Nx, Ny, Nz = size(r.grid) + Nt = length(r.times) + sz_str = string(Nx, "×", Ny, "×", Nz, "×", Nt) + return string(sz_str, " PrescribedRadiation") +end + +function Base.show(io::IO, r::PrescribedRadiation) + print(io, summary(r), " on ", grid_name(r.grid), ":", '\n') + print(io, "├── times: ", prettysummary(r.times), '\n') + print(io, "├── stefan_boltzmann_constant: ", prettysummary(r.stefan_boltzmann_constant), '\n') + print(io, "└── surface_properties: ", keys(r.surface_properties)) +end + +# Filter out surfaces whose property kwarg was passed as `nothing`. +@inline function _filter_surface_properties(; ocean=nothing, sea_ice=nothing, snow=nothing, land=nothing) + pairs = () + isnothing(ocean) || (pairs = (pairs..., :ocean => ocean)) + isnothing(sea_ice) || (pairs = (pairs..., :sea_ice => sea_ice)) + isnothing(snow) || (pairs = (pairs..., :snow => snow)) + isnothing(land) || (pairs = (pairs..., :land => land)) + return NamedTuple(pairs) +end + +""" + PrescribedRadiation(downwelling_shortwave, downwelling_longwave; + clock = nothing, + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0), + snow_surface = nothing, + land_surface = nothing, + stefan_boltzmann_constant = default_stefan_boltzmann_constant) + +Construct a `PrescribedRadiation` component from `FieldTimeSeries` of +downwelling shortwave and longwave radiation. Grid + times are inferred from +the shortwave FTS. + +Pass `*_surface = nothing` to omit that surface from `surface_properties`. +""" +function PrescribedRadiation(downwelling_shortwave::FieldTimeSeries, + downwelling_longwave::FieldTimeSeries; + clock = nothing, + ocean_surface = SurfaceRadiationProperties(0.05, 0.97), + sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0), + snow_surface = nothing, + land_surface = nothing, + stefan_boltzmann_constant = default_stefan_boltzmann_constant) + + grid = downwelling_shortwave.grid + times = downwelling_shortwave.times + FT = eltype(downwelling_shortwave) + + if isnothing(clock) + clock = Clock{FT}(time = 0) + end + + surface_properties = _filter_surface_properties(ocean = ocean_surface, + sea_ice = sea_ice_surface, + snow = snow_surface, + land = land_surface) + + radiation = PrescribedRadiation(grid, + clock, + downwelling_shortwave, + downwelling_longwave, + surface_properties, + convert(FT, stefan_boltzmann_constant), + nothing, # interface_fluxes — populated at ESM construction + times) + update_state!(radiation) + return radiation +end + +""" + PrescribedRadiation(grid, times = [zero(grid)]; kwargs...) + +Construct a `PrescribedRadiation` with zero downwelling shortwave and longwave +fields on `grid`. Useful for emission-only mode (the surface radiates ϵσT⁴ but +no incoming radiation is absorbed). All other keyword arguments are forwarded +to the FTS-form constructor. +""" +function PrescribedRadiation(grid, times = [zero(eltype(grid))]; kwargs...) + sw = FieldTimeSeries{Center, Center, Nothing}(grid, times) + lw = FieldTimeSeries{Center, Center, Nothing}(grid, times) + return PrescribedRadiation(sw, lw; kwargs...) +end + +@inline function update_state!(radiation::PrescribedRadiation) + time = Time(radiation.clock.time) + ftses = extract_field_time_series(radiation) + + for fts in ftses + update_field_time_series!(fts, time) + end + return nothing +end + +@inline function time_step!(radiation::PrescribedRadiation, Δt) + tick!(radiation.clock, Δt) + update_state!(radiation) + return nothing +end + +# Prescribed radiation has no internal state to update from net fluxes. +update_net_fluxes!(coupled_model, ::PrescribedRadiation) = nothing + +""" + allocate_interface_fluxes!(radiation, exchange_grid, surfaces) + +Populate `radiation.interface_fluxes` with one `InterfaceRadiationFlux` +per surface present in the model. +""" +function allocate_interface_fluxes!(radiation::PrescribedRadiation, exchange_grid, surfaces) + pairs = (surface => InterfaceRadiationFlux(exchange_grid) for surface in surfaces) + radiation.interface_fluxes = NamedTuple(pairs) + return nothing +end + +##### +##### Checkpointing +##### + +import Oceananigans: prognostic_state, restore_prognostic_state! + +function prognostic_state(radiation::PrescribedRadiation) + return (; clock = prognostic_state(radiation.clock)) +end + +function restore_prognostic_state!(radiation::PrescribedRadiation, state) + restore_prognostic_state!(radiation.clock, state.clock) + update_state!(radiation) + return radiation +end + +restore_prognostic_state!(radiation::PrescribedRadiation, ::Nothing) = radiation diff --git a/src/Radiations/prescribed_radiation_regridder.jl b/src/Radiations/prescribed_radiation_regridder.jl new file mode 100644 index 000000000..5316fb5eb --- /dev/null +++ b/src/Radiations/prescribed_radiation_regridder.jl @@ -0,0 +1,51 @@ +function ComponentExchanger(radiation::PrescribedRadiation, grid) + + regridder = radiation_regridder(radiation, grid) + + state = (; ℐꜜˢʷ = Field{Center, Center, Nothing}(grid), + ℐꜜˡʷ = Field{Center, Center, Nothing}(grid)) + + return ComponentExchanger(state, regridder) +end + +function radiation_regridder(radiation::PrescribedRadiation, exchange_grid) + rad_grid = radiation.grid + arch = architecture(exchange_grid) + + FT = eltype(rad_grid) + TX, TY, TZ = topology(exchange_grid) + fi = TX() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FT) + fj = TY() isa Flat ? nothing : Field{Center, Center, Nothing}(exchange_grid, FT) + return (i = fi, j = fj) +end + +function initialize!(exchanger::ComponentExchanger, grid, radiation::PrescribedRadiation) + frac_indices = exchanger.regridder + # Skip horizontal regridding when both fractional-index buffers are + # absent (purely Flat horizontal exchange grid). + if isnothing(frac_indices.i) && isnothing(frac_indices.j) + return nothing + end + rad_grid = radiation.grid + kernel_parameters = interface_kernel_parameters(grid) + launch!(architecture(grid), grid, kernel_parameters, + _compute_radiation_fractional_indices!, frac_indices, grid, rad_grid) + return nothing +end + +@kernel function _compute_radiation_fractional_indices!(indices_tuple, exchange_grid, rad_grid) + i, j = @index(Global, NTuple) + kᴺ = size(exchange_grid, 3) + X = _node(i, j, kᴺ + 1, exchange_grid, Center(), Center(), Face()) + fractional_indices_ij = FractionalIndices(X, rad_grid, Center(), Center(), Center()) + fi = indices_tuple.i + fj = indices_tuple.j + @inbounds begin + if !isnothing(fi) + fi[i, j, 1] = fractional_indices_ij.i + end + if !isnothing(fj) + fj[i, j, 1] = fractional_indices_ij.j + end + end +end diff --git a/src/Radiations/radiation_kernels.jl b/src/Radiations/radiation_kernels.jl new file mode 100644 index 000000000..439d8277e --- /dev/null +++ b/src/Radiations/radiation_kernels.jl @@ -0,0 +1,24 @@ +@inline hack_cosd(φ) = cos(π * φ / 180) +@inline hack_sind(φ) = sin(π * φ / 180) + +const CCC = (Center, Center, Center) + +@inline function emitted_longwave_radiation(i, j, k, grid, time, T, σ, ϵ) + ϵi = stateindex(ϵ, i, j, k, grid, time, CCC) + return σ * ϵi * T^4 +end + +@inline function absorbed_longwave_radiation(i, j, k, grid, time, ϵ, ℐꜜˡʷ) + ϵi = stateindex(ϵ, i, j, k, grid, time, CCC) + return - ϵi * ℐꜜˡʷ +end + +@inline function transmitted_shortwave_radiation(i, j, k, grid, time, α, ℐꜜˢʷ) + αi = stateindex(α, i, j, k, grid, time, CCC, ℐꜜˢʷ) + return - (1 - αi) * ℐꜜˢʷ +end + +# Inside the solver we lose both spatial and temporal information, but the +# radiative properties have already been computed correctly +@inline net_absorbed_interface_radiation(ℐꜜˢʷ, ℐꜜˡʷ, α, ϵ) = - (1 - α) * ℐꜜˢʷ - ϵ * ℐꜜˡʷ +@inline emitted_longwave_radiation(T, σ, ϵ) = σ * ϵ * T^4 diff --git a/src/Radiations/surface_radiation_properties.jl b/src/Radiations/surface_radiation_properties.jl new file mode 100644 index 000000000..0b069189d --- /dev/null +++ b/src/Radiations/surface_radiation_properties.jl @@ -0,0 +1,29 @@ +struct SurfaceRadiationProperties{A, E} + albedo :: A + emissivity :: E +end + +""" + SurfaceRadiationProperties(albedo, emissivity) + +Bundle the radiative properties of a single surface (ocean, sea ice, snow, land) +that participate in radiative flux computation: shortwave reflectivity (`albedo`) +and longwave emissivity (`emissivity`). + +`albedo` may be a `Number`, a `LatitudeDependentAlbedo`, a `TabulatedAlbedo`, or +any other object for which `stateindex` is defined. `emissivity` may be a `Number` +or any other `stateindex`-able object. +""" +SurfaceRadiationProperties(; albedo, emissivity = 0.97) = SurfaceRadiationProperties(albedo, emissivity) + +Adapt.adapt_structure(to, s::SurfaceRadiationProperties) = + SurfaceRadiationProperties(Adapt.adapt(to, s.albedo), + Adapt.adapt(to, s.emissivity)) + +Base.summary(::SurfaceRadiationProperties) = "SurfaceRadiationProperties" + +function Base.show(io::IO, s::SurfaceRadiationProperties) + print(io, summary(s), ":", '\n') + print(io, "├── albedo: ", prettysummary(s.albedo), '\n') + print(io, "└── emissivity: ", prettysummary(s.emissivity)) +end diff --git a/src/EarthSystemModels/InterfaceComputations/tabulated_albedo.jl b/src/Radiations/tabulated_albedo.jl similarity index 100% rename from src/EarthSystemModels/InterfaceComputations/tabulated_albedo.jl rename to src/Radiations/tabulated_albedo.jl diff --git a/src/SeaIces/assemble_net_sea_ice_fluxes.jl b/src/SeaIces/assemble_net_sea_ice_fluxes.jl index 17d1e10c0..0021cbffb 100644 --- a/src/SeaIces/assemble_net_sea_ice_fluxes.jl +++ b/src/SeaIces/assemble_net_sea_ice_fluxes.jl @@ -1,9 +1,6 @@ using NumericalEarth.EarthSystemModels.InterfaceComputations: computed_fluxes, interface_kernel_parameters, - convert_to_kelvin, - emitted_longwave_radiation, - absorbed_longwave_radiation, - transmitted_shortwave_radiation + convert_to_kelvin update_net_fluxes!(coupled_model, ::FreezingLimitedOceanTemperature) = nothing @@ -18,21 +15,12 @@ function update_net_fluxes!(coupled_model, sea_ice::Simulation{<:SeaIceModel}) sea_ice_ocean_fluxes = computed_fluxes(coupled_model.interfaces.sea_ice_ocean_interface) atmosphere_sea_ice_fluxes = computed_fluxes(coupled_model.interfaces.atmosphere_sea_ice_interface) - # Simplify NamedTuple to reduce parameter space consumption. - # See https://github.com/CliMA/NumericalEarth.jl/issues/116. atmosphere_fields = coupled_model.interfaces.exchanger.atmosphere.state - - downwelling_radiation = (ℐꜜˢʷ = atmosphere_fields.ℐꜜˢʷ.data, - ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data) - freshwater_flux = atmosphere_fields.Jᶜ.data - atmos_sea_ice_properties = coupled_model.interfaces.atmosphere_sea_ice_interface.properties sea_ice_properties = coupled_model.interfaces.sea_ice_properties - - sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature ice_concentration = sea_ice_concentration(sea_ice) - + launch!(arch, grid, :xy, _assemble_net_sea_ice_fluxes!, top_fluxes, @@ -43,10 +31,7 @@ function update_net_fluxes!(coupled_model, sea_ice::Simulation{<:SeaIceModel}) sea_ice_ocean_fluxes, freshwater_flux, ice_concentration, - sea_ice_surface_temperature, - downwelling_radiation, - sea_ice_properties, - atmos_sea_ice_properties) + sea_ice_properties) return nothing end @@ -59,24 +44,15 @@ end sea_ice_ocean_fluxes, freshwater_flux, # Where do we add this one? ice_concentration, - surface_temperature, - downwelling_radiation, - sea_ice_properties, - atmos_sea_ice_properties) + sea_ice_properties) i, j = @index(Global, NTuple) kᴺ = size(grid, 3) - time = Time(clock.time) @inbounds begin - Ts = surface_temperature[i, j, kᴺ] - Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts) ℵi = ice_concentration[i, j, 1] - - ℐꜜˢʷ = downwelling_radiation.ℐꜜˢʷ[i, j, 1] - ℐꜜˡʷ = downwelling_radiation.ℐꜜˡʷ[i, j, 1] - 𝒬ᵀ = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible heat flux - 𝒬ᵛ = atmosphere_sea_ice_fluxes.latent_heat[i, j, 1] # latent heat flux + 𝒬ᵀ = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible heat flux + 𝒬ᵛ = atmosphere_sea_ice_fluxes.latent_heat[i, j, 1] # latent heat flux 𝒬ᶠʳᶻ = sea_ice_ocean_fluxes.frazil_heat[i, j, 1] # frazil heat flux 𝒬ⁱⁿᵗ = sea_ice_ocean_fluxes.interface_heat[i, j, 1] # interfacial heat flux end @@ -84,15 +60,8 @@ end ρτˣ = atmosphere_sea_ice_fluxes.x_momentum # zonal momentum flux ρτʸ = atmosphere_sea_ice_fluxes.y_momentum # meridional momentum flux - # Compute radiation fluxes - σ = atmos_sea_ice_properties.radiation.σ - α = atmos_sea_ice_properties.radiation.α - ϵ = atmos_sea_ice_properties.radiation.ϵ - ℐꜛˡʷ = emitted_longwave_radiation(i, j, kᴺ, grid, time, Ts, σ, ϵ) - ℐₜˢʷ = transmitted_shortwave_radiation(i, j, kᴺ, grid, time, α, ℐꜜˢʷ) - ℐₐˡʷ = absorbed_longwave_radiation(i, j, kᴺ, grid, time, ϵ, ℐꜜˡʷ) - - ΣQt = (ℐₜˢʷ + ℐₐˡʷ + ℐꜛˡʷ + 𝒬ᵀ + 𝒬ᵛ) * (ℵi > 0) # If ℵi == 0 there is no heat flux from the top! + # Turbulent contributions only (radiation added later by apply_air_sea_ice_radiative_fluxes!) + ΣQt = (𝒬ᵀ + 𝒬ᵛ) * (ℵi > 0) ΣQb = 𝒬ᶠʳᶻ + 𝒬ⁱⁿᵗ # Mask fluxes over land for convenience diff --git a/src/SeaIces/freezing_limited_ocean_temperature.jl b/src/SeaIces/freezing_limited_ocean_temperature.jl index d721b0b95..2860962f0 100644 --- a/src/SeaIces/freezing_limited_ocean_temperature.jl +++ b/src/SeaIces/freezing_limited_ocean_temperature.jl @@ -24,7 +24,7 @@ The melting temperature is a function of salinity and is controlled by the `liqu FreezingLimitedOceanTemperature(FT::DataType=Oceananigans.defaults.FloatType; liquidus=LinearLiquidus(FT)) = FreezingLimitedOceanTemperature(liquidus) -const FreezingLimitedEarthSystemModel = EarthSystemModel{<:FreezingLimitedOceanTemperature, A, L, O, <:NoSeaIceInterface} where {A, L, O} +const FreezingLimitedEarthSystemModel = EarthSystemModel{R, A, L, <:FreezingLimitedOceanTemperature, O, <:NoSeaIceInterface} where {R, A, L, O} # Extend interface methods to work with a `FreezingLimitedOceanTemperature` sea_ice_concentration(::FreezingLimitedOceanTemperature) = ZeroField() @@ -50,8 +50,8 @@ InterfaceComputations.sea_ice_ocean_interface(grid, ::FreezingLimitedOceanTemper InterfaceComputations.net_fluxes(::FreezingLimitedOceanTemperature) = nothing -const OnlyOceanwithFreezingLimited = EarthSystemModel{<:FreezingLimitedOceanTemperature, <:Nothing, <:Any, <:Any} -const OnlyAtmospherewithFreezingLimited = EarthSystemModel{<:FreezingLimitedOceanTemperature, <:Any, <:Any, <:Nothing} +const OnlyOceanwithFreezingLimited = EarthSystemModel{<:Any, <:Nothing, <:Any, <:FreezingLimitedOceanTemperature, <:Any} +const OnlyAtmospherewithFreezingLimited = EarthSystemModel{<:Any, <:Any, <:Any, <:FreezingLimitedOceanTemperature, <:Nothing} const SingleComponentPlusFreezingLimited = Union{OnlyAtmospherewithFreezingLimited, OnlyOceanwithFreezingLimited} # Also for the ocean nothing really happens here diff --git a/test/runtests.jl b/test/runtests.jl index ef9d1bc50..7a88b6ccb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -79,6 +79,9 @@ function __init__() try atmosphere = JRA55PrescribedAtmosphere(backend=JRA55NetCDFBackend(2)) + # Touch the radiation variables (rlds/rsds) too, so a corrupted cached + # download is caught by the same fallback path. + radiation = JRA55PrescribedRadiation(backend=JRA55NetCDFBackend(2)) catch e @warn "Original JRA55 download failed, trying NumericalEarthArtifacts fallback..." exception=(e, catch_backtrace()) emit_ci_warning("Broken JRA55 download", "Original source failed during init") @@ -87,6 +90,7 @@ function __init__() download_from_artifacts(metadata_path(datum)) end atmosphere = JRA55PrescribedAtmosphere(backend=JRA55NetCDFBackend(2)) + radiation = JRA55PrescribedRadiation(backend=JRA55NetCDFBackend(2)) end ##### diff --git a/test/test_checkpointer.jl b/test/test_checkpointer.jl index 3d058af42..ab13917b4 100644 --- a/test/test_checkpointer.jl +++ b/test/test_checkpointer.jl @@ -30,8 +30,9 @@ using Oceananigans.OutputWriters: Checkpointer backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) land = JRA55PrescribedLand(arch; backend) + radiation = JRA55PrescribedRadiation(arch; backend) - return OceanSeaIceModel(ocean, sea_ice; atmosphere, land) + return OceanSeaIceModel(sea_ice, ocean; atmosphere, land, radiation) end # Reference run: 3 iterations, then continue to 6 diff --git a/test/test_diagnostics_2.jl b/test/test_diagnostics_2.jl index 59e6d1c68..28fec3701 100644 --- a/test/test_diagnostics_2.jl +++ b/test/test_diagnostics_2.jl @@ -23,7 +23,7 @@ for arch in test_architectures sea_ice = sea_ice_simulation(grid, ocean) atmosphere = PrescribedAtmosphere(grid, [0.0]) - esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation = Radiation()) + esm = OceanSeaIceModel(sea_ice, ocean; atmosphere) T_flux = ocean.model.tracers.T.boundary_conditions.top.condition S_flux = ocean.model.tracers.S.boundary_conditions.top.condition diff --git a/test/test_ecco_atmosphere.jl b/test/test_ecco_atmosphere.jl index 70ac2b8ca..12d120930 100644 --- a/test/test_ecco_atmosphere.jl +++ b/test/test_ecco_atmosphere.jl @@ -2,8 +2,9 @@ include("runtests_setup.jl") include("download_utils.jl") using Statistics: median -using NumericalEarth.Atmospheres: PrescribedAtmosphere, TwoBandDownwellingRadiation -using NumericalEarth.ECCO: ECCOPrescribedAtmosphere, ECCO4Monthly +using NumericalEarth.Atmospheres: PrescribedAtmosphere +using NumericalEarth.Radiations: PrescribedRadiation +using NumericalEarth.ECCO: ECCOPrescribedAtmosphere, ECCOPrescribedRadiation, ECCO4Monthly using NumericalEarth.DataWrangling: download_dataset, metadata_path, higher_bound # Pre-download ECCO4Monthly atmospheric forcing variables through the artifacts @@ -34,7 +35,14 @@ end end_date, time_indices_in_memory = 2) + radiation = ECCOPrescribedRadiation(arch; + dataset, + start_date, + end_date, + time_indices_in_memory = 2) + @test atmosphere isa PrescribedAtmosphere + @test radiation isa PrescribedRadiation # Test that all expected fields are present @test haskey(atmosphere.velocities, :u) @@ -42,12 +50,11 @@ end @test haskey(atmosphere.tracers, :T) @test haskey(atmosphere.tracers, :q) @test !isnothing(atmosphere.pressure) - @test !isnothing(atmosphere.downwelling_radiation) @test haskey(atmosphere.freshwater_flux, :rain) # Test downwelling radiation components - ℐꜜˢʷ = atmosphere.downwelling_radiation.shortwave - ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave + ℐꜜˢʷ = radiation.downwelling_shortwave + ℐꜜˡʷ = radiation.downwelling_longwave @test ℐꜜˢʷ isa FieldTimeSeries @test ℐꜜˡʷ isa FieldTimeSeries diff --git a/test/test_ocean_only_model.jl b/test/test_ocean_only_model.jl index 6254f9e63..09c5e4c4f 100644 --- a/test/test_ocean_only_model.jl +++ b/test/test_ocean_only_model.jl @@ -19,7 +19,7 @@ using Oceananigans.OrthogonalSphericalShellGrids add_callback!(ocean, pushdata) backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) coupled_model = OceanOnlyModel(ocean; atmosphere, radiation) Δt = 60 for n = 1:3 @@ -50,7 +50,7 @@ using Oceananigans.OrthogonalSphericalShellGrids backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) # Fluxes are computed when the model is constructed, so we just test that this works. @test begin diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl index 4188435ec..abfc3200e 100644 --- a/test/test_ocean_sea_ice_model.jl +++ b/test/test_ocean_sea_ice_model.jl @@ -59,12 +59,12 @@ using ClimaSeaIce.Rheologies backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) # Fluxes are computed when the model is constructed, so we just test that this works. # And that we can time step with sea ice @test begin - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation) time_step!(coupled_model, 1) true end @@ -79,7 +79,7 @@ using ClimaSeaIce.Rheologies sea_ice_with_land = sea_ice_simulation(grid, ocean_with_land; advection=WENO(order=7)) above_freezing_ocean_temperature!(ocean_with_land, grid, sea_ice_with_land) - coupled_model = OceanSeaIceModel(ocean_with_land, sea_ice_with_land; atmosphere, land, radiation) + coupled_model = OceanSeaIceModel(sea_ice_with_land, ocean_with_land; atmosphere, land, radiation) @test !isnothing(coupled_model.interfaces.exchanger.land) time_step!(coupled_model, 1) true diff --git a/test/test_ospapa.jl b/test/test_ospapa.jl index 81b2cccf9..15a1af718 100644 --- a/test/test_ospapa.jl +++ b/test/test_ospapa.jl @@ -2,6 +2,7 @@ include("runtests_setup.jl") using NumericalEarth.OSPapa using NumericalEarth.Atmospheres: PrescribedAtmosphere +using NumericalEarth.Radiations: PrescribedRadiation using Oceananigans.BoundaryConditions: BoundaryCondition, Flux, getbc using Oceananigans.Units: minutes using CUDA: @allowscalar @@ -18,7 +19,12 @@ const OSPAPA_TEST_END = DateTime(2012, 10, 3) start_date = OSPAPA_TEST_START, end_date = OSPAPA_TEST_END) + radiation = OSPapaPrescribedRadiation(arch; + start_date = OSPAPA_TEST_START, + end_date = OSPAPA_TEST_END) + @test atmosphere isa PrescribedAtmosphere + @test radiation isa PrescribedRadiation # All expected fields are present @test haskey(atmosphere.velocities, :u) @@ -26,12 +32,11 @@ const OSPAPA_TEST_END = DateTime(2012, 10, 3) @test haskey(atmosphere.tracers, :T) @test haskey(atmosphere.tracers, :q) @test !isnothing(atmosphere.pressure) - @test !isnothing(atmosphere.downwelling_radiation) @test haskey(atmosphere.freshwater_flux, :rain) # Radiation sanity checks - ℐꜜˢʷ = atmosphere.downwelling_radiation.shortwave - ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave + ℐꜜˢʷ = radiation.downwelling_shortwave + ℐꜜˡʷ = radiation.downwelling_longwave @allowscalar begin sw_data = interior(ℐꜜˢʷ) @@ -216,8 +221,11 @@ end atmosphere = OSPapaPrescribedAtmosphere(arch; start_date = OSPAPA_TEST_START, end_date = OSPAPA_TEST_END) + radiation = OSPapaPrescribedRadiation(arch; + start_date = OSPAPA_TEST_START, + end_date = OSPAPA_TEST_END) - coupled_model = OceanOnlyModel(ocean; atmosphere, radiation=Radiation(arch)) + coupled_model = OceanOnlyModel(ocean; atmosphere, radiation) simulation = Simulation(coupled_model; Δt=ocean.Δt, stop_iteration=2) @test begin diff --git a/test/test_radiations.jl b/test/test_radiations.jl new file mode 100644 index 000000000..ac651db89 --- /dev/null +++ b/test/test_radiations.jl @@ -0,0 +1,82 @@ +include("runtests_setup.jl") + +using NumericalEarth.Radiations: PrescribedRadiation, + SurfaceRadiationProperties, + InterfaceRadiationFlux + +@testset "PrescribedRadiation construction" begin + for arch in test_architectures + A = typeof(arch) + + # Form B: grid-only constructor (zero downwelling, surface properties only) + @info "Testing PrescribedRadiation(grid) on $A..." + grid = RectilinearGrid(arch, size = 10, z = (-100, 0), topology = (Flat, Flat, Bounded)) + rad = PrescribedRadiation(grid) + @test rad isa PrescribedRadiation + @test rad.surface_properties isa NamedTuple + @test haskey(rad.surface_properties, :ocean) + @test haskey(rad.surface_properties, :sea_ice) + @test rad.surface_properties.ocean isa SurfaceRadiationProperties + @test rad.surface_properties.ocean.albedo == 0.05 + @test rad.surface_properties.ocean.emissivity == 0.97 + @test rad.surface_properties.sea_ice.albedo == 0.7 + @test rad.surface_properties.sea_ice.emissivity == 1.0 + @test rad.stefan_boltzmann_constant ≈ 5.67e-8 atol=1e-10 + @test isnothing(rad.interface_fluxes) + + # Surfaces can be omitted + @info "Testing PrescribedRadiation(grid; sea_ice_surface=nothing) on $A..." + rad_ocean_only = PrescribedRadiation(grid; sea_ice_surface = nothing) + @test haskey(rad_ocean_only.surface_properties, :ocean) + @test !haskey(rad_ocean_only.surface_properties, :sea_ice) + + # Custom surface properties + custom_ocean = SurfaceRadiationProperties(0.1, 0.95) + rad_custom = PrescribedRadiation(grid; ocean_surface = custom_ocean) + @test rad_custom.surface_properties.ocean.albedo == 0.1 + @test rad_custom.surface_properties.ocean.emissivity == 0.95 + + # time_step! works + @info "Testing time_step!(::PrescribedRadiation) on $A..." + rad2 = PrescribedRadiation(grid) + time_step!(rad2, 60.0) + @test rad2.clock.time == 60.0 + end +end + +@testset "PrescribedRadiation paired with model" begin + for arch in test_architectures + A = typeof(arch) + + @info "Testing OceanOnlyModel with PrescribedRadiation on $A..." + + grid = RectilinearGrid(arch, size = 10, z = (-100, 0), topology = (Flat, Flat, Bounded)) + ocean = ocean_simulation(grid) + radiation = PrescribedRadiation(grid) + model = OceanOnlyModel(ocean; radiation) + + # interface_fluxes are allocated for present surfaces (ocean + sea_ice + # via FreezingLimitedOceanTemperature). + @test !isnothing(model.radiation.interface_fluxes) + @test model.radiation.interface_fluxes.ocean isa InterfaceRadiationFlux + @test model.radiation.interface_fluxes.sea_ice isa InterfaceRadiationFlux + + time_step!(model, 60) + @test iteration(model) == 1 + end +end + +@testset "JRA55PrescribedRadiation" begin + for arch in test_architectures + A = typeof(arch) + @info "Testing JRA55PrescribedRadiation on $A..." + + backend = JRA55NetCDFBackend(2) + radiation = JRA55PrescribedRadiation(arch; backend) + + @test radiation isa PrescribedRadiation + @test radiation.downwelling_shortwave isa FieldTimeSeries + @test radiation.downwelling_longwave isa FieldTimeSeries + @test radiation.surface_properties.ocean isa SurfaceRadiationProperties + end +end diff --git a/test/test_reactant.jl b/test/test_reactant.jl index c37963121..fc32d8dff 100644 --- a/test/test_reactant.jl +++ b/test/test_reactant.jl @@ -36,8 +36,7 @@ end atmos_times = range(0, 360Oceananigans.Units.days, length=10) atmosphere = PrescribedAtmosphere(atmos_grid, atmos_times) - radiation = Radiation(arch) - coupled_model = OceanOnlyModel(ocean; atmosphere, radiation) + coupled_model = OceanOnlyModel(ocean; atmosphere) # Test that Reactant does _not_ initialize in the constructor for EarthSystemModel exchanger = coupled_model.interfaces.exchanger.atmosphere diff --git a/test/test_sea_ice_ocean_heat_fluxes.jl b/test/test_sea_ice_ocean_heat_fluxes.jl index 7506d3d78..b20601ab8 100644 --- a/test/test_sea_ice_ocean_heat_fluxes.jl +++ b/test/test_sea_ice_ocean_heat_fluxes.jl @@ -201,14 +201,14 @@ end backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) for sea_ice_ocean_heat_flux in [IceBathHeatFlux(), ThreeEquationHeatFlux()] @testset "Salt flux with $(nameof(typeof(sea_ice_ocean_heat_flux)))" begin interfaces = ComponentInterfaces(atmosphere, ocean, sea_ice; radiation, sea_ice_ocean_heat_flux) - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation, interfaces) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation, interfaces) # Test melting conditions: warm ocean above freezing # Freezing point at S=35 is about -1.9°C @@ -259,14 +259,14 @@ end backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) for sea_ice_ocean_heat_flux in [IceBathHeatFlux(), ThreeEquationHeatFlux()] @testset "Flux magnitude with $(nameof(typeof(sea_ice_ocean_heat_flux)))" begin interfaces = ComponentInterfaces(atmosphere, ocean, sea_ice; radiation, sea_ice_ocean_heat_flux) - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation, interfaces) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation, interfaces) # Set up melting conditions set!(ocean.model, T=2.0, S=35.0) # Warm ocean @@ -403,14 +403,14 @@ end backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) for sea_ice_ocean_heat_flux in [IceBathHeatFlux(), ThreeEquationHeatFlux()] @testset "Frazil with $(nameof(typeof(sea_ice_ocean_heat_flux)))" begin interfaces = ComponentInterfaces(atmosphere, ocean, sea_ice; radiation, sea_ice_ocean_heat_flux) - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation, interfaces) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation, interfaces) # Set up conditions where frazil might form: # Cold ocean near freezing with ice present @@ -454,11 +454,11 @@ end backend = JRA55NetCDFBackend(4) atmosphere = JRA55PrescribedAtmosphere(arch; backend) - radiation = Radiation(arch) + radiation = JRA55PrescribedRadiation(arch; backend) # Test with ThreeEquationHeatFlux (default) @test begin - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation) flux_form = coupled_model.interfaces.sea_ice_ocean_interface.flux_formulation flux_form isa ThreeEquationHeatFlux end @@ -469,7 +469,7 @@ end interfaces = ComponentInterfaces(atmosphere, ocean, sea_ice; radiation, sea_ice_ocean_heat_flux = flux) - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation, interfaces) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation, interfaces) flux_form = coupled_model.interfaces.sea_ice_ocean_interface.flux_formulation flux_form isa IceBathHeatFlux end @@ -482,7 +482,7 @@ end interfaces = ComponentInterfaces(atmosphere, ocean, sea_ice; radiation, sea_ice_ocean_heat_flux) - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation, interfaces) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation, interfaces) @test begin time_step!(coupled_model, 60) true diff --git a/test/test_speedy_coupling.jl b/test/test_speedy_coupling.jl index 59cf8c58b..413f030d5 100644 --- a/test/test_speedy_coupling.jl +++ b/test/test_speedy_coupling.jl @@ -15,8 +15,7 @@ Oceananigans.set!(ocean.model, T=EN4Metadatum(:temperature), S=EN4Metadatum(:sal atmos = NumericalEarth.atmosphere_simulation(spectral_grid) -radiation = Radiation(ocean_emissivity=0.0, sea_ice_emissivity=0.0) -earth_model = EarthSystemModel(atmos, ocean, default_sea_ice(); radiation) +earth_model = EarthSystemModel(; atmosphere=atmos, sea_ice=default_sea_ice(), ocean) Qca = atmos.variables.parameterizations.ocean.sensible_heat_flux.data Mva = atmos.variables.parameterizations.ocean.surface_humidity_flux.data diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index 3ffdd7b60..8f2cf51b6 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -75,15 +75,12 @@ end ρᵃᵗ = Thermodynamics.air_density(ℂᵃᵗ, Tᵃᵗ, pᵃᵗ, qᵃᵗ) ℰv = Thermodynamics.latent_heat_vapor(ℂᵃᵗ, Tᵃᵗ) - # No radiation equivalent - radiation = Radiation(ocean_emissivity=0, ocean_albedo=1) - - # turbulent fluxes that force a specific humidity at the ocean's surface + # No radiation: pass `radiation = nothing` to disable radiative + # contributions wholesale. for atmosphere_ocean_interface_temperature in (BulkTemperature(), SkinTemperature(DiffusiveFlux(1, 1e-2))) @info " Testing zero fluxes with $(atmosphere_ocean_interface_temperature)..." interfaces = ComponentInterfaces(atmosphere, ocean; - radiation, atmosphere_ocean_interface_specific_humidity, atmosphere_ocean_interface_temperature) @@ -246,7 +243,7 @@ end # Always cooling! fill!(atmosphere.tracers.T, 273.15 - 20) - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere) # Test that the temperature has snapped up to freezing @test minimum(ocean.model.tracers.T) == 0 @@ -286,7 +283,7 @@ end fill!(ocean.model.tracers.T, -2.0) # Test that we populate the sea-ice ocean stress - earth = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation=Radiation()) + earth = OceanSeaIceModel(sea_ice, ocean; atmosphere) τˣ = earth.interfaces.sea_ice_ocean_interface.fluxes.x_momentum τʸ = earth.interfaces.sea_ice_ocean_interface.fluxes.y_momentum @@ -330,7 +327,7 @@ end # radiation = Radiation(ocean_albedo=0.1, ocean_emissivity=1.0) # sea_ice = nothing -# coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +# coupled_model = OceanSeaIceModel(sea_ice, ocean; atmosphere, radiation) # times = 0:1hours:1days # Ntimes = length(times)