Refactor PrescribedAtmosphere and add an idealized atmosphere simulat…

…ion (#251) * Refactor PrescribedAtmosphere and add an idealized atmosphere simulation * Add arctic night case * Convert to JRA55PrescribedAtmosphere * Import thermodynamics parameters * bugfixes * Fix bug with Default coriolis * Fix syntax * Update src/OceanSeaIceModels/PrescribedAtmospheres.jl Co-authored-by: Simone Silvestri <[email protected]> --------- Co-authored-by: Navid C. Constantinou <[email protected]> Co-authored-by: Simone Silvestri <[email protected]>
CliMA · Dec 11, 2024 · c8dc984 · c8dc984
1 parent cd71fed
commit c8dc984
Show file tree

Hide file tree

Showing 24 changed files with 401 additions and 113 deletions.
diff --git a/README.md b/README.md
@@ -64,7 +64,7 @@ set!(ocean.model, T = ClimaOcean.ECCOMetadata(:temperature; date),
                   S = ClimaOcean.ECCOMetadata(:salinity; date))
 
 # Build and run an OceanSeaIceModel (with no sea ice component) forced by JRA55 reanalysis
-atmosphere = ClimaOcean.JRA55_prescribed_atmosphere(arch)
+atmosphere = ClimaOcean.JRA55PrescribedAtmosphere(arch)
 coupled_model = ClimaOcean.OceanSeaIceModel(ocean; atmosphere)
 simulation = Simulation(coupled_model, Δt=5minutes, stop_time=30days)
 run!(simulation)

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -54,7 +54,7 @@ set!(ocean.model, T = ClimaOcean.ECCOMetadata(:temperature; date),
                   S = ClimaOcean.ECCOMetadata(:salinity; date))
 
 # Build and run an OceanSeaIceModel (with no sea ice component) forced by JRA55 reanalysis
-atmosphere = ClimaOcean.JRA55_prescribed_atmosphere(arch)
+atmosphere = ClimaOcean.JRA55PrescribedAtmosphere(arch)
 coupled_model = ClimaOcean.OceanSeaIceModel(ocean; atmosphere)
 simulation = Simulation(coupled_model, Δt=5minutes, stop_time=30days)
 run!(simulation)

diff --git a/examples/arctic_night.jl b/examples/arctic_night.jl
@@ -0,0 +1,118 @@
+using ClimaOcean
+using Oceananigans
+using Oceananigans.Units
+using SeawaterPolynomials
+
+# Ocean state parameters
+T₀ = 0   # Surface temperature, ᵒC
+S₀ = 35   # Surface salinity
+N² = 1e-5 # Buoyancy gradient due to temperature stratification
+f = 0     # Coriolis parameter
+
+# Atmospheric state parameters
+Tₐ = 273.15 - 5 # Kelvin
+u₁₀ = 10 # wind at 10 m, m/s
+qₐ = 0.0 # specific humidity
+Qℓ = 0 # shortwave radiation (W m⁻², positive means heating right now)
+
+# Build the atmosphere
+atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+atmosphere_times = range(0, 1day, length=3)
+atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times)
+
+parent(atmosphere.tracers.T) .= Tₐ     # K
+parent(atmosphere.velocities.u) .= u₁₀ # m/s
+parent(atmosphere.tracers.q) .= qₐ     # mass ratio
+parent(atmosphere.downwelling_radiation.shortwave) .= Qℓ # W
+
+# Build ocean model at rest with initial temperature stratification
+grid = RectilinearGrid(size=400, z=(-400, 0), topology=(Flat, Flat, Bounded))
+ocean = ocean_simulation(grid, coriolis=FPlane(; f))
+
+eos = ocean.model.buoyancy.model.equation_of_state
+g = ocean.model.buoyancy.model.gravitational_acceleration
+α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos)
+dTdz = N² / (α * g)
+Tᵢ(z) = max(-1.8, T₀ + dTdz * z)
+set!(ocean.model, T=Tᵢ, S=S₀)
+
+radiation = Radiation(ocean_albedo=0.1)
+model = OceanSeaIceModel(ocean; atmosphere, radiation)
+simulation = Simulation(model, Δt=10.0, stop_time=3days)
+
+Q = model.fluxes.total.ocean.heat
+Ql = model.fluxes.turbulent.fields.latent_heat
+Qs = model.fluxes.turbulent.fields.sensible_heat
+τx = model.fluxes.turbulent.fields.x_momentum
+τy = model.fluxes.turbulent.fields.y_momentum
+Fv = model.fluxes.turbulent.fields.water_vapor
+
+fluxes = (; Q, Ql, Qs, τx, τy, Fv)
+outputs = merge(ocean.model.velocities, ocean.model.tracers, fluxes)
+ow = JLD2OutputWriter(ocean.model, outputs,
+                      schedule = TimeInterval(10minutes),
+                      filename = "idealized_atmosphere.jld2",
+                      overwrite_existing = true)
+
+simulation.output_writers[:ow] = ow
+
+run!(simulation)
+
+ut = FieldTimeSeries("idealized_atmosphere.jld2", "u")
+vt = FieldTimeSeries("idealized_atmosphere.jld2", "v")
+Tt = FieldTimeSeries("idealized_atmosphere.jld2", "T")
+St = FieldTimeSeries("idealized_atmosphere.jld2", "S")
+τx = FieldTimeSeries("idealized_atmosphere.jld2", "τx")
+Q = FieldTimeSeries("idealized_atmosphere.jld2", "Q")
+Nt = length(St)
+
+using GLMakie
+
+fig = Figure(size=(1000, 800))
+
+axτ = Axis(fig[1, 1:3], xlabel="Time (days)", ylabel="Zonal momentum \n flux (N m⁻²)", xaxisposition=:top)
+axQ = Axis(fig[2, 1:3], xlabel="Time (days)", ylabel="Heat flux \n (W m⁻²)")
+axT = Axis(fig[3, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)")
+axS = Axis(fig[3, 2], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)")
+axu = Axis(fig[3, 3], xlabel="Velocities (m s⁻¹)", ylabel="z (m)", yaxisposition=:right)
+
+slider = Slider(fig[4, 1:3], startvalue=1, range=1:Nt)
+n = slider.value
+
+un = @lift ut[$n]
+vn = @lift vt[$n]
+Tn = @lift Tt[$n]
+Sn = @lift St[$n]
+
+lines!(axT, Tn)
+lines!(axS, Sn)
+lines!(axu, un, label="u")
+lines!(axu, vn, label="v")
+axislegend(axu)
+
+t = τx.times ./ days
+tn = @lift t[$n]
+
+lines!(axτ, t, τx[:], label="τx")
+vlines!(axτ, tn)
+
+lines!(axQ, t, Q[:])
+vlines!(axQ, tn)
+
+rowsize!(fig.layout, 1, Relative(0.2))
+rowsize!(fig.layout, 2, Relative(0.2))
+
+xlims!(axS, 34.9, 35.1)
+xlims!(axu, -0.1, 4.0)
+hideydecorations!(axS, grid=false)
+hidespines!(axT, :t, :r)
+hidespines!(axS, :t, :l, :r)
+hidespines!(axu, :t, :l)
+hidespines!(axτ, :b, :r)
+hidespines!(axQ, :t, :r)
+
+record(fig, "idealized_atmosphere.mp4", 1:Nt, framerate=24) do nn
+    @info "Drawing frame $nn of $Nt..."
+    n[] = nn
+end
+
diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl
@@ -44,7 +44,7 @@ save("ECCO_continents.png", fig) #hide
 # We invoke the constructor with only the first two time indices, corresponding to 
 # January 1st (at 00:00 AM and 03:00 AM).
 
-atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory())
+atmosphere = JRA55PrescribedAtmosphere(1:2; backend = InMemory())
 ocean = ocean_simulation(grid)
 
 # Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity

diff --git a/examples/idealized_single_column_simulation.jl b/examples/idealized_single_column_simulation.jl
@@ -0,0 +1,118 @@
+using ClimaOcean
+using Oceananigans
+using Oceananigans.Units
+using SeawaterPolynomials
+
+# Ocean state parameters
+T₀ = 20   # Surface temperature, ᵒC
+S₀ = 35   # Surface salinity
+N² = 1e-4 # Buoyancy gradient due to temperature stratification
+f = 0     # Coriolis parameter
+
+# Atmospheric state parameters
+Tₐ = 273.15 + 20 # Kelvin
+u₁₀ = 10 # wind at 10 m, m/s
+qₐ = 0.01 # specific humidity
+Qℓ = 400 # shortwave radiation (W m⁻², positive means heating right now)
+
+# Build the atmosphere
+atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+atmosphere_times = range(0, 1day, length=3)
+atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times)
+
+parent(atmosphere.tracers.T) .= Tₐ     # K
+parent(atmosphere.velocities.u) .= u₁₀ # m/s
+parent(atmosphere.tracers.q) .= qₐ     # mass ratio
+parent(atmosphere.downwelling_radiation.shortwave) .= Qℓ # W
+
+# Build ocean model at rest with initial temperature stratification
+grid = RectilinearGrid(size=100, z=(-400, 0), topology=(Flat, Flat, Bounded))
+ocean = ocean_simulation(grid, coriolis=FPlane(; f))
+
+eos = ocean.model.buoyancy.model.equation_of_state
+g = ocean.model.buoyancy.model.gravitational_acceleration
+α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos)
+dTdz = N² / (α * g)
+Tᵢ(z) = T₀ + dTdz * z
+set!(ocean.model, T=Tᵢ, S=S₀)
+
+radiation = Radiation(ocean_albedo=0.1)
+model = OceanSeaIceModel(ocean; atmosphere, radiation)
+simulation = Simulation(model, Δt=5minutes, stop_time=30days)
+
+Q = model.fluxes.total.ocean.heat
+Ql = model.fluxes.turbulent.fields.latent_heat
+Qs = model.fluxes.turbulent.fields.sensible_heat
+τx = model.fluxes.turbulent.fields.x_momentum
+τy = model.fluxes.turbulent.fields.y_momentum
+Fv = model.fluxes.turbulent.fields.water_vapor
+
+fluxes = (; Q, Ql, Qs, τx, τy, Fv)
+outputs = merge(ocean.model.velocities, ocean.model.tracers, fluxes)
+ow = JLD2OutputWriter(ocean.model, outputs,
+                      schedule = TimeInterval(1hour),
+                      filename = "idealized_atmosphere.jld2",
+                      overwrite_existing = true)
+
+simulation.output_writers[:ow] = ow
+
+run!(simulation)
+
+ut = FieldTimeSeries("idealized_atmosphere.jld2", "u")
+vt = FieldTimeSeries("idealized_atmosphere.jld2", "v")
+Tt = FieldTimeSeries("idealized_atmosphere.jld2", "T")
+St = FieldTimeSeries("idealized_atmosphere.jld2", "S")
+τx = FieldTimeSeries("idealized_atmosphere.jld2", "τx")
+Q = FieldTimeSeries("idealized_atmosphere.jld2", "Q")
+Nt = length(St)
+
+using GLMakie
+
+fig = Figure(size=(1000, 800))
+
+axτ = Axis(fig[1, 1:3], xlabel="Time (days)", ylabel="Zonal momentum \n flux (N m⁻²)", xaxisposition=:top)
+axQ = Axis(fig[2, 1:3], xlabel="Time (days)", ylabel="Heat flux \n (W m⁻²)")
+axT = Axis(fig[3, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)")
+axS = Axis(fig[3, 2], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)")
+axu = Axis(fig[3, 3], xlabel="Velocities (m s⁻¹)", ylabel="z (m)", yaxisposition=:right)
+
+slider = Slider(fig[4, 1:3], startvalue=1, range=1:Nt)
+n = slider.value
+
+un = @lift ut[$n]
+vn = @lift vt[$n]
+Tn = @lift Tt[$n]
+Sn = @lift St[$n]
+
+lines!(axT, Tn)
+lines!(axS, Sn)
+lines!(axu, un, label="u")
+lines!(axu, vn, label="v")
+axislegend(axu)
+
+t = τx.times ./ days
+tn = @lift t[$n]
+
+lines!(axτ, t, τx[:], label="τx")
+vlines!(axτ, tn)
+
+lines!(axQ, t, Q[:])
+vlines!(axQ, tn)
+
+rowsize!(fig.layout, 1, Relative(0.2))
+rowsize!(fig.layout, 2, Relative(0.2))
+
+xlims!(axS, 34.9, 35.1)
+xlims!(axu, -0.1, 4.0)
+hideydecorations!(axS, grid=false)
+hidespines!(axT, :t, :r)
+hidespines!(axS, :t, :l, :r)
+hidespines!(axu, :t, :l)
+hidespines!(axτ, :b, :r)
+hidespines!(axQ, :t, :r)
+
+record(fig, "idealized_atmosphere.mp4", 1:Nt, framerate=24) do nn
+    @info "Drawing frame $nn of $Nt..."
+    n[] = nn
+end
+
diff --git a/examples/inspect_JRA55_data.jl b/examples/inspect_JRA55_data.jl
@@ -43,6 +43,7 @@ label = @lift string("Repeat-year JRA55 forcing on year-day ",
 Label(fig[0, 1:2], label, fontsize=24)
 
 sf = surface!(axsw, x, y, z, color=Qswn, colorrange=(0, 1200))
+
 Colorbar(fig[2, 1], sf,
          vertical = false,
          width = Relative(0.5),

diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl
@@ -107,7 +107,7 @@ radiation = Radiation(arch)
 # The number of snapshots that are loaded into memory is determined by
 # the `backend`. Here, we load 41 snapshots at a time into memory.
 
-atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(41))
+atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
 
 # ## The coupled simulation
 

diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl
@@ -59,10 +59,10 @@ set!(ocean.model, T=ECCOMetadata(:temperature), S=ECCOMetadata(:salinity))
 simulation_days = 31
 snapshots_per_day = 8 # corresponding to JRA55's 3-hour frequency
 last_time = simulation_days * snapshots_per_day
-atmosphere = JRA55_prescribed_atmosphere(1:last_time;
-                                         longitude = λ★,
-                                         latitude = φ★,
-                                         backend = InMemory())
+atmosphere = JRA55PrescribedAtmosphere(1:last_time;
+                                       longitude = λ★,
+                                       latitude = φ★,
+                                       backend = InMemory())
 
 # This builds a representation of the atmosphere on the small grid
 

diff --git a/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl b/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl
@@ -116,7 +116,7 @@ ocean.output_writers[:checkpoint] = Checkpointer(ocean.model,
 #####
 
 backend = JRA55NetCDFBackend(4) 
-atmosphere = JRA55_prescribed_atmosphere(arch; backend)
+atmosphere = JRA55PrescribedAtmosphere(arch; backend)
 radiation = Radiation(arch, ocean_albedo=LatitudeDependentAlbedo())
 coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
 

diff --git a/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl b/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl
@@ -15,7 +15,7 @@ using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: Radiation, SimilarityTheory
 using ClimaOcean.VerticalGrids: exponential_z_faces
 using ClimaOcean.JRA55
 using ClimaOcean.ECCO
-using ClimaOcean.JRA55: JRA55NetCDFBackend, JRA55_prescribed_atmosphere
+using ClimaOcean.JRA55: JRA55NetCDFBackend, JRA55PrescribedAtmosphere
 using ClimaOcean.ECCO: ECCO_restoring_forcing, ECCO4Monthly, ECCO2Daily, ECCOMetadata
 using ClimaOcean.Bathymetry
 using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: LatitudeDependentAlbedo
@@ -119,7 +119,7 @@ initial_date = dates[1]
 #####
 
 backend = JRA55NetCDFBackend(4) 
-atmosphere = JRA55_prescribed_atmosphere(arch; backend)
+atmosphere = JRA55PrescribedAtmosphere(arch; backend)
 radiation = Radiation(arch)
 coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
 

diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl
@@ -88,7 +88,7 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
 #####
 
 radiation  = Radiation(arch)
-atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20))
+atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20))
 
 #####
 ##### Coupled simulation

diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl
@@ -15,7 +15,8 @@ export
     SimilarityTheoryTurbulentFluxes,
     SkinTemperature,
     BulkTemperature,
-    JRA55_prescribed_atmosphere,
+    PrescribedAtmosphere,
+    JRA55PrescribedAtmosphere,
     JRA55NetCDFBackend,
     regrid_bathymetry,
     retrieve_bathymetry,
@@ -87,7 +88,9 @@ using .InitialConditions
 using .OceanSeaIceModels
 using .OceanSimulations
 using .DataWrangling: JRA55, ECCO
-using ClimaOcean.DataWrangling.JRA55: JRA55_prescribed_atmosphere, JRA55NetCDFBackend
+
+using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere
+using ClimaOcean.DataWrangling.JRA55: JRA55PrescribedAtmosphere, JRA55NetCDFBackend
 using ClimaOcean.DataWrangling.ECCO
 
 end # module