From 3114f74310c6de1985b7140a675ee0c79f54360e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 14:15:23 -0800 Subject: [PATCH 01/48] this works best? --- .../one_degree_simulation/one_degree_simulation.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 08577fed..17b2dc32 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -46,9 +46,8 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) catke = ClimaOcean.OceanSimulations.default_ocean_closure() -viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000) -closure = (gm, catke, viscous_closure) +closure = (gm, catke) ##### ##### Restoring @@ -72,7 +71,9 @@ forcing = (T=FT, S=FS) ##### Ocean simulation ##### -momentum_advection = VectorInvariant() +momentum_advection = WENOVectorInvariant(vorticity_order=5, + upwinding=OnlySelfUpwinding(cross_scheme=Centered())) + tracer_advection = Centered(order=2) # Should we add a side drag since this is at a coarser resolution? @@ -95,7 +96,7 @@ atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) ##### coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days) +simulation = Simulation(coupled_model; Δt=2minutes, stop_time=30days) ##### ##### Run it! @@ -127,3 +128,7 @@ end add_callback!(simulation, progress, IterationInterval(10)) run!(simulation) + +simulation.Δt = 25minutes + +run!(simulation) From 9101d5791c94764cb04d477da3a20fa3ba48c6fe Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 14:18:02 -0800 Subject: [PATCH 02/48] One degree simulation take 4 --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 70c54fca..660f3321 100644 --- a/Project.toml +++ b/Project.toml @@ -42,9 +42,9 @@ JLD2 = "0.4, 0.5" KernelAbstractions = "0.9" MPI = "0.20" NCDatasets = "0.12, 0.13, 0.14" -Oceananigans = "0.94.2 - 0.99" +Oceananigans = "0.94.3 - 0.99" # This needs to be 0.94.4 when PolarBoundaryCondition is implemented OffsetArrays = "1.14" -OrthogonalSphericalShellGrids = "0.1.8" +OrthogonalSphericalShellGrids = "0.1.9" Scratch = "1" SeawaterPolynomials = "0.3.4" StaticArrays = "1" From 5833d1945fdc0a12003dd8c5cd607b6c362ef1e2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 15:24:30 -0800 Subject: [PATCH 03/48] add a biharmonic viscosity --- experiments/one_degree_simulation/one_degree_simulation.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 17b2dc32..c551a051 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -46,8 +46,9 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) catke = ClimaOcean.OceanSimulations.default_ocean_closure() +viscous_closure = HorizontalScalarBiharmonicDiffusivity(ν = 1e11) -closure = (gm, catke) +closure = (gm, catke, viscous_closure) ##### ##### Restoring @@ -71,9 +72,7 @@ forcing = (T=FT, S=FS) ##### Ocean simulation ##### -momentum_advection = WENOVectorInvariant(vorticity_order=5, - upwinding=OnlySelfUpwinding(cross_scheme=Centered())) - +momentum_advection = VectorInvariant() tracer_advection = Centered(order=2) # Should we add a side drag since this is at a coarser resolution? From 0aa74dcd863f85ff2a1b2d76158a29a7a9e0558b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sat, 23 Nov 2024 06:27:25 -0500 Subject: [PATCH 04/48] infer maximum delta t --- src/OceanSimulations/OceanSimulations.jl | 29 +++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 8c95965b..f87ffe4b 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -42,10 +42,33 @@ default_or_override(override, alternative_default=nothing) = override # Some defaults default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) -# 70 substeps is a safe rule of thumb for an ocean at 1/4 - 1/10th of a degree -# TODO: pass the cfl and a given Δt to calculate the number of substeps? +function infer_maximum_Δt(grid) + Δx = mean(xspacings(grid)) + Δy = mean(yspacings(grid)) + Δθ = rad2deg(mean([Δx, Δy])) / grid.radius + + # The maximum Δt is roughly 40minutes / Δθ, giving: + # - 40 minutes for a 1 degree ocean + # - 20 minutes for a 1/4 degree ocean + # - 10 minutes for a 1/8 degree ocean + # - 5 minutes for a 1/16 degree ocean + # - 2.5 minutes for a 1/32 degree ocean + + return 40minutes / Δθ +end + +# Some defaults +default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) + const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} -default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70) + +function default_free_surface(grid::TripolarOfSomeKind; + fixed_Δt = infer_maximum_Δt(grid), + cfl = 0.8) + free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) + @info "Using a $(free_surface)" + return free_surface +end function default_ocean_closure() mixing_length = CATKEMixingLength(Cᵇ=0.01) From 7c87f2b22ba43e473e13f1d228621c7eef862245 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 10:16:17 +0100 Subject: [PATCH 05/48] changes --- examples/generate_bathymetry.jl | 2 +- src/OceanSimulations/OceanSimulations.jl | 7 ++----- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/examples/generate_bathymetry.jl b/examples/generate_bathymetry.jl index e063f379..5da3e435 100644 --- a/examples/generate_bathymetry.jl +++ b/examples/generate_bathymetry.jl @@ -38,7 +38,7 @@ grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1), # one interpolation passes and no restrictions on connected regions. # - `h_smooth` shows the output of the function with 40 interpolation passes, which results # in a smoother bathymetry. -# - `h_no_connected_regions` shows the output of the function with `connected_regions_allowed = 0`, which +# - `h_no_connected_regions` shows the output of the function with `major_basins = 1`, which # means that the function does not allow connected regions in the bathymetry (e.g., lakes) # and fills them with land. diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index f87ffe4b..7322b375 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -57,14 +57,11 @@ function infer_maximum_Δt(grid) return 40minutes / Δθ end -# Some defaults -default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) - const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} function default_free_surface(grid::TripolarOfSomeKind; fixed_Δt = infer_maximum_Δt(grid), - cfl = 0.8) + cfl = 0.7) free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) @info "Using a $(free_surface)" return free_surface @@ -98,7 +95,7 @@ default_tracer_advection() = FluxFormAdvection(WENO(order=7), # TODO: Specify the grid to a grid on the sphere; otherwise we can provide a different # function that requires latitude and longitude etc for computing coriolis=FPlane... function ocean_simulation(grid; - Δt = 5minutes, + Δt = infer_maximum_Δt(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), From 3e3b77d31054b9acf3d94787c2888f9ed535b8ae Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 11:20:30 +0100 Subject: [PATCH 06/48] this sohuld work --- .../one_degree_simulation.jl | 6 +- .../CrossRealmFluxes/CrossRealmFluxes.jl | 1 + .../atmosphere_ocean_fluxes.jl | 22 +++- .../similarity_theory_turbulent_fluxes.jl | 113 +++++++++++------- .../CrossRealmFluxes/surface_temperature.jl | 58 +++++++++ 5 files changed, 152 insertions(+), 48 deletions(-) create mode 100644 src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index c551a051..22fd8fec 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -10,7 +10,7 @@ using CUDA: @allowscalar, device! using Oceananigans.Grids: znode -arch = GPU() +arch = CPU() ##### ##### Grid and Bathymetry @@ -18,7 +18,7 @@ arch = GPU() Nx = 360 Ny = 180 -Nz = 100 +Nz = 50 z_faces = exponential_z_faces(; Nz, depth=5000, h=34) @@ -48,7 +48,7 @@ gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1 catke = ClimaOcean.OceanSimulations.default_ocean_closure() viscous_closure = HorizontalScalarBiharmonicDiffusivity(ν = 1e11) -closure = (gm, catke, viscous_closure) +closure = (catke, viscous_closure) ##### ##### Restoring diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl index fde3a09b..9cac3060 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl @@ -31,6 +31,7 @@ include("tabulated_albedo.jl") include("roughness_lengths.jl") include("stability_functions.jl") include("seawater_saturation_specific_humidity.jl") +include("surface_temperature.jl") include("similarity_theory_turbulent_fluxes.jl") include("ocean_sea_ice_surface_fluxes.jl") include("atmosphere_ocean_fluxes.jl") diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index 29981870..16da0e41 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -143,6 +143,7 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) ocean_state, coupled_model.fluxes.ocean_temperature_units, surface_atmosphere_state, + radiation_properties, atmosphere.reference_height, # height at which the state is known atmosphere.boundary_layer_height, atmosphere.thermodynamics_parameters) @@ -266,6 +267,7 @@ end ocean_state, ocean_temperature_units, surface_atmos_state, + radiation, atmosphere_reference_height, atmosphere_boundary_layer_height, atmos_thermodynamics_parameters) @@ -279,6 +281,8 @@ end Tₐ = surface_atmos_state.T[i, j, 1] pₐ = surface_atmos_state.p[i, j, 1] qₐ = surface_atmos_state.q[i, j, 1] + Rs = surface_atmos_state.Qs[i, j, 1] + Rℓ = surface_atmos_state.Qℓ[i, j, 1] # Extract state variables at cell centers # Ocean state @@ -322,11 +326,17 @@ end inactive = inactive_node(i, j, kᴺ, grid, c, c, c) maxiter = ifelse(inactive, 1, similarity_theory.maxiter) - turbulent_fluxes = compute_similarity_theory_fluxes(similarity_theory, - dynamic_ocean_state, - dynamic_atmos_state, - atmosphere_boundary_layer_height, - ℂₐ, g, ϰ, maxiter) + prescribed_heat_fluxes = net_downwelling_radiation(i, j, grid, time, radiation, Rs, Rℓ) + stefan_boltzmann_constant = radiation.stefan_boltzmann_constant + ocean_emissivity = stateindex(radiation.emission.ocean, i, j, 1, grid, time) + + turbulent_fluxes, surface_temperature = compute_similarity_theory_fluxes(similarity_theory, + dynamic_ocean_state, + dynamic_atmos_state, + prescribed_heat_fluxes, + (; stefan_boltzmann_constant, ocean_emissivity), + atmosphere_boundary_layer_height, + ℂₐ, g, ϰ, maxiter) # Store fluxes Qv = similarity_theory.fields.latent_heat @@ -334,6 +344,7 @@ end Fv = similarity_theory.fields.water_vapor ρτx = similarity_theory.fields.x_momentum ρτy = similarity_theory.fields.y_momentum + Ts = similarity_theory.fields.surface_temperature @inbounds begin # +0: cooling, -0: heating @@ -342,6 +353,7 @@ end Fv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.water_vapor) ρτx[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.x_momentum) ρτy[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.y_momentum) + Ts[i, j, 1] = surface_temperature end end diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index cdf1d8c1..4221e68c 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -40,6 +40,7 @@ struct SimilarityTheoryTurbulentFluxes{FT, UF, TP, S, W, R, B, V, F} water_mole_fraction :: W # mole fraction of H₂O in seawater roughness_lengths :: R # parameterization for turbulent fluxes similarity_profile_type :: B # similarity profile relating atmosphere to surface state + surface_temperature :: T # surface temperature either diagnostic or prescribed bulk_velocity :: V # bulk velocity scale for turbulent fluxes tolerance :: FT # solver option maxiter :: Int # solver option @@ -65,6 +66,7 @@ Adapt.adapt_structure(to, fluxes::STTF) = SimilarityTheoryTurbulentFluxes(adapt( adapt(to, fluxes.water_mole_fraction), adapt(to, fluxes.roughness_lengths), adapt(to, fluxes.similarity_profile_type), + adapt(to, fluxes.surface_temperature), adapt(to, fluxes.bulk_velocity), fluxes.tolerance, fluxes.maxiter, @@ -83,16 +85,17 @@ end function Base.show(io::IO, fluxes::SimilarityTheoryTurbulentFluxes) print(io, summary(fluxes), '\n', - "├── gravitational_acceleration: ", prettysummary(fluxes.gravitational_acceleration), '\n', - "├── von_karman_constant: ", prettysummary(fluxes.von_karman_constant), '\n', - "├── turbulent_prandtl_number: ", prettysummary(fluxes.turbulent_prandtl_number), '\n', - "├── gustiness_parameter: ", prettysummary(fluxes.gustiness_parameter), '\n', - "├── stability_functions: ", summary(fluxes.stability_functions), '\n', - "├── water_mole_fraction: ", summary(fluxes.water_mole_fraction), '\n', - "├── water_vapor_saturation: ", summary(fluxes.water_vapor_saturation), '\n', - "├── roughness_lengths: ", summary(fluxes.roughness_lengths), '\n', - "├── similarity_profile_type: ", summary(fluxes.similarity_profile_type), '\n', - "└── thermodynamics_parameters: ", summary(fluxes.thermodynamics_parameters)) + "├── gravitational_acceleration: ", prettysummary(fluxes.gravitational_acceleration), '\n', + "├── von_karman_constant: ", prettysummary(fluxes.von_karman_constant), '\n', + "├── turbulent_prandtl_number: ", prettysummary(fluxes.turbulent_prandtl_number), '\n', + "├── gustiness_parameter: ", prettysummary(fluxes.gustiness_parameter), '\n', + "├── stability_functions: ", summary(fluxes.stability_functions), '\n', + "├── water_mole_fraction: ", summary(fluxes.water_mole_fraction), '\n', + "├── water_vapor_saturation: ", summary(fluxes.water_vapor_saturation), '\n', + "├── roughness_lengths: ", summary(fluxes.roughness_lengths), '\n', + "├── similarity_profile_type: ", summary(fluxes.similarity_profile_type), '\n', + "├── surface_temperature: ", summary(fluxes.surface_temperature), '\n', + "└── thermodynamics_parameters: ", summary(fluxes.thermodynamics_parameters)) end const PATP = PrescribedAtmosphereThermodynamicsParameters @@ -115,6 +118,7 @@ struct RelativeVelocity end water_mole_fraction = convert(FT, 0.98), roughness_lengths = default_roughness_lengths(FT), similarity_profile_type = LogarithmicSimilarityProfile(), + surface_temperature = PrescribedSurfaceTemperature(), bulk_velocity = RelativeVelocity(), tolerance = 1e-8, maxiter = 100, @@ -159,6 +163,7 @@ function SimilarityTheoryTurbulentFluxes(FT::DataType = Float64; water_mole_fraction = convert(FT, 0.98), roughness_lengths = default_roughness_lengths(FT), similarity_profile_type = LogarithmicSimilarityProfile(), + surface_temperature = PrescribedSurfaceTemperature(), bulk_velocity = RelativeVelocity(), tolerance = 1e-8, maxiter = 100, @@ -174,20 +179,24 @@ function SimilarityTheoryTurbulentFluxes(FT::DataType = Float64; water_mole_fraction, roughness_lengths, similarity_profile_type, + surface_temperature, bulk_velocity, convert(FT, tolerance), maxiter, fields) end -function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; kw...) +function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature = PrescribedSurfaceTemperature(), kw...) water_vapor = Field{Center, Center, Nothing}(grid) latent_heat = Field{Center, Center, Nothing}(grid) sensible_heat = Field{Center, Center, Nothing}(grid) x_momentum = Field{Center, Center, Nothing}(grid) y_momentum = Field{Center, Center, Nothing}(grid) + T_surface = Field{Center, Center, Nothing}(grid) - fields = (; latent_heat, sensible_heat, water_vapor, x_momentum, y_momentum) + fields = (; latent_heat, sensible_heat, water_vapor, x_momentum, y_momentum, T_surface) + + surface_temperature = regularize_surface_temperature(surface_temperature, grid) return SimilarityTheoryTurbulentFluxes(eltype(grid); kw..., fields) end @@ -232,6 +241,8 @@ struct COARELogarithmicSimilarityProfile end surface_state, atmos_state, atmos_boundary_layer_height, + prescribed_heat_fluxes, # Possibly use in state_differences + radiation, thermodynamics_parameters, gravitational_acceleration, von_karman_constant, @@ -239,23 +250,19 @@ struct COARELogarithmicSimilarityProfile end # Prescribed difference between two states ℂₐ = thermodynamics_parameters - Δh, Δu, Δv, Δθ, Δq = state_differences(ℂₐ, - atmos_state, - surface_state, - gravitational_acceleration, - similarity_theory.bulk_velocity) - differences = (; u=Δu, v=Δv, θ=Δθ, q=Δq, h=Δh) # Initial guess for the characteristic scales u★, θ★, q★. # Does not really matter if we are sophisticated or not, it converges # in about 10 iterations no matter what... u★ = convert(eltype(Δh), 1e-4) Σ★ = SimilarityScales(u★, u★, u★) + Δu, Δv = velocity_differences(atmos_state, surface_state, bulk_velocity) - # The inital velocity scale assumes that - # the gustiness velocity `Uᴳ` is equal to 0.5 ms⁻¹. - # That will be refined later on. + # The inital velocity scale assumes that the gustiness velocity `Uᴳ` is equal to 0.5 ms⁻¹. + # The initial surface temperature is the same as the ocean temperature (in Δθ). + # These will be refined later on. + θs = AtmosphericThermodynamics.air_temperature(ℂₐ, surface_state.ts.temperature) FT = eltype(Δh) Uᴳᵢ² = convert(FT, 0.5^2) ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ²) @@ -266,16 +273,19 @@ struct COARELogarithmicSimilarityProfile end while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) Σ₀ = Σ★ - # Refine both the characteristic scale and the effective - # velocity difference ΔU, including gustiness. - Σ★, ΔU = refine_similarity_variables(Σ★, ΔU, - similarity_theory, - surface_state, - differences, - atmos_boundary_layer_height, - thermodynamics_parameters, - gravitational_acceleration, - von_karman_constant) + # Refine both the characteristic scale, the effective + # velocity difference ΔU, including gustiness, and the surface + # state temperature. + Σ★, θs, ΔU = refine_similarity_variables(Σ★, θs, ΔU, + similarity_theory, + atmos_state, + surface_state, + atmos_boundary_layer_height, + thermodynamics_parameters, + prescribed_heat_fluxes, + radiation, + gravitational_acceleration, + von_karman_constant) iteration += 1 end @@ -303,8 +313,8 @@ struct COARELogarithmicSimilarityProfile end x_momentum = + ρₐ * τx, y_momentum = + ρₐ * τy, ) - - return fluxes + + return fluxes, θs end # Iterating condition for the characteristic scales solvers @@ -360,19 +370,28 @@ end @inline velocity_differences(𝒰₁, 𝒰₀, ::RelativeVelocity) = @inbounds 𝒰₁.u[1] - 𝒰₀.u[1], 𝒰₁.u[2] - 𝒰₀.u[2] @inline velocity_differences(𝒰₁, 𝒰₀, ::WindVelocity) = @inbounds 𝒰₁.u[1], 𝒰₁.u[2] -@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, g, bulk_velocity) +@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, Σ★, g, surface_temperature, + prescribed_heat_fluxes, + radiation, + bulk_velocity) z₁ = 𝒰₁.z z₀ = 𝒰₀.z Δh = z₁ - z₀ Δu, Δv = velocity_differences(𝒰₁, 𝒰₀, bulk_velocity) - + # Thermodynamic state 𝒬₁ = 𝒰₁.ts 𝒬₀ = 𝒰₀.ts + ρₐ = AtmosphericThermodynamics.air_density(ℂₐ, 𝒬ₐ) + cₚ = AtmosphericThermodynamics.cp_m(ℂₐ, 𝒬ₐ) # moist heat capacity + ℰv = AtmosphericThermodynamics.latent_heat_vapor(ℂₐ, 𝒬ₐ) + + θ₀ = retrieve_temperature(surface_temperature, θ₀, ℂ, 𝒬₀, 𝒬₁, ρₐ, cₚ, ℰv, Σ★, + prescribed_heat_fluxes, + radiation) + θ₁ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₁) - θ₀ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) - cₚ = AtmosphericThermodynamics.cp_m(ℂ, 𝒬₁) # moist heat capacity # Temperature difference including the ``lapse rate'' `α = g / cₚ` Δθ = θ₁ - θ₀ + g / cₚ * Δh @@ -381,19 +400,33 @@ end q₀ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₀) Δq = q₁ - q₀ - return Δh, Δu, Δv, Δθ, Δq + return Δh, Δu, Δv, Δθ, Δq, θ₀ end @inline function refine_similarity_variables(estimated_characteristic_scales, + surface_temperature, velocity_scale, similarity_theory, + atmos_state, surface_state, - differences, atmos_boundary_layer_height, thermodynamics_parameters, + prescribed_heat_fluxes, + radiation, gravitational_acceleration, von_karman_constant) + Δh, Δu, Δv, Δθ, Δq, θ₀ = state_differences(thermodynamics_parameters, + atmos_state, + surface_state, + surface_temperature, + estimated_characteristic_scales, + gravitational_acceleration, + similarity_theory.surface_temperature, + prescribed_heat_fluxes, + radiation, + similarity_theory.bulk_velocity) + # "initial" scales because we will recompute them u★ = estimated_characteristic_scales.momentum θ★ = estimated_characteristic_scales.temperature @@ -452,5 +485,5 @@ end # New velocity difference accounting for gustiness ΔU = sqrt(Δu^2 + Δv^2 + Uᴳ^2) - return SimilarityScales(u★, θ★, q★), ΔU + return SimilarityScales(u★, θ★, q★), θ₀, ΔU end diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl new file mode 100644 index 00000000..43857976 --- /dev/null +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -0,0 +1,58 @@ +import Thermodynamics as AtmosphericThermodynamics + +struct PrescribedSurfaceTemperature end + +regularize_surface_temperature(::PrescribedSurfaceTemperature, grid) = PrescribedSurfaceTemperature() + +struct DiagnosticSurfaceTemperature{I} + internal_flux :: I +end + +struct DiffusiveFlux{Δz, K} + Δz :: Z + κ :: K # diffusivity in W m⁻² K⁻¹ +end + +function DiffusiveFlux(grid; κ = 0.2) + Δz = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) + return DiffusiveFlux(Δz, κ) +end + +regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = + DiagnosticSurfaceTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ)) + +@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo + Qn / F.κ * F.Δz + +# Do nothing (just copy the temperature) +@inline retrieve_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ + +# Change 𝒬₀ as a function of incoming and outgoing fluxes. The flaw here is that +# the ocean emissivity and albedo are fixed, but they might be a function of the +# surface temperature, so we might need to pass actually the radiation and the +# albedo and emissivity as arguments. +@inline function retrieve_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, + ρₐ, cₚ, ℰv, Σ★, + prescribed_heat_fluxes, + radiation_properties) + + Rd = prescribed_heat_fluxes # net downwelling radiation + ϵ = radiation_properties.ocean_emissivity + σ = radiation_properties.stefan_boltzmann_constant + + Ru = ϵ * σ * θ₀^4 # upwelling radiation (calculated explicitly) + Rn = Ru - Rd # net radiation + + u★ = Σ★.momentum + θ★ = Σ★.temperature + q★ = Σ★.water_vapor + + # sensible heat flux + latent heat flux + Qs = - ρₐ * u★ * (cₚ * θ★ + q★ * ℰv) + + # Net heat flux (positive out of the ocean) + Qn = Qs + Rn + θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) + + # surface temperature calculated as a balance of fluxes + return flux_balance_temperature(st.internal_flux, θo, Qn) +end From d7ac6fa4a0f5de3ee32cff165f02ec3029dce8ae Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 11:40:25 +0100 Subject: [PATCH 07/48] using CUDA --- src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 43857976..ce780cbf 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -1,3 +1,4 @@ +using CUDA: @allowscalar import Thermodynamics as AtmosphericThermodynamics struct PrescribedSurfaceTemperature end @@ -8,7 +9,7 @@ struct DiagnosticSurfaceTemperature{I} internal_flux :: I end -struct DiffusiveFlux{Δz, K} +struct DiffusiveFlux{Z, K} Δz :: Z κ :: K # diffusivity in W m⁻² K⁻¹ end From 9eaa668ab1baa5c64c1b8b9e22f53f34715fd5ca Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 11:57:02 +0100 Subject: [PATCH 08/48] tests NaN --- .../atmosphere_ocean_fluxes.jl | 7 +-- .../CrossRealmFluxes/radiation.jl | 9 +++ .../similarity_theory_turbulent_fluxes.jl | 61 ++++++++----------- .../CrossRealmFluxes/surface_temperature.jl | 22 ++++--- src/OceanSimulations/OceanSimulations.jl | 1 + 5 files changed, 55 insertions(+), 45 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index 16da0e41..43272d42 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -327,14 +327,13 @@ end maxiter = ifelse(inactive, 1, similarity_theory.maxiter) prescribed_heat_fluxes = net_downwelling_radiation(i, j, grid, time, radiation, Rs, Rℓ) - stefan_boltzmann_constant = radiation.stefan_boltzmann_constant - ocean_emissivity = stateindex(radiation.emission.ocean, i, j, 1, grid, time) + radiative_properties = local_radiation_properties(i, j, 1, grid, radiation, time) turbulent_fluxes, surface_temperature = compute_similarity_theory_fluxes(similarity_theory, dynamic_ocean_state, dynamic_atmos_state, prescribed_heat_fluxes, - (; stefan_boltzmann_constant, ocean_emissivity), + radiative_properties, atmosphere_boundary_layer_height, ℂₐ, g, ϰ, maxiter) @@ -344,7 +343,7 @@ end Fv = similarity_theory.fields.water_vapor ρτx = similarity_theory.fields.x_momentum ρτy = similarity_theory.fields.y_momentum - Ts = similarity_theory.fields.surface_temperature + Ts = similarity_theory.fields.T_surface @inbounds begin # +0: cooling, -0: heating diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl index 58b86595..52176fd8 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl @@ -91,3 +91,12 @@ function Base.show(io::IO, properties::SurfaceProperties) print(io, "└── sea_ice: ", summary(properties.sea_ice)) end +@inline local_radiation_properties(i, j, k, grid, ::Nothing, time) = nothing + +@inline function local_radiation_properties(i, j, k, grid, r::Radiation, time) + σ = r.stefan_boltzmann_constant + ϵ = stateindex(r.emission.ocean, i, j, k, grid, time) + α = stateindex(r.reflection.ocean, i, j, k, grid, time) + + return (; ϵ, α, σ) +end diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 4221e68c..968fb659 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -29,7 +29,7 @@ import SurfaceFluxes.Parameters: ##### Bulk turbulent fluxes based on similarity theory ##### -struct SimilarityTheoryTurbulentFluxes{FT, UF, TP, S, W, R, B, V, F} +struct SimilarityTheoryTurbulentFluxes{FT, UF, TP, S, W, R, B, T, V, F} gravitational_acceleration :: FT # parameter von_karman_constant :: FT # parameter turbulent_prandtl_number :: FT # parameter @@ -40,7 +40,7 @@ struct SimilarityTheoryTurbulentFluxes{FT, UF, TP, S, W, R, B, V, F} water_mole_fraction :: W # mole fraction of H₂O in seawater roughness_lengths :: R # parameterization for turbulent fluxes similarity_profile_type :: B # similarity profile relating atmosphere to surface state - surface_temperature :: T # surface temperature either diagnostic or prescribed + surface_temperature_type :: T # surface temperature either diagnostic or prescribed bulk_velocity :: V # bulk velocity scale for turbulent fluxes tolerance :: FT # solver option maxiter :: Int # solver option @@ -66,7 +66,7 @@ Adapt.adapt_structure(to, fluxes::STTF) = SimilarityTheoryTurbulentFluxes(adapt( adapt(to, fluxes.water_mole_fraction), adapt(to, fluxes.roughness_lengths), adapt(to, fluxes.similarity_profile_type), - adapt(to, fluxes.surface_temperature), + adapt(to, fluxes.surface_temperature_type), adapt(to, fluxes.bulk_velocity), fluxes.tolerance, fluxes.maxiter, @@ -94,7 +94,7 @@ function Base.show(io::IO, fluxes::SimilarityTheoryTurbulentFluxes) "├── water_vapor_saturation: ", summary(fluxes.water_vapor_saturation), '\n', "├── roughness_lengths: ", summary(fluxes.roughness_lengths), '\n', "├── similarity_profile_type: ", summary(fluxes.similarity_profile_type), '\n', - "├── surface_temperature: ", summary(fluxes.surface_temperature), '\n', + "├── surface_temperature: ", summary(fluxes.surface_temperature_type), '\n', "└── thermodynamics_parameters: ", summary(fluxes.thermodynamics_parameters)) end @@ -118,7 +118,7 @@ struct RelativeVelocity end water_mole_fraction = convert(FT, 0.98), roughness_lengths = default_roughness_lengths(FT), similarity_profile_type = LogarithmicSimilarityProfile(), - surface_temperature = PrescribedSurfaceTemperature(), + surface_temperature_type = PrescribedSurfaceTemperature(), bulk_velocity = RelativeVelocity(), tolerance = 1e-8, maxiter = 100, @@ -163,7 +163,7 @@ function SimilarityTheoryTurbulentFluxes(FT::DataType = Float64; water_mole_fraction = convert(FT, 0.98), roughness_lengths = default_roughness_lengths(FT), similarity_profile_type = LogarithmicSimilarityProfile(), - surface_temperature = PrescribedSurfaceTemperature(), + surface_temperature_type = PrescribedSurfaceTemperature(), bulk_velocity = RelativeVelocity(), tolerance = 1e-8, maxiter = 100, @@ -179,14 +179,14 @@ function SimilarityTheoryTurbulentFluxes(FT::DataType = Float64; water_mole_fraction, roughness_lengths, similarity_profile_type, - surface_temperature, + surface_temperature_type, bulk_velocity, convert(FT, tolerance), maxiter, fields) end -function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature = PrescribedSurfaceTemperature(), kw...) +function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature_type = PrescribedSurfaceTemperature(), kw...) water_vapor = Field{Center, Center, Nothing}(grid) latent_heat = Field{Center, Center, Nothing}(grid) sensible_heat = Field{Center, Center, Nothing}(grid) @@ -196,9 +196,9 @@ function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature fields = (; latent_heat, sensible_heat, water_vapor, x_momentum, y_momentum, T_surface) - surface_temperature = regularize_surface_temperature(surface_temperature, grid) + surface_temperature_type = regularize_surface_temperature(surface_temperature_type, grid) - return SimilarityTheoryTurbulentFluxes(eltype(grid); kw..., fields) + return SimilarityTheoryTurbulentFluxes(eltype(grid); surface_temperature_type, kw..., fields) end ##### @@ -242,7 +242,7 @@ struct COARELogarithmicSimilarityProfile end atmos_state, atmos_boundary_layer_height, prescribed_heat_fluxes, # Possibly use in state_differences - radiation, + radiative_properties, thermodynamics_parameters, gravitational_acceleration, von_karman_constant, @@ -250,22 +250,21 @@ struct COARELogarithmicSimilarityProfile end # Prescribed difference between two states ℂₐ = thermodynamics_parameters - + FT = eltype(ℂₐ) # Initial guess for the characteristic scales u★, θ★, q★. # Does not really matter if we are sophisticated or not, it converges # in about 10 iterations no matter what... - u★ = convert(eltype(Δh), 1e-4) + u★ = convert(FT, 1e-4) Σ★ = SimilarityScales(u★, u★, u★) - Δu, Δv = velocity_differences(atmos_state, surface_state, bulk_velocity) + Δu, Δv = velocity_differences(atmos_state, surface_state, similarity_theory.bulk_velocity) # The inital velocity scale assumes that the gustiness velocity `Uᴳ` is equal to 0.5 ms⁻¹. - # The initial surface temperature is the same as the ocean temperature (in Δθ). + # The initial surface temperature is the same as the ocean temperature. # These will be refined later on. - θs = AtmosphericThermodynamics.air_temperature(ℂₐ, surface_state.ts.temperature) - FT = eltype(Δh) + θs = AtmosphericThermodynamics.air_temperature(ℂₐ, surface_state.ts) Uᴳᵢ² = convert(FT, 0.5^2) - ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ²) + ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ²) # Initialize the solver iteration = 0 @@ -283,7 +282,7 @@ struct COARELogarithmicSimilarityProfile end atmos_boundary_layer_height, thermodynamics_parameters, prescribed_heat_fluxes, - radiation, + radiative_properties, gravitational_acceleration, von_karman_constant) iteration += 1 @@ -370,7 +369,7 @@ end @inline velocity_differences(𝒰₁, 𝒰₀, ::RelativeVelocity) = @inbounds 𝒰₁.u[1] - 𝒰₀.u[1], 𝒰₁.u[2] - 𝒰₀.u[2] @inline velocity_differences(𝒰₁, 𝒰₀, ::WindVelocity) = @inbounds 𝒰₁.u[1], 𝒰₁.u[2] -@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, Σ★, g, surface_temperature, +@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, Σ★, g, surface_temperature_type, prescribed_heat_fluxes, radiation, bulk_velocity) @@ -383,11 +382,11 @@ end 𝒬₁ = 𝒰₁.ts 𝒬₀ = 𝒰₀.ts - ρₐ = AtmosphericThermodynamics.air_density(ℂₐ, 𝒬ₐ) - cₚ = AtmosphericThermodynamics.cp_m(ℂₐ, 𝒬ₐ) # moist heat capacity - ℰv = AtmosphericThermodynamics.latent_heat_vapor(ℂₐ, 𝒬ₐ) + ρₐ = AtmosphericThermodynamics.air_density(ℂ, 𝒬₁) + cₚ = AtmosphericThermodynamics.cp_m(ℂ, 𝒬₁) # moist heat capacity + ℰv = AtmosphericThermodynamics.latent_heat_vapor(ℂ, 𝒬₁) - θ₀ = retrieve_temperature(surface_temperature, θ₀, ℂ, 𝒬₀, 𝒬₁, ρₐ, cₚ, ℰv, Σ★, + θ₀ = retrieve_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, 𝒬₁, ρₐ, cₚ, ℰv, Σ★, prescribed_heat_fluxes, radiation) @@ -422,7 +421,7 @@ end surface_temperature, estimated_characteristic_scales, gravitational_acceleration, - similarity_theory.surface_temperature, + similarity_theory.surface_temperature_type, prescribed_heat_fluxes, radiation, similarity_theory.bulk_velocity) @@ -444,7 +443,6 @@ end ℓq = similarity_theory.roughness_lengths.water_vapor β = similarity_theory.gustiness_parameter - h = differences.h ℂ = thermodynamics_parameters g = gravitational_acceleration 𝒬ₒ = surface_state.ts # thermodynamic state @@ -463,14 +461,9 @@ end # Transfer coefficients at height `h` profile_type = similarity_theory.similarity_profile_type - χu = ϰ / similarity_profile(profile_type, ψu, h, ℓu₀, L★) - χθ = ϰ / similarity_profile(profile_type, ψθ, h, ℓθ₀, L★) - χq = ϰ / similarity_profile(profile_type, ψq, h, ℓq₀, L★) - - Δu = differences.u - Δv = differences.v - Δθ = differences.θ - Δq = differences.q + χu = ϰ / similarity_profile(profile_type, ψu, Δh, ℓu₀, L★) + χθ = ϰ / similarity_profile(profile_type, ψθ, Δh, ℓθ₀, L★) + χq = ϰ / similarity_profile(profile_type, ψq, Δh, ℓq₀, L★) # u★ including gustiness u★ = χu * uτ diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index ce780cbf..392a0d2d 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -1,10 +1,21 @@ using CUDA: @allowscalar import Thermodynamics as AtmosphericThermodynamics +#### +#### Prescribed surface temperature +#### + struct PrescribedSurfaceTemperature end regularize_surface_temperature(::PrescribedSurfaceTemperature, grid) = PrescribedSurfaceTemperature() +# Do nothing (just copy the temperature) +@inline retrieve_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ + +#### +#### Diagnostic surface temperature calculated as a flux balance +#### + struct DiagnosticSurfaceTemperature{I} internal_flux :: I end @@ -24,9 +35,6 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, @inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo + Qn / F.κ * F.Δz -# Do nothing (just copy the temperature) -@inline retrieve_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ - # Change 𝒬₀ as a function of incoming and outgoing fluxes. The flaw here is that # the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass actually the radiation and the @@ -37,10 +45,7 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, radiation_properties) Rd = prescribed_heat_fluxes # net downwelling radiation - ϵ = radiation_properties.ocean_emissivity - σ = radiation_properties.stefan_boltzmann_constant - - Ru = ϵ * σ * θ₀^4 # upwelling radiation (calculated explicitly) + Ru = upwelling_radiation(θ₀, radiation_properties) # upwelling radiation (calculated explicitly) Rn = Ru - Rd # net radiation u★ = Σ★.momentum @@ -57,3 +62,6 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, # surface temperature calculated as a balance of fluxes return flux_balance_temperature(st.internal_flux, θo, Qn) end + +@inline upwelling_radiation(θ₀, ::Nothing) = zero(θ₀) +@inline upwelling_radiation(θ₀, r) = r.σ * r.ϵ * θ₀^4 \ No newline at end of file diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 7322b375..b9098120 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -16,6 +16,7 @@ using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEEquation using SeawaterPolynomials.TEOS10: TEOS10EquationOfState +using Statistics: mean using Oceananigans.BuoyancyModels: g_Earth using Oceananigans.Coriolis: Ω_Earth From 6d094b19394ca627ea8f00178901ee68d81cecff Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 13:40:54 +0100 Subject: [PATCH 09/48] test pass --- .../similarity_theory_turbulent_fluxes.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 968fb659..170a1b99 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -240,9 +240,9 @@ struct COARELogarithmicSimilarityProfile end @inline function compute_similarity_theory_fluxes(similarity_theory, surface_state, atmos_state, - atmos_boundary_layer_height, prescribed_heat_fluxes, # Possibly use in state_differences radiative_properties, + atmos_boundary_layer_height, thermodynamics_parameters, gravitational_acceleration, von_karman_constant, @@ -269,13 +269,13 @@ struct COARELogarithmicSimilarityProfile end # Initialize the solver iteration = 0 Σ₀ = Σ★ - + while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) Σ₀ = Σ★ # Refine both the characteristic scale, the effective # velocity difference ΔU, including gustiness, and the surface # state temperature. - Σ★, θs, ΔU = refine_similarity_variables(Σ★, θs, ΔU, + Σ★, θs, ΔU = refine_similarity_variables(Σ★, θs, ΔU, similarity_theory, atmos_state, surface_state, @@ -425,12 +425,12 @@ end prescribed_heat_fluxes, radiation, similarity_theory.bulk_velocity) - + # "initial" scales because we will recompute them u★ = estimated_characteristic_scales.momentum θ★ = estimated_characteristic_scales.temperature q★ = estimated_characteristic_scales.water_vapor - uτ = velocity_scale + ΔU = velocity_scale # Similarity functions from Edson et al. (2013) ψu = similarity_theory.stability_functions.momentum @@ -453,7 +453,7 @@ end # Monin-Obhukov characteristic length scale and non-dimensional height ϰ = von_karman_constant L★ = ifelse(b★ == 0, zero(b★), - u★^2 / (ϰ * b★)) - + # Compute roughness length scales ℓu₀ = roughness_length(ℓu, u★, 𝒬ₒ, ℂ) ℓq₀ = roughness_length(ℓq, ℓu₀, u★, 𝒬ₒ, ℂ) @@ -466,7 +466,7 @@ end χq = ϰ / similarity_profile(profile_type, ψq, Δh, ℓq₀, L★) # u★ including gustiness - u★ = χu * uτ + u★ = χu * ΔU θ★ = χθ * Δθ q★ = χq * Δq From 51328c8fc2ff0cdcc45f9811a7765c7d27f68993 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 13:43:09 +0100 Subject: [PATCH 10/48] regression test pass! --- .../CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 170a1b99..de5cc0f8 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -270,6 +270,9 @@ struct COARELogarithmicSimilarityProfile end iteration = 0 Σ₀ = Σ★ + # break the cycle if Δu == Δv == gustiness_parameter == 0 + + # Iterate until convergence while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) Σ₀ = Σ★ # Refine both the characteristic scale, the effective From 33f69abc721cdee1bde44911dfc05a111fe7a74b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 13:44:52 +0100 Subject: [PATCH 11/48] regression tests pass! --- .../CrossRealmFluxes/surface_temperature.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 392a0d2d..1c838282 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -25,6 +25,14 @@ struct DiffusiveFlux{Z, K} κ :: K # diffusivity in W m⁻² K⁻¹ end +# A default constructor for DiagnosticSurfaceTemperature +function DiagnosticSurfaceTemperature() + internal_flux = DiffusiveFlux() + return DiagnosticSurfaceTemperature(internal_flux) +end + +DiffusiveFlux(; κ = 0.2) = DiffusiveFlux(nothing, κ) + function DiffusiveFlux(grid; κ = 0.2) Δz = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) return DiffusiveFlux(Δz, κ) From d1a3caf754b9b602a6e9b6dda7a5dc1ea0f4d319 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 15:13:24 +0100 Subject: [PATCH 12/48] now it works --- .../os_papa_surface_temperature.jl | 301 ++++++++++++++++++ .../atmosphere_ocean_fluxes.jl | 6 + .../similarity_theory_turbulent_fluxes.jl | 33 +- .../CrossRealmFluxes/surface_temperature.jl | 19 +- 4 files changed, 339 insertions(+), 20 deletions(-) create mode 100644 experiments/single_column_experiments/os_papa_surface_temperature.jl diff --git a/experiments/single_column_experiments/os_papa_surface_temperature.jl b/experiments/single_column_experiments/os_papa_surface_temperature.jl new file mode 100644 index 00000000..3261d611 --- /dev/null +++ b/experiments/single_column_experiments/os_papa_surface_temperature.jl @@ -0,0 +1,301 @@ +# This simulation tests the difference between using a +# `PrescribedSurfaceTemperature` and a `DiagnosticSurfaceTemperature` + +using ClimaOcean +using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: SimilarityTheoryTurbulentFluxes, DiagnosticSurfaceTemperature +using Oceananigans +using Oceananigans.Units +using Oceananigans.BuoyancyModels: buoyancy_frequency +using Oceananigans.Units: Time +using Printf + +# Ocean station papa location +location_name = "ocean_station_papa" +λ★, φ★ = 35.1, 50.1 + +grid = RectilinearGrid(size = 200, + x = λ★, + y = φ★, + z = (-400, 0), + topology = (Flat, Flat, Bounded)) + +ocean1 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) +# ocean2 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) + +# We set initial conditions from ECCO: +# set!(ocean1.model.tracers.T, ECCOMetadata(:temperature)) +# set!(ocean1.model.tracers.S, ECCOMetadata(:salinity)) +set!(ocean2.model.tracers.T, ECCOMetadata(:temperature)) +set!(ocean2.model.tracers.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()) + +# We continue constructing a simulation. + +radiation = Radiation() +coupled_model_prescribed = OceanSeaIceModel(ocean1; atmosphere, radiation) + +similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = DiagnosticSurfaceTemperature()) +coupled_model_diagnostic = OceanSeaIceModel(ocean2; atmosphere, radiation, similarity_theory) + +for (coupled_model, suffix) in zip([coupled_model_prescribed, coupled_model_diagnostic], + ["diagnostic", "prescribed"]) + + simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) + + wall_clock = Ref(time_ns()) + + function progress(sim) + msg = "Ocean Station Papa" + msg *= string(", iter: ", iteration(sim), ", time: ", prettytime(sim)) + + elapsed = 1e-9 * (time_ns() - wall_clock[]) + msg *= string(", wall time: ", prettytime(elapsed)) + wall_clock[] = time_ns() + + u, v, w = sim.model.ocean.model.velocities + msg *= @sprintf(", max|u|: (%.2e, %.2e)", maximum(abs, u), maximum(abs, v)) + + T = sim.model.ocean.model.tracers.T + S = sim.model.ocean.model.tracers.S + e = sim.model.ocean.model.tracers.e + + τx = first(sim.model.fluxes.total.ocean.momentum.u) + τy = first(sim.model.fluxes.total.ocean.momentum.v) + Q = first(sim.model.fluxes.total.ocean.heat) + + u★ = sqrt(sqrt(τx^2 + τy^2)) + + Nz = size(T, 3) + msg *= @sprintf(", u★: %.2f m s⁻¹", u★) + msg *= @sprintf(", Q: %.2f W m⁻²", Q) + msg *= @sprintf(", T₀: %.2f ᵒC", first(interior(T, 1, 1, Nz))) + msg *= @sprintf(", extrema(T): (%.2f, %.2f) ᵒC", minimum(T), maximum(T)) + msg *= @sprintf(", S₀: %.2f g/kg", first(interior(S, 1, 1, Nz))) + msg *= @sprintf(", e₀: %.2e m² s⁻²", first(interior(e, 1, 1, Nz))) + + @info msg + end + + simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) + + # Build flux outputs + τx = coupled_model.fluxes.total.ocean.momentum.u + τy = coupled_model.fluxes.total.ocean.momentum.v + JT = coupled_model.fluxes.total.ocean.tracers.T + Js = coupled_model.fluxes.total.ocean.tracers.S + E = coupled_model.fluxes.turbulent.fields.water_vapor + Qc = coupled_model.fluxes.turbulent.fields.sensible_heat + Qv = coupled_model.fluxes.turbulent.fields.latent_heat + Ts = coupled_model.fluxes.turbulent.fields.T_surface + ρₒ = coupled_model.fluxes.ocean_reference_density + cₚ = coupled_model.fluxes.ocean_heat_capacity + + Q = ρₒ * cₚ * JT + ρτx = ρₒ * τx + ρτy = ρₒ * τy + N² = buoyancy_frequency(ocean.model) + κc = ocean.model.diffusivity_fields.κc + + fluxes = (; ρτx, ρτy, E, Js, Qv, Qc, Ts) + auxiliary_fields = (; N², κc) + fields = merge(ocean.model.velocities, ocean.model.tracers, auxiliary_fields) + + # Slice fields at the surface + outputs = merge(fields, fluxes) + + filename = "single_column_omip_$(location_name)_$(suffix)" + + simulation.output_writers[:jld2] = JLD2OutputWriter(ocean.model, outputs; filename, + schedule = TimeInterval(3hours), + overwrite_existing = true) + + run!(simulation) +end + + +##### +##### Visualization +##### + +filename_prescribed = "single_column_omip_$(location_name)_prescribed" +filename_diagnostic = "single_column_omip_$(location_name)_diagnostic" + +# Diagnosed +ud = FieldTimeSeries(filename_diagnostic, "u") +vd = FieldTimeSeries(filename_diagnostic, "v") +ed = FieldTimeSeries(filename_diagnostic, "e") +Nd² = FieldTimeSeries(filename_diagnostic, "N²") +κd = FieldTimeSeries(filename_diagnostic, "κc") + +Tsd = FieldTimeSeries(filename_diagnostic, "Ts") +Qvd = FieldTimeSeries(filename_diagnostic, "Qv") +Qcd = FieldTimeSeries(filename_diagnostic, "Qc") +Jsd = FieldTimeSeries(filename_diagnostic, "Js") +Evd = FieldTimeSeries(filename_diagnostic, "E") +ρτxd = FieldTimeSeries(filename_diagnostic, "ρτx") +ρτyd = FieldTimeSeries(filename_diagnostic, "ρτy") + +# Prescribed +up = FieldTimeSeries(filename_diagnostic, "u") +vp = FieldTimeSeries(filename_diagnostic, "v") +ep = FieldTimeSeries(filename_diagnostic, "e") +Np² = FieldTimeSeries(filename_diagnostic, "N²") +κp = FieldTimeSeries(filename_diagnostic, "κc") + +Tsp = FieldTimeSeries(filename_diagnostic, "Ts") +Qvp = FieldTimeSeries(filename_diagnostic, "Qv") +Qcp = FieldTimeSeries(filename_diagnostic, "Qc") +Jsp = FieldTimeSeries(filename_diagnostic, "Js") +Evp = FieldTimeSeries(filename_diagnostic, "E") +ρτxp = FieldTimeSeries(filename_diagnostic, "ρτx") +ρτyp = FieldTimeSeries(filename_diagnostic, "ρτy") + +Nz = size(T, 3) +times = Qc.times + +fig = Figure(size=(1800, 1800)) + +axτ = Axis(fig[1, 1:3], xlabel="Days since Oct 1 1992", ylabel="Wind stress (N m⁻²)") +axQ = Axis(fig[1, 4:6], xlabel="Days since Oct 1 1992", ylabel="Heat flux (W m⁻²)") +axu = Axis(fig[2, 1:3], xlabel="Days since Oct 1 1992", ylabel="Velocities (m s⁻¹)") +axT = Axis(fig[2, 4:6], xlabel="Days since Oct 1 1992", ylabel="Surface temperature (ᵒC)") +axF = Axis(fig[3, 1:3], xlabel="Days since Oct 1 1992", ylabel="Freshwater volume flux (m s⁻¹)") +axS = Axis(fig[3, 4:6], xlabel="Days since Oct 1 1992", ylabel="Surface salinity (g kg⁻¹)") + +axuz = Axis(fig[4:5, 1:2], xlabel="Velocities (m s⁻¹)", ylabel="z (m)") +axTz = Axis(fig[4:5, 3:4], xlabel="Temperature (ᵒC)", ylabel="z (m)") +axSz = Axis(fig[4:5, 5:6], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)") +axNz = Axis(fig[6:7, 1:2], xlabel="Buoyancy frequency (s⁻²)", ylabel="z (m)") +axκz = Axis(fig[6:7, 3:4], xlabel="Eddy diffusivity (m² s⁻¹)", ylabel="z (m)", xscale=log10) +axez = Axis(fig[6:7, 5:6], xlabel="Turbulent kinetic energy (m² s⁻²)", ylabel="z (m)", xscale=log10) + +title = @sprintf("Single-column simulation at %.2f, %.2f", φ★, λ★) +Label(fig[0, 1:6], title) + +n = Observable(1) + +times = (times .- times[1]) ./days +Nt = length(times) +tn = @lift times[$n] + +colors = Makie.wong_colors() + +ρₒ = coupled_model_prescribed.fluxes.ocean_reference_density +τxp = interior(ρτxp, 1, 1, 1, :) ./ ρₒ +τyp = interior(ρτyp, 1, 1, 1, :) ./ ρₒ +up★ = @. (τxp^2 + τyp^2)^(1/4) +τxd = interior(ρτxd, 1, 1, 1, :) ./ ρₒ +τyd = interior(ρτyd, 1, 1, 1, :) ./ ρₒ +ud★ = @. (τxd^2 + τyd^2)^(1/4) + +lines!(axu, times, interior(ud, 1, 1, Nz, :), color=colors[1], label="Zonal") +lines!(axu, times, interior(vd, 1, 1, Nz, :), color=colors[2], label="Meridional") +lines!(axu, times, interior(up, 1, 1, Nz, :), color=colors[1], linestyle = :dashed, label="Zonal") +lines!(axu, times, interior(vp, 1, 1, Nz, :), color=colors[2], linestyle = :dashed, label="Meridional") +lines!(axu, times, ud★, color=colors[3], label="Ocean-side u★") +lines!(axu, times, up★, color=colors[3], label="Ocean-side u★", linestyle = :dashed) +vlines!(axu, tn, linewidth=4, color=(:black, 0.5)) +axislegend(axu) + +lines!(axτ, times, interior(ρτxd, 1, 1, 1, :), label="Zonal") +lines!(axτ, times, interior(ρτyd, 1, 1, 1, :), label="Meridional") +lines!(axτ, times, interior(ρτxp, 1, 1, 1, :), linestyle=:dashed, label="Zonal") +lines!(axτ, times, interior(ρτyp, 1, 1, 1, :), linestyle=:dashed, label="Meridional") +vlines!(axτ, tn, linewidth=4, color=(:black, 0.5)) +axislegend(axτ) + +lines!(axT, times, interior(Tsd, 1, 1, 1, :), color=colors[2], linewidth=4, label="Ocean surface temperature") +lines!(axT, times, interior(Tsp, 1, 1, 1, :), color=colors[2], linewidth=4, linestyle = :dashed, label="Ocean surface temperature") +vlines!(axT, tn, linewidth=4, color=(:black, 0.5)) +axislegend(axT) + +lines!(axQ, times, interior(Qvd, 1, 1, 1, 1:Nt), color=colors[2], label="Sensible", linewidth=2) +lines!(axQ, times, interior(Qcd, 1, 1, 1, 1:Nt), color=colors[3], label="Latent", linewidth=2) +lines!(axQ, times, - interior(Qswd, 1, 1, 1, 1:Nt), color=colors[4], label="Shortwave", linewidth=2) +lines!(axQ, times, - interior(Qlwd, 1, 1, 1, 1:Nt), color=colors[5], label="Longwave", linewidth=2) +lines!(axQ, times, interior(Qvp, 1, 1, 1, 1:Nt), color=colors[2], linestyle=:dashed, label="Sensible", linewidth=2) +lines!(axQ, times, interior(Qcp, 1, 1, 1, 1:Nt), color=colors[3], linestyle=:dashed, label="Latent", linewidth=2) +lines!(axQ, times, - interior(Qswp, 1, 1, 1, 1:Nt), color=colors[4], linestyle=:dashed, label="Shortwave", linewidth=2) +lines!(axQ, times, - interior(Qlwp, 1, 1, 1, 1:Nt), color=colors[5], linestyle=:dashed, label="Longwave", linewidth=2) +vlines!(axQ, tn, linewidth=4, color=(:black, 0.5)) +axislegend(axQ) + +lines!(axF, times, Ptd[1:Nt], label="Prescribed freshwater flux") +lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), label="Evaporation") +lines!(axF, times, Ptd[1:Nt], linestyle=:dashed, label="Prescribed freshwater flux") +lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), linestyle=:dashed, label="Evaporation") +vlines!(axF, tn, linewidth=4, color=(:black, 0.5)) +axislegend(axF) + +lines!(axS, times, interior(S, 1, 1, Nz, :)) +vlines!(axS, tn, linewidth=4, color=(:black, 0.5)) + +zc = znodes(T) +zf = znodes(κ) +upn = @lift interior(up[$n], 1, 1, :) +vpn = @lift interior(vp[$n], 1, 1, :) +Tpn = @lift interior(Tp[$n], 1, 1, :) +Spn = @lift interior(Sp[$n], 1, 1, :) +κpn = @lift interior(κp[$n], 1, 1, :) +epn = @lift interior(ep[$n], 1, 1, :) +Np²n = @lift interior(Np²[$n], 1, 1, :) + +udn = @lift interior(ud[$n], 1, 1, :) +vdn = @lift interior(vd[$n], 1, 1, :) +Tdn = @lift interior(Td[$n], 1, 1, :) +Sdn = @lift interior(Sd[$n], 1, 1, :) +κdn = @lift interior(κd[$n], 1, 1, :) +edn = @lift interior(ed[$n], 1, 1, :) +Nd²n = @lift interior(Nd²[$n], 1, 1, :) + +scatterlines!(axuz, udn, zc, label="u") +scatterlines!(axuz, vdn, zc, label="v") +scatterlines!(axTz, Tdn, zc) +scatterlines!(axSz, Sdn, zc) +scatterlines!(axez, edn, zc) +scatterlines!(axNz, Nd²n, zf) +scatterlines!(axκz, κdn, zf) +scatterlines!(axuz, upn, zc, linestyle=:dashed, label="u") +scatterlines!(axuz, vpn, zc, linestyle=:dashed, label="v") +scatterlines!(axTz, Tpn, zc, linestyle=:dashed) +scatterlines!(axSz, Spn, zc, linestyle=:dashed) +scatterlines!(axez, epn, zc, linestyle=:dashed) +scatterlines!(axNz, Np²n, zf, linestyle=:dashed) +scatterlines!(axκz, κpn, zf, linestyle=:dashed) + +axislegend(axuz) + +ulim = max(maximum(abs, up), maximum(abs, vp)) +xlims!(axuz, -ulim, ulim) + +Tmax = maximum(interior(Tp)) +Tmin = minimum(interior(Tp)) +xlims!(axTz, Tmin - 0.1, Tmax + 0.1) + +Nmax = maximum(interior(Np²)) +xlims!(axNz, -Nmax/10, Nmax * 1.05) + +κmax = maximum(interior(κp)) +xlims!(axκz, 1e-9, κmax * 1.1) + +emax = maximum(interior(ep)) +xlims!(axez, 1e-11, emax * 1.1) + +Smax = maximum(interior(Sp)) +Smin = minimum(interior(Sp)) +xlims!(axSz, Smin - 0.2, Smax + 0.2) + +record(fig, "single_column_profiles.mp4", 1:Nt, framerate=24) do nn + @info "Drawing frame $nn of $Nt..." + n[] = nn +end +nothing #hide + +# ![](single_column_profiles.mp4) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index 43272d42..422bb9fb 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -141,6 +141,8 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) grid, clock, ocean_state, + coupled_model.fluxes.ocean_reference_density, + coupled_model.fluxes.ocean_heat_capacity, coupled_model.fluxes.ocean_temperature_units, surface_atmosphere_state, radiation_properties, @@ -265,6 +267,8 @@ end grid, clock, ocean_state, + ocean_density, + ocean_heat_capacity, ocean_temperature_units, surface_atmos_state, radiation, @@ -334,6 +338,8 @@ end dynamic_atmos_state, prescribed_heat_fluxes, radiative_properties, + ocean_density, + ocean_heat_capacity, atmosphere_boundary_layer_height, ℂₐ, g, ϰ, maxiter) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index de5cc0f8..de641ad7 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -24,7 +24,6 @@ import SurfaceFluxes.Parameters: universal_func_type, grav - ##### ##### Bulk turbulent fluxes based on similarity theory ##### @@ -242,6 +241,8 @@ struct COARELogarithmicSimilarityProfile end atmos_state, prescribed_heat_fluxes, # Possibly use in state_differences radiative_properties, + ocean_density, + ocean_heat_capacity, atmos_boundary_layer_height, thermodynamics_parameters, gravitational_acceleration, @@ -265,12 +266,14 @@ struct COARELogarithmicSimilarityProfile end θs = AtmosphericThermodynamics.air_temperature(ℂₐ, surface_state.ts) Uᴳᵢ² = convert(FT, 0.5^2) ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ²) + + # break the cycle if Δu == Δv == gustiness_parameter == 0 since that would mean + # that u★ == 0 so there is no turbulent transfer and the solver will not converge, leading to NaNs. + zero_shear_velocity = (Δu == 0) && (Δv == 0) && (similarity_theory.gustiness_parameter == 0) # Initialize the solver - iteration = 0 - Σ₀ = Σ★ - - # break the cycle if Δu == Δv == gustiness_parameter == 0 + iteration = ifelse(zero_shear_velocity, maxiter+1, 0) + Σ₀ = ifelse(zero_shear_velocity, SimilarityScales(0, 0, 0), Σ★) # Iterate until convergence while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) @@ -282,6 +285,8 @@ struct COARELogarithmicSimilarityProfile end similarity_theory, atmos_state, surface_state, + ocean_density, + ocean_heat_capacity, atmos_boundary_layer_height, thermodynamics_parameters, prescribed_heat_fluxes, @@ -372,9 +377,9 @@ end @inline velocity_differences(𝒰₁, 𝒰₀, ::RelativeVelocity) = @inbounds 𝒰₁.u[1] - 𝒰₀.u[1], 𝒰₁.u[2] - 𝒰₀.u[2] @inline velocity_differences(𝒰₁, 𝒰₀, ::WindVelocity) = @inbounds 𝒰₁.u[1], 𝒰₁.u[2] -@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, Σ★, g, surface_temperature_type, +@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, Σ★, g, ρₒ, cpₒ, surface_temperature_type, prescribed_heat_fluxes, - radiation, + radiative_properties, bulk_velocity) z₁ = 𝒰₁.z z₀ = 𝒰₀.z @@ -389,10 +394,10 @@ end cₚ = AtmosphericThermodynamics.cp_m(ℂ, 𝒬₁) # moist heat capacity ℰv = AtmosphericThermodynamics.latent_heat_vapor(ℂ, 𝒬₁) - θ₀ = retrieve_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, 𝒬₁, ρₐ, cₚ, ℰv, Σ★, + θ₀ = retrieve_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, prescribed_heat_fluxes, - radiation) - + radiative_properties) + θ₁ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₁) # Temperature difference including the ``lapse rate'' `α = g / cₚ` @@ -411,10 +416,12 @@ end similarity_theory, atmos_state, surface_state, + ocean_density, + ocean_heat_capacity, atmos_boundary_layer_height, thermodynamics_parameters, prescribed_heat_fluxes, - radiation, + radiative_properties, gravitational_acceleration, von_karman_constant) @@ -423,10 +430,12 @@ end surface_state, surface_temperature, estimated_characteristic_scales, + ocean_density, + ocean_heat_capacity, gravitational_acceleration, similarity_theory.surface_temperature_type, prescribed_heat_fluxes, - radiation, + radiative_properties, similarity_theory.bulk_velocity) # "initial" scales because we will recompute them diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 1c838282..f3b3fd08 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -22,7 +22,7 @@ end struct DiffusiveFlux{Z, K} Δz :: Z - κ :: K # diffusivity in W m⁻² K⁻¹ + κ :: K # diffusivity in m²s⁻¹ end # A default constructor for DiagnosticSurfaceTemperature @@ -31,9 +31,9 @@ function DiagnosticSurfaceTemperature() return DiagnosticSurfaceTemperature(internal_flux) end -DiffusiveFlux(; κ = 0.2) = DiffusiveFlux(nothing, κ) +DiffusiveFlux(; κ = 1e-2) = DiffusiveFlux(nothing, κ) -function DiffusiveFlux(grid; κ = 0.2) +function DiffusiveFlux(grid; κ = 1e-2) Δz = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) return DiffusiveFlux(Δz, κ) end @@ -48,13 +48,16 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, # surface temperature, so we might need to pass actually the radiation and the # albedo and emissivity as arguments. @inline function retrieve_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, - ρₐ, cₚ, ℰv, Σ★, + ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, prescribed_heat_fluxes, radiation_properties) - Rd = prescribed_heat_fluxes # net downwelling radiation - Ru = upwelling_radiation(θ₀, radiation_properties) # upwelling radiation (calculated explicitly) - Rn = Ru - Rd # net radiation + Rd = prescribed_heat_fluxes # net downwelling radiation (positive out of the ocean) + + # upwelling radiation is calculated explicitly + # TODO: we could calculate it semi-implicitly as ϵσTⁿ⁺¹Tⁿ³ + Ru = upwelling_radiation(θ₀, radiation_properties) + Rn = Ru + Rd # net radiation u★ = Σ★.momentum θ★ = Σ★.temperature @@ -64,7 +67,7 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, Qs = - ρₐ * u★ * (cₚ * θ★ + q★ * ℰv) # Net heat flux (positive out of the ocean) - Qn = Qs + Rn + Qn = (Qs + Rn) / ρₒ / cpₒ θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) # surface temperature calculated as a balance of fluxes From 1f545253dfe3a799b37efa90727b792741babad1 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 15:37:02 +0100 Subject: [PATCH 13/48] reverted to false? --- .../os_papa_surface_temperature.jl | 14 +++++++------- src/DataWrangling/ECCO/ECCO.jl | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/experiments/single_column_experiments/os_papa_surface_temperature.jl b/experiments/single_column_experiments/os_papa_surface_temperature.jl index 3261d611..095da987 100644 --- a/experiments/single_column_experiments/os_papa_surface_temperature.jl +++ b/experiments/single_column_experiments/os_papa_surface_temperature.jl @@ -20,11 +20,11 @@ grid = RectilinearGrid(size = 200, topology = (Flat, Flat, Bounded)) ocean1 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) -# ocean2 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) +ocean2 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) # We set initial conditions from ECCO: -# set!(ocean1.model.tracers.T, ECCOMetadata(:temperature)) -# set!(ocean1.model.tracers.S, ECCOMetadata(:salinity)) +set!(ocean1.model.tracers.T, ECCOMetadata(:temperature)) +set!(ocean1.model.tracers.S, ECCOMetadata(:salinity)) set!(ocean2.model.tracers.T, ECCOMetadata(:temperature)) set!(ocean2.model.tracers.S, ECCOMetadata(:salinity)) @@ -44,10 +44,10 @@ coupled_model_prescribed = OceanSeaIceModel(ocean1; atmosphere, radiation) similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = DiagnosticSurfaceTemperature()) coupled_model_diagnostic = OceanSeaIceModel(ocean2; atmosphere, radiation, similarity_theory) -for (coupled_model, suffix) in zip([coupled_model_prescribed, coupled_model_diagnostic], +for (coupled_model, suffix) in zip([coupled_model_diagnostic, coupled_model_prescribed], ["diagnostic", "prescribed"]) - simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) + simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=10days) wall_clock = Ref(time_ns()) @@ -71,12 +71,12 @@ for (coupled_model, suffix) in zip([coupled_model_prescribed, coupled_model_diag Q = first(sim.model.fluxes.total.ocean.heat) u★ = sqrt(sqrt(τx^2 + τy^2)) + Ts = first(interior(sim.model.fluxes.turbulent.fields.T_surface, 1, 1, 1)) Nz = size(T, 3) msg *= @sprintf(", u★: %.2f m s⁻¹", u★) msg *= @sprintf(", Q: %.2f W m⁻²", Q) - msg *= @sprintf(", T₀: %.2f ᵒC", first(interior(T, 1, 1, Nz))) - msg *= @sprintf(", extrema(T): (%.2f, %.2f) ᵒC", minimum(T), maximum(T)) + msg *= @sprintf(", T₀: %.2f ᵒC", Ts) msg *= @sprintf(", S₀: %.2f g/kg", first(interior(S, 1, 1, Nz))) msg *= @sprintf(", e₀: %.2e m² s⁻²", first(interior(e, 1, 1, Nz))) diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl index bdef42ce..a02bd07a 100644 --- a/src/DataWrangling/ECCO/ECCO.jl +++ b/src/DataWrangling/ECCO/ECCO.jl @@ -141,7 +141,7 @@ function ECCO_field(metadata::ECCOMetadata; inpainting = NearestNeighborInpainting(Inf), mask = nothing, horizontal_halo = (7, 7), - cache_inpainted_data = false) + cache_inpainted_data = true) field = empty_ECCO_field(metadata; architecture, horizontal_halo) inpainted_path = inpainted_metadata_path(metadata) From 3513848c58332ae02467c55c5ba0a7e05c498a22 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 15:53:21 +0100 Subject: [PATCH 14/48] the complete thing --- .../os_papa_surface_temperature.jl | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/experiments/single_column_experiments/os_papa_surface_temperature.jl b/experiments/single_column_experiments/os_papa_surface_temperature.jl index 095da987..69164be7 100644 --- a/experiments/single_column_experiments/os_papa_surface_temperature.jl +++ b/experiments/single_column_experiments/os_papa_surface_temperature.jl @@ -47,7 +47,7 @@ coupled_model_diagnostic = OceanSeaIceModel(ocean2; atmosphere, radiation, simil for (coupled_model, suffix) in zip([coupled_model_diagnostic, coupled_model_prescribed], ["diagnostic", "prescribed"]) - simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=10days) + simulation = Simulation(coupled_model, Δt=ocean1.Δt, stop_time=10days) wall_clock = Ref(time_ns()) @@ -100,37 +100,37 @@ for (coupled_model, suffix) in zip([coupled_model_diagnostic, coupled_model_pres Q = ρₒ * cₚ * JT ρτx = ρₒ * τx ρτy = ρₒ * τy - N² = buoyancy_frequency(ocean.model) - κc = ocean.model.diffusivity_fields.κc + N² = buoyancy_frequency(coupled_model.ocean.model) + κc = coupled_model.ocean.model.diffusivity_fields.κc fluxes = (; ρτx, ρτy, E, Js, Qv, Qc, Ts) auxiliary_fields = (; N², κc) - fields = merge(ocean.model.velocities, ocean.model.tracers, auxiliary_fields) + fields = merge(coupled_model.ocean.model.velocities, coupled_model.ocean.model.tracers, auxiliary_fields) # Slice fields at the surface outputs = merge(fields, fluxes) filename = "single_column_omip_$(location_name)_$(suffix)" - simulation.output_writers[:jld2] = JLD2OutputWriter(ocean.model, outputs; filename, + simulation.output_writers[:jld2] = JLD2OutputWriter(coupled_model.ocean.model, outputs; filename, schedule = TimeInterval(3hours), overwrite_existing = true) run!(simulation) end - ##### ##### Visualization ##### -filename_prescribed = "single_column_omip_$(location_name)_prescribed" -filename_diagnostic = "single_column_omip_$(location_name)_diagnostic" +filename_prescribed = "single_column_omip_$(location_name)_prescribed.jld2" +filename_diagnostic = "single_column_omip_$(location_name)_diagnostic.jld2" # Diagnosed ud = FieldTimeSeries(filename_diagnostic, "u") vd = FieldTimeSeries(filename_diagnostic, "v") -ed = FieldTimeSeries(filename_diagnostic, "e") +Td = FieldTimeSeries(filename_diagnostic, "v") +ed = FieldTimeSeries(filename_diagnostic, "T") Nd² = FieldTimeSeries(filename_diagnostic, "N²") κd = FieldTimeSeries(filename_diagnostic, "κc") @@ -145,6 +145,7 @@ Evd = FieldTimeSeries(filename_diagnostic, "E") # Prescribed up = FieldTimeSeries(filename_diagnostic, "u") vp = FieldTimeSeries(filename_diagnostic, "v") +Tp = FieldTimeSeries(filename_diagnostic, "T") ep = FieldTimeSeries(filename_diagnostic, "e") Np² = FieldTimeSeries(filename_diagnostic, "N²") κp = FieldTimeSeries(filename_diagnostic, "κc") @@ -157,8 +158,8 @@ Evp = FieldTimeSeries(filename_diagnostic, "E") ρτxp = FieldTimeSeries(filename_diagnostic, "ρτx") ρτyp = FieldTimeSeries(filename_diagnostic, "ρτy") -Nz = size(T, 3) -times = Qc.times +Nz = size(ud, 3) +times = Qcd.times fig = Figure(size=(1800, 1800)) @@ -292,10 +293,10 @@ Smax = maximum(interior(Sp)) Smin = minimum(interior(Sp)) xlims!(axSz, Smin - 0.2, Smax + 0.2) -record(fig, "single_column_profiles.mp4", 1:Nt, framerate=24) do nn - @info "Drawing frame $nn of $Nt..." - n[] = nn -end -nothing #hide +# record(fig, "single_column_profiles.mp4", 1:Nt, framerate=24) do nn +# @info "Drawing frame $nn of $Nt..." +# n[] = nn +# end +# nothing #hide # ![](single_column_profiles.mp4) From f754f7e3c54c8f09aabfe7d77084356018592972 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 16:13:57 +0100 Subject: [PATCH 15/48] chnage defaults --- .../os_papa_surface_temperature.jl | 92 +++++++++---------- .../CrossRealmFluxes/surface_temperature.jl | 4 +- 2 files changed, 46 insertions(+), 50 deletions(-) diff --git a/experiments/single_column_experiments/os_papa_surface_temperature.jl b/experiments/single_column_experiments/os_papa_surface_temperature.jl index 69164be7..cf428d82 100644 --- a/experiments/single_column_experiments/os_papa_surface_temperature.jl +++ b/experiments/single_column_experiments/os_papa_surface_temperature.jl @@ -41,7 +41,7 @@ atmosphere = JRA55_prescribed_atmosphere(1:last_time; radiation = Radiation() coupled_model_prescribed = OceanSeaIceModel(ocean1; atmosphere, radiation) -similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = DiagnosticSurfaceTemperature()) +similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = DiagnosticSurfaceTemperature(; κ = 0.5)) coupled_model_diagnostic = OceanSeaIceModel(ocean2; atmosphere, radiation, similarity_theory) for (coupled_model, suffix) in zip([coupled_model_diagnostic, coupled_model_prescribed], @@ -131,6 +131,7 @@ ud = FieldTimeSeries(filename_diagnostic, "u") vd = FieldTimeSeries(filename_diagnostic, "v") Td = FieldTimeSeries(filename_diagnostic, "v") ed = FieldTimeSeries(filename_diagnostic, "T") +Sd = FieldTimeSeries(filename_diagnostic, "S") Nd² = FieldTimeSeries(filename_diagnostic, "N²") κd = FieldTimeSeries(filename_diagnostic, "κc") @@ -143,20 +144,21 @@ Evd = FieldTimeSeries(filename_diagnostic, "E") ρτyd = FieldTimeSeries(filename_diagnostic, "ρτy") # Prescribed -up = FieldTimeSeries(filename_diagnostic, "u") -vp = FieldTimeSeries(filename_diagnostic, "v") -Tp = FieldTimeSeries(filename_diagnostic, "T") -ep = FieldTimeSeries(filename_diagnostic, "e") -Np² = FieldTimeSeries(filename_diagnostic, "N²") -κp = FieldTimeSeries(filename_diagnostic, "κc") - -Tsp = FieldTimeSeries(filename_diagnostic, "Ts") -Qvp = FieldTimeSeries(filename_diagnostic, "Qv") -Qcp = FieldTimeSeries(filename_diagnostic, "Qc") -Jsp = FieldTimeSeries(filename_diagnostic, "Js") -Evp = FieldTimeSeries(filename_diagnostic, "E") -ρτxp = FieldTimeSeries(filename_diagnostic, "ρτx") -ρτyp = FieldTimeSeries(filename_diagnostic, "ρτy") +up = FieldTimeSeries(filename_prescribed, "u") +vp = FieldTimeSeries(filename_prescribed, "v") +Tp = FieldTimeSeries(filename_prescribed, "T") +Sp = FieldTimeSeries(filename_prescribed, "S") +ep = FieldTimeSeries(filename_prescribed, "e") +Np² = FieldTimeSeries(filename_prescribed, "N²") +κp = FieldTimeSeries(filename_prescribed, "κc") + +Tsp = FieldTimeSeries(filename_prescribed, "Ts") +Qvp = FieldTimeSeries(filename_prescribed, "Qv") +Qcp = FieldTimeSeries(filename_prescribed, "Qc") +Jsp = FieldTimeSeries(filename_prescribed, "Js") +Evp = FieldTimeSeries(filename_prescribed, "E") +ρτxp = FieldTimeSeries(filename_prescribed, "ρτx") +ρτyp = FieldTimeSeries(filename_prescribed, "ρτy") Nz = size(ud, 3) times = Qcd.times @@ -198,48 +200,43 @@ ud★ = @. (τxd^2 + τyd^2)^(1/4) lines!(axu, times, interior(ud, 1, 1, Nz, :), color=colors[1], label="Zonal") lines!(axu, times, interior(vd, 1, 1, Nz, :), color=colors[2], label="Meridional") -lines!(axu, times, interior(up, 1, 1, Nz, :), color=colors[1], linestyle = :dashed, label="Zonal") -lines!(axu, times, interior(vp, 1, 1, Nz, :), color=colors[2], linestyle = :dashed, label="Meridional") +lines!(axu, times, interior(up, 1, 1, Nz, :), color=colors[1], linestyle = :dash, label="Zonal") +lines!(axu, times, interior(vp, 1, 1, Nz, :), color=colors[2], linestyle = :dash, label="Meridional") lines!(axu, times, ud★, color=colors[3], label="Ocean-side u★") -lines!(axu, times, up★, color=colors[3], label="Ocean-side u★", linestyle = :dashed) +lines!(axu, times, up★, color=colors[3], label="Ocean-side u★", linestyle = :dash) vlines!(axu, tn, linewidth=4, color=(:black, 0.5)) axislegend(axu) lines!(axτ, times, interior(ρτxd, 1, 1, 1, :), label="Zonal") lines!(axτ, times, interior(ρτyd, 1, 1, 1, :), label="Meridional") -lines!(axτ, times, interior(ρτxp, 1, 1, 1, :), linestyle=:dashed, label="Zonal") -lines!(axτ, times, interior(ρτyp, 1, 1, 1, :), linestyle=:dashed, label="Meridional") +lines!(axτ, times, interior(ρτxp, 1, 1, 1, :), linestyle=:dash, label="Zonal") +lines!(axτ, times, interior(ρτyp, 1, 1, 1, :), linestyle=:dash, label="Meridional") vlines!(axτ, tn, linewidth=4, color=(:black, 0.5)) axislegend(axτ) lines!(axT, times, interior(Tsd, 1, 1, 1, :), color=colors[2], linewidth=4, label="Ocean surface temperature") -lines!(axT, times, interior(Tsp, 1, 1, 1, :), color=colors[2], linewidth=4, linestyle = :dashed, label="Ocean surface temperature") +lines!(axT, times, interior(Tsp, 1, 1, 1, :), color=colors[2], linewidth=4, linestyle = :dash, label="Ocean surface temperature") vlines!(axT, tn, linewidth=4, color=(:black, 0.5)) axislegend(axT) -lines!(axQ, times, interior(Qvd, 1, 1, 1, 1:Nt), color=colors[2], label="Sensible", linewidth=2) -lines!(axQ, times, interior(Qcd, 1, 1, 1, 1:Nt), color=colors[3], label="Latent", linewidth=2) -lines!(axQ, times, - interior(Qswd, 1, 1, 1, 1:Nt), color=colors[4], label="Shortwave", linewidth=2) -lines!(axQ, times, - interior(Qlwd, 1, 1, 1, 1:Nt), color=colors[5], label="Longwave", linewidth=2) -lines!(axQ, times, interior(Qvp, 1, 1, 1, 1:Nt), color=colors[2], linestyle=:dashed, label="Sensible", linewidth=2) -lines!(axQ, times, interior(Qcp, 1, 1, 1, 1:Nt), color=colors[3], linestyle=:dashed, label="Latent", linewidth=2) -lines!(axQ, times, - interior(Qswp, 1, 1, 1, 1:Nt), color=colors[4], linestyle=:dashed, label="Shortwave", linewidth=2) -lines!(axQ, times, - interior(Qlwp, 1, 1, 1, 1:Nt), color=colors[5], linestyle=:dashed, label="Longwave", linewidth=2) +lines!(axQ, times, interior(Qvd, 1, 1, 1, 1:Nt), color=colors[2], label="Sensible", linewidth=2) +lines!(axQ, times, interior(Qcd, 1, 1, 1, 1:Nt), color=colors[3], label="Latent", linewidth=2) +lines!(axQ, times, interior(Qvp, 1, 1, 1, 1:Nt), color=colors[2], linestyle=:dash, label="Sensible", linewidth=2) +lines!(axQ, times, interior(Qcp, 1, 1, 1, 1:Nt), color=colors[3], linestyle=:dash, label="Latent", linewidth=2) vlines!(axQ, tn, linewidth=4, color=(:black, 0.5)) axislegend(axQ) -lines!(axF, times, Ptd[1:Nt], label="Prescribed freshwater flux") lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), label="Evaporation") -lines!(axF, times, Ptd[1:Nt], linestyle=:dashed, label="Prescribed freshwater flux") -lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), linestyle=:dashed, label="Evaporation") +lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), linestyle=:dash, label="Evaporation") vlines!(axF, tn, linewidth=4, color=(:black, 0.5)) axislegend(axF) -lines!(axS, times, interior(S, 1, 1, Nz, :)) +lines!(axS, times, interior(Sd, 1, 1, Nz, :)) +lines!(axS, times, interior(Sp, 1, 1, Nz, :), linestyle=:dash) vlines!(axS, tn, linewidth=4, color=(:black, 0.5)) -zc = znodes(T) -zf = znodes(κ) +zc = znodes(Tp) +zf = znodes(κp) upn = @lift interior(up[$n], 1, 1, :) vpn = @lift interior(vp[$n], 1, 1, :) Tpn = @lift interior(Tp[$n], 1, 1, :) @@ -263,13 +260,13 @@ scatterlines!(axSz, Sdn, zc) scatterlines!(axez, edn, zc) scatterlines!(axNz, Nd²n, zf) scatterlines!(axκz, κdn, zf) -scatterlines!(axuz, upn, zc, linestyle=:dashed, label="u") -scatterlines!(axuz, vpn, zc, linestyle=:dashed, label="v") -scatterlines!(axTz, Tpn, zc, linestyle=:dashed) -scatterlines!(axSz, Spn, zc, linestyle=:dashed) -scatterlines!(axez, epn, zc, linestyle=:dashed) -scatterlines!(axNz, Np²n, zf, linestyle=:dashed) -scatterlines!(axκz, κpn, zf, linestyle=:dashed) +scatterlines!(axuz, upn, zc, linestyle=:dash, label="u") +scatterlines!(axuz, vpn, zc, linestyle=:dash, label="v") +scatterlines!(axTz, Tpn, zc, linestyle=:dash) +scatterlines!(axSz, Spn, zc, linestyle=:dash) +scatterlines!(axez, epn, zc, linestyle=:dash) +scatterlines!(axNz, Np²n, zf, linestyle=:dash) +scatterlines!(axκz, κpn, zf, linestyle=:dash) axislegend(axuz) @@ -293,10 +290,9 @@ Smax = maximum(interior(Sp)) Smin = minimum(interior(Sp)) xlims!(axSz, Smin - 0.2, Smax + 0.2) -# record(fig, "single_column_profiles.mp4", 1:Nt, framerate=24) do nn -# @info "Drawing frame $nn of $Nt..." -# n[] = nn -# end -# nothing #hide - +record(fig, "single_column_profiles.mp4", 1:Nt, framerate=24) do nn + @info "Drawing frame $nn of $Nt..." + n[] = nn +end +nothing #hide # ![](single_column_profiles.mp4) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index f3b3fd08..d33d5615 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -26,8 +26,8 @@ struct DiffusiveFlux{Z, K} end # A default constructor for DiagnosticSurfaceTemperature -function DiagnosticSurfaceTemperature() - internal_flux = DiffusiveFlux() +function DiagnosticSurfaceTemperature(; κ = 0.1) + internal_flux = DiffusiveFlux(; κ) return DiagnosticSurfaceTemperature(internal_flux) end From e6d8829eb6ba5a53463f6c737f182c7655940807 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 17:16:38 +0100 Subject: [PATCH 16/48] switched a couple of variables --- .../similarity_theory_turbulent_fluxes.jl | 2 +- .../CrossRealmFluxes/surface_temperature.jl | 24 +++++++++++++++---- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index de641ad7..2f7130ca 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -430,9 +430,9 @@ end surface_state, surface_temperature, estimated_characteristic_scales, + gravitational_acceleration, ocean_density, ocean_heat_capacity, - gravitational_acceleration, similarity_theory.surface_temperature_type, prescribed_heat_fluxes, radiative_properties, diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index d33d5615..c9850383 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -41,7 +41,15 @@ end regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = DiagnosticSurfaceTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ)) -@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo + Qn / F.κ * F.Δz +# The balance solved is +# +# Θo - Θₛ +# κ ------- = Qn (all fluxes positive upwards) +# Δz +# +# Where the LHS is the internal diffusive flux inside the ocean and the RHS are the +# atmospheric and radiative fluxes (provided explicitly and iterated upon). +@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo - Qn / F.κ * F.Δz # Change 𝒬₀ as a function of incoming and outgoing fluxes. The flaw here is that # the ocean emissivity and albedo are fixed, but they might be a function of the @@ -57,20 +65,28 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, # upwelling radiation is calculated explicitly # TODO: we could calculate it semi-implicitly as ϵσTⁿ⁺¹Tⁿ³ Ru = upwelling_radiation(θ₀, radiation_properties) - Rn = Ru + Rd # net radiation + Rn = Rd + Ru u★ = Σ★.momentum θ★ = Σ★.temperature q★ = Σ★.water_vapor - - # sensible heat flux + latent heat flux + + # sensible heat flux + latent heat flux (positive out of the ocean) Qs = - ρₐ * u★ * (cₚ * θ★ + q★ * ℰv) # Net heat flux (positive out of the ocean) Qn = (Qs + Rn) / ρₒ / cpₒ + θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) # surface temperature calculated as a balance of fluxes + # t0 = + + # if θ₀ < 0 + # @show Qn, Qs, Rn, θ₀, u★, θ★, q★, ℰv + # throw(ArgumentError("Negative temperature")) + # end + return flux_balance_temperature(st.internal_flux, θo, Qn) end From 4e8faead9563a1eca7fa4471c469aa4d1a033fdb Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 17:29:02 +0100 Subject: [PATCH 17/48] remove comment --- .../CrossRealmFluxes/surface_temperature.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index c9850383..5402fb38 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -79,14 +79,6 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) - # surface temperature calculated as a balance of fluxes - # t0 = - - # if θ₀ < 0 - # @show Qn, Qs, Rn, θ₀, u★, θ★, q★, ℰv - # throw(ArgumentError("Negative temperature")) - # end - return flux_balance_temperature(st.internal_flux, θo, Qn) end From af28e4e8dbef3fe19b4cb42af43be428672363e2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 18:50:15 +0100 Subject: [PATCH 18/48] done better --- .../CrossRealmFluxes/surface_temperature.jl | 29 +++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 5402fb38..2bcfa2a9 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -2,12 +2,20 @@ using CUDA: @allowscalar import Thermodynamics as AtmosphericThermodynamics #### -#### Prescribed surface temperature +#### Utilities #### -struct PrescribedSurfaceTemperature end +@inline upwelling_radiation(θ₀, ::Nothing) = zero(θ₀) +@inline upwelling_radiation(θ₀, r) = r.σ * r.ϵ * θ₀^4 + +# For any surface temperture type that does not depend on the grid +regularize_surface_temperature(surface_temperature_type, grid) = surface_temperature_type + +#### +#### Prescribed surface temperature (the easiest case) +#### -regularize_surface_temperature(::PrescribedSurfaceTemperature, grid) = PrescribedSurfaceTemperature() +struct PrescribedSurfaceTemperature end # Do nothing (just copy the temperature) @inline retrieve_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ @@ -21,8 +29,8 @@ struct DiagnosticSurfaceTemperature{I} end struct DiffusiveFlux{Z, K} - Δz :: Z - κ :: K # diffusivity in m²s⁻¹ + δ :: Z # Boundary layer thickness, as a first guess we will use half the grid spacing + κ :: K # diffusivity in m²s⁻¹ end # A default constructor for DiagnosticSurfaceTemperature @@ -45,18 +53,18 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, # # Θo - Θₛ # κ ------- = Qn (all fluxes positive upwards) -# Δz +# δ # # Where the LHS is the internal diffusive flux inside the ocean and the RHS are the # atmospheric and radiative fluxes (provided explicitly and iterated upon). -@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo - Qn / F.κ * F.Δz +@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo - Qn / F.κ * F.δ # Change 𝒬₀ as a function of incoming and outgoing fluxes. The flaw here is that # the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass actually the radiation and the # albedo and emissivity as arguments. @inline function retrieve_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, - ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, + ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, prescribed_heat_fluxes, radiation_properties) @@ -80,7 +88,4 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) return flux_balance_temperature(st.internal_flux, θo, Qn) -end - -@inline upwelling_radiation(θ₀, ::Nothing) = zero(θ₀) -@inline upwelling_radiation(θ₀, r) = r.σ * r.ϵ * θ₀^4 \ No newline at end of file +end \ No newline at end of file From 9b32122ff3b06b6d1411694f92d3c84b08762427 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 18:50:58 +0100 Subject: [PATCH 19/48] better --- .../CrossRealmFluxes/surface_temperature.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 2bcfa2a9..c136aeba 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -41,9 +41,9 @@ end DiffusiveFlux(; κ = 1e-2) = DiffusiveFlux(nothing, κ) -function DiffusiveFlux(grid; κ = 1e-2) - Δz = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) - return DiffusiveFlux(Δz, κ) +function DiffusiveFlux(grid; κ = 0.1) + δ = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) + return DiffusiveFlux(δ, κ) end regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = From d2f3feaf33e7d9d930ceb5b1f5c8cd0b0a8e2a8d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 19:10:02 +0100 Subject: [PATCH 20/48] pass also gravity --- .../similarity_theory_turbulent_fluxes.jl | 6 ++-- .../CrossRealmFluxes/surface_temperature.jl | 35 +++++++++++-------- 2 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 2f7130ca..688c36c1 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -394,10 +394,10 @@ end cₚ = AtmosphericThermodynamics.cp_m(ℂ, 𝒬₁) # moist heat capacity ℰv = AtmosphericThermodynamics.latent_heat_vapor(ℂ, 𝒬₁) - θ₀ = retrieve_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, + θ₀ = retrieve_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, prescribed_heat_fluxes, radiative_properties) - + θ₁ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₁) # Temperature difference including the ``lapse rate'' `α = g / cₚ` @@ -437,7 +437,7 @@ end prescribed_heat_fluxes, radiative_properties, similarity_theory.bulk_velocity) - + # "initial" scales because we will recompute them u★ = estimated_characteristic_scales.momentum θ★ = estimated_characteristic_scales.temperature diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index c136aeba..cb9b9d6e 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -49,20 +49,26 @@ end regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = DiagnosticSurfaceTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ)) -# The balance solved is +# The flux balance could be solved either # -# Θo - Θₛ -# κ ------- = Qn (all fluxes positive upwards) -# δ +# Θₒ - Θₛⁿ⁺¹ +# κ ---------- = Qn (all fluxes positive upwards) +# δ # -# Where the LHS is the internal diffusive flux inside the ocean and the RHS are the -# atmospheric and radiative fluxes (provided explicitly and iterated upon). -@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = θo - Qn / F.κ * F.δ - -# Change 𝒬₀ as a function of incoming and outgoing fluxes. The flaw here is that -# the ocean emissivity and albedo are fixed, but they might be a function of the -# surface temperature, so we might need to pass actually the radiation and the -# albedo and emissivity as arguments. +# Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) +# and the RHS are the atmospheric and radiative fluxes are provided explicitly, or +# +# Θₒ - Θₛⁿ⁺¹ +# κ ---------- - σϵ θₛⁿ⁺¹θₛⁿ³ = Qn (all fluxes positive upwards) +# δ +# +# Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) +# plus the (semi-implicit) outgoing longwave flux and the RHS are the remaining atmospheric and radiative fluxes +# provided explicitly. +@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = (θo - Qn / F.κ * F.δ) + +# he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the +# surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. @inline function retrieve_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, prescribed_heat_fluxes, @@ -71,10 +77,9 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, Rd = prescribed_heat_fluxes # net downwelling radiation (positive out of the ocean) # upwelling radiation is calculated explicitly - # TODO: we could calculate it semi-implicitly as ϵσTⁿ⁺¹Tⁿ³ Ru = upwelling_radiation(θ₀, radiation_properties) Rn = Rd + Ru - + u★ = Σ★.momentum θ★ = Σ★.temperature q★ = Σ★.water_vapor @@ -83,7 +88,7 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, Qs = - ρₐ * u★ * (cₚ * θ★ + q★ * ℰv) # Net heat flux (positive out of the ocean) - Qn = (Qs + Rn) / ρₒ / cpₒ + Qn = (Qs + Rn) / ρₒ / cpₒ θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) From 552fd53580fe60ec88254cc3f189864ca925c781 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 08:30:03 +0100 Subject: [PATCH 21/48] change water vapor saturation --- .../atmosphere_ocean_fluxes.jl | 5 ++-- .../similarity_theory_turbulent_fluxes.jl | 28 ++++++++++++++----- .../CrossRealmFluxes/surface_temperature.jl | 28 ++++++++++--------- 3 files changed, 39 insertions(+), 22 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index 422bb9fb..efc655e1 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -334,10 +334,11 @@ end radiative_properties = local_radiation_properties(i, j, 1, grid, radiation, time) turbulent_fluxes, surface_temperature = compute_similarity_theory_fluxes(similarity_theory, - dynamic_ocean_state, - dynamic_atmos_state, + dynamic_ocean_state, + dynamic_atmos_state, prescribed_heat_fluxes, radiative_properties, + Sₒ, ocean_density, ocean_heat_capacity, atmosphere_boundary_layer_height, diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 688c36c1..170d0e28 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -195,7 +195,7 @@ function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature fields = (; latent_heat, sensible_heat, water_vapor, x_momentum, y_momentum, T_surface) - surface_temperature_type = regularize_surface_temperature(surface_temperature_type, grid) + surface_temperature_type = regularize_surface_temperature_type(surface_temperature_type, grid) return SimilarityTheoryTurbulentFluxes(eltype(grid); surface_temperature_type, kw..., fields) end @@ -241,6 +241,7 @@ struct COARELogarithmicSimilarityProfile end atmos_state, prescribed_heat_fluxes, # Possibly use in state_differences radiative_properties, + ocean_salinity, ocean_density, ocean_heat_capacity, atmos_boundary_layer_height, @@ -285,6 +286,7 @@ struct COARELogarithmicSimilarityProfile end similarity_theory, atmos_state, surface_state, + ocean_salinity, ocean_density, ocean_heat_capacity, atmos_boundary_layer_height, @@ -377,7 +379,10 @@ end @inline velocity_differences(𝒰₁, 𝒰₀, ::RelativeVelocity) = @inbounds 𝒰₁.u[1] - 𝒰₀.u[1], 𝒰₁.u[2] - 𝒰₀.u[2] @inline velocity_differences(𝒰₁, 𝒰₀, ::WindVelocity) = @inbounds 𝒰₁.u[1], 𝒰₁.u[2] -@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, Σ★, g, ρₒ, cpₒ, surface_temperature_type, +@inline function state_differences(ℂ, 𝒰₁, 𝒰₀, θ₀, S₀, Σ★, g, ρₒ, cpₒ, + water_mole_fraction, + water_vapor_saturation, + surface_temperature_type, prescribed_heat_fluxes, radiative_properties, bulk_velocity) @@ -394,17 +399,22 @@ end cₚ = AtmosphericThermodynamics.cp_m(ℂ, 𝒬₁) # moist heat capacity ℰv = AtmosphericThermodynamics.latent_heat_vapor(ℂ, 𝒬₁) - θ₀ = retrieve_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, - prescribed_heat_fluxes, - radiative_properties) + θ₀ = compute_surface_temperature(surface_temperature_type, θ₀, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, + prescribed_heat_fluxes, + radiative_properties) θ₁ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₁) # Temperature difference including the ``lapse rate'' `α = g / cₚ` Δθ = θ₁ - θ₀ + g / cₚ * Δh - + + T₀ = θ₀ - celsius_to_kelvin q₁ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₁) - q₀ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₀) + q₀ = seawater_saturation_specific_humidity(ℂ, T₀, S₀, 𝒬ₐ, + water_mole_fraction, + water_vapor_saturation, + AtmosphericThermodynamics.Liquid()) + Δq = q₁ - q₀ return Δh, Δu, Δv, Δθ, Δq, θ₀ @@ -416,6 +426,7 @@ end similarity_theory, atmos_state, surface_state, + surface_salinity, ocean_density, ocean_heat_capacity, atmos_boundary_layer_height, @@ -429,10 +440,13 @@ end atmos_state, surface_state, surface_temperature, + surface_salinity, estimated_characteristic_scales, gravitational_acceleration, ocean_density, ocean_heat_capacity, + similarity_theory.water_mole_fraction, + similarity_theory.water_vapor_saturation, similarity_theory.surface_temperature_type, prescribed_heat_fluxes, radiative_properties, diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index cb9b9d6e..a9d577cf 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -9,7 +9,7 @@ import Thermodynamics as AtmosphericThermodynamics @inline upwelling_radiation(θ₀, r) = r.σ * r.ϵ * θ₀^4 # For any surface temperture type that does not depend on the grid -regularize_surface_temperature(surface_temperature_type, grid) = surface_temperature_type +regularize_surface_temperature_type(surface_temperature_type, grid) = surface_temperature_type #### #### Prescribed surface temperature (the easiest case) @@ -18,7 +18,7 @@ regularize_surface_temperature(surface_temperature_type, grid) = surface_tempera struct PrescribedSurfaceTemperature end # Do nothing (just copy the temperature) -@inline retrieve_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ +@inline compute_surface_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ #### #### Diagnostic surface temperature calculated as a flux balance @@ -34,20 +34,22 @@ struct DiffusiveFlux{Z, K} end # A default constructor for DiagnosticSurfaceTemperature -function DiagnosticSurfaceTemperature(; κ = 0.1) - internal_flux = DiffusiveFlux(; κ) +function DiagnosticSurfaceTemperature(; κ = 0.1, δ = nothing) + internal_flux = DiffusiveFlux(; κ, δ) return DiagnosticSurfaceTemperature(internal_flux) end -DiffusiveFlux(; κ = 1e-2) = DiffusiveFlux(nothing, κ) +DiffusiveFlux(; κ = 1e-2, δ = nothing) = DiffusiveFlux(δ, κ) -function DiffusiveFlux(grid; κ = 0.1) - δ = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) +function DiffusiveFlux(grid; κ = 0.1, δ = nothing) + if isnothing(δ) + δ = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) + end return DiffusiveFlux(δ, κ) end -regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = - DiagnosticSurfaceTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ)) +regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = + DiagnosticSurfaceTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ, δ = T.internal_flux.δ)) # The flux balance could be solved either # @@ -69,10 +71,10 @@ regularize_surface_temperature(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. -@inline function retrieve_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, - ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, - prescribed_heat_fluxes, - radiation_properties) +@inline function compute_surface_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, + ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, + prescribed_heat_fluxes, + radiation_properties) Rd = prescribed_heat_fluxes # net downwelling radiation (positive out of the ocean) From 03c2849a23bc6f8242e3703925ac6aeed3d3f9a2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 08:33:17 +0100 Subject: [PATCH 22/48] change nomenclature --- .../CrossRealmFluxes/surface_temperature.jl | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index a9d577cf..dde422e7 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -67,11 +67,11 @@ regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveF # Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) # plus the (semi-implicit) outgoing longwave flux and the RHS are the remaining atmospheric and radiative fluxes # provided explicitly. -@inline flux_balance_temperature(F::DiffusiveFlux, θo, Qn) = (θo - Qn / F.κ * F.δ) +@inline flux_balance_temperature(F::DiffusiveFlux, θo, Jᵀ) = (θₒ - Jᵀ / F.κ * F.δ) # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. -@inline function compute_surface_temperature(st::DiagnosticSurfaceTemperature, θ₀, ℂ, 𝒬₀, +@inline function compute_surface_temperature(st::DiagnosticSurfaceTemperature, θₛ, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, prescribed_heat_fluxes, radiation_properties) @@ -79,20 +79,20 @@ regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveF Rd = prescribed_heat_fluxes # net downwelling radiation (positive out of the ocean) # upwelling radiation is calculated explicitly - Ru = upwelling_radiation(θ₀, radiation_properties) - Rn = Rd + Ru + Ru = upwelling_radiation(θₛ, radiation_properties) + Rn = Rd + Ru # Net radiation (positive out of the ocean) u★ = Σ★.momentum θ★ = Σ★.temperature q★ = Σ★.water_vapor - # sensible heat flux + latent heat flux (positive out of the ocean) - Qs = - ρₐ * u★ * (cₚ * θ★ + q★ * ℰv) + # Turbulent heat fluxes, sensible + latent (positive out of the ocean) + Qt = - ρₐ * u★ * (cₚ * θ★ + q★ * ℰv) - # Net heat flux (positive out of the ocean) - Qn = (Qs + Rn) / ρₒ / cpₒ + # Net temperature flux (positive out of the ocean) + Jᵀ = (Qt + Rn) / ρₒ / cpₒ - θo = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) + θₒ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) - return flux_balance_temperature(st.internal_flux, θo, Qn) + return flux_balance_temperature(st.internal_flux, θₒ, Jᵀ) end \ No newline at end of file From 6022a6feaf77aa2d73c5414604e82da29300fa0c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 08:34:32 +0100 Subject: [PATCH 23/48] finish changing nomenclature --- .../CrossRealmFluxes/surface_temperature.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index dde422e7..915caa59 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -54,20 +54,20 @@ regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveF # The flux balance could be solved either # # Θₒ - Θₛⁿ⁺¹ -# κ ---------- = Qn (all fluxes positive upwards) +# κ ---------- = Jᵀ (all fluxes positive upwards) # δ # # Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) # and the RHS are the atmospheric and radiative fluxes are provided explicitly, or # -# Θₒ - Θₛⁿ⁺¹ -# κ ---------- - σϵ θₛⁿ⁺¹θₛⁿ³ = Qn (all fluxes positive upwards) -# δ +# Θₒ - Θₛⁿ⁺¹ σϵ θₛⁿ⁺¹θₛⁿ³ +# κ ---------- - ------------ = Jᵀ (all fluxes positive upwards) +# δ ρₒ cpₒ # # Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) # plus the (semi-implicit) outgoing longwave flux and the RHS are the remaining atmospheric and radiative fluxes # provided explicitly. -@inline flux_balance_temperature(F::DiffusiveFlux, θo, Jᵀ) = (θₒ - Jᵀ / F.κ * F.δ) +@inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = (θₒ - Jᵀ / F.κ * F.δ) # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. From c261e368cbe9eebb24a114968083020b40493c6f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 08:35:24 +0100 Subject: [PATCH 24/48] Update src/OceanSimulations/OceanSimulations.jl Co-authored-by: Gregory L. Wagner --- src/OceanSimulations/OceanSimulations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index b9098120..41c8b2b2 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -61,7 +61,7 @@ end const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} function default_free_surface(grid::TripolarOfSomeKind; - fixed_Δt = infer_maximum_Δt(grid), + fixed_Δt = compute_maximum_Δt(grid), cfl = 0.7) free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) @info "Using a $(free_surface)" From b479f79453454423194f0293888927c47dc59fb9 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 09:59:45 +0100 Subject: [PATCH 25/48] correct temperature --- .../similarity_theory_turbulent_fluxes.jl | 8 ++------ .../CrossRealmFluxes/surface_temperature.jl | 7 ++++--- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 170d0e28..ba3c704a 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -21,7 +21,6 @@ import SurfaceFluxes.Parameters: thermodynamics_params, uf_params, von_karman_const, - universal_func_type, grav ##### @@ -53,8 +52,6 @@ const STTF = SimilarityTheoryTurbulentFluxes @inline grav(fluxes::STTF) = fluxes.gravitational_acceleration @inline molmass_ratio(fluxes::STTF) = molmass_ratio(fluxes.thermodynamics_parameters) -@inline universal_func_type(::STTF{<:Any, <:Any, <:BusingerParams}) = BusingerType() - Adapt.adapt_structure(to, fluxes::STTF) = SimilarityTheoryTurbulentFluxes(adapt(to, fluxes.gravitational_acceleration), adapt(to, fluxes.von_karman_constant), adapt(to, fluxes.turbulent_prandtl_number), @@ -407,10 +404,9 @@ end # Temperature difference including the ``lapse rate'' `α = g / cₚ` Δθ = θ₁ - θ₀ + g / cₚ * Δh - - T₀ = θ₀ - celsius_to_kelvin + q₁ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₁) - q₀ = seawater_saturation_specific_humidity(ℂ, T₀, S₀, 𝒬ₐ, + q₀ = seawater_saturation_specific_humidity(ℂ, θ₀, S₀, 𝒬₁, water_mole_fraction, water_vapor_saturation, AtmosphericThermodynamics.Liquid()) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 915caa59..172930fd 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -65,9 +65,9 @@ regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveF # δ ρₒ cpₒ # # Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) -# plus the (semi-implicit) outgoing longwave flux and the RHS are the remaining atmospheric and radiative fluxes +# plus the (semi-implicit) outgoing longwave radiation and the RHS are the remaining atmospheric and radiative fluxes # provided explicitly. -@inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = (θₒ - Jᵀ / F.κ * F.δ) +@inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = θₒ - Jᵀ / F.κ * F.δ # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. @@ -93,6 +93,7 @@ regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveF Jᵀ = (Qt + Rn) / ρₒ / cpₒ θₒ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬₀) + θₛ = flux_balance_temperature(st.internal_flux, θₒ, Jᵀ) - return flux_balance_temperature(st.internal_flux, θₒ, Jᵀ) + return θₛ end \ No newline at end of file From 846fda5db027e8cf493d607111802d2d61a3111f Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 10:31:58 +0100 Subject: [PATCH 26/48] try new one degree --- experiments/one_degree_simulation/one_degree_simulation.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 22fd8fec..3651d4ed 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -10,7 +10,7 @@ using CUDA: @allowscalar, device! using Oceananigans.Grids: znode -arch = CPU() +arch = GPU() ##### ##### Grid and Bathymetry @@ -48,7 +48,7 @@ gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1 catke = ClimaOcean.OceanSimulations.default_ocean_closure() viscous_closure = HorizontalScalarBiharmonicDiffusivity(ν = 1e11) -closure = (catke, viscous_closure) +closure = (gm, catke, viscous_closure) ##### ##### Restoring @@ -87,6 +87,8 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), ##### Atmospheric forcing ##### +surface_temperature_type = DiagnosticSurfaceTemperature(κ=0.01, δ=1.0) +similarity_theory = SimilarityTheoryTurbulentFluxes(grid, surface_temperature_type) radiation = Radiation(arch) atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) From c20b865f4e16680ec21a55bcf123e40a96266598 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 10:32:23 +0100 Subject: [PATCH 27/48] export --- experiments/one_degree_simulation/one_degree_simulation.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 3651d4ed..279d5ab4 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -3,6 +3,7 @@ using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting using OrthogonalSphericalShellGrids using Oceananigans using Oceananigans.Units +using Oceananigans.OceanSeaIceModels.CrossRealmFluxes: DiagnosticSurfaceTemperature, DiffusiveFlux using CFTime using Dates using Printf From 54226bfdfe444d08d2d7f435629cf173479ce7e2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 11:36:23 +0100 Subject: [PATCH 28/48] rebuilding the steps --- .../CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index ba3c704a..7f99c9ca 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -406,11 +406,16 @@ end Δθ = θ₁ - θ₀ + g / cₚ * Δh q₁ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬₁) + + # Recomputing the saturation specific humidity at the surface based on the new temperature q₀ = seawater_saturation_specific_humidity(ℂ, θ₀, S₀, 𝒬₁, water_mole_fraction, water_vapor_saturation, AtmosphericThermodynamics.Liquid()) + 𝒬ₛ = AtmosphericThermodynamics.PhaseEquil_pTq(ℂ, 𝒬₀.p, θ₀, q₀) + q₀ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬ₛ) + Δq = q₁ - q₀ return Δh, Δu, Δv, Δθ, Δq, θ₀ From 6848c54dd9978aec68a3455e51edf271ab942dc3 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 11:23:24 +0100 Subject: [PATCH 29/48] change names to the types --- examples/generate_surface_fluxes.jl | 74 +++++++++++++------ .../one_degree_simulation.jl | 4 +- .../os_papa_surface_temperature.jl | 6 +- .../similarity_theory_turbulent_fluxes.jl | 6 +- .../CrossRealmFluxes/surface_temperature.jl | 18 ++--- test/test_surface_fluxes.jl | 5 +- 6 files changed, 73 insertions(+), 40 deletions(-) diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 946a8be5..25aa19e0 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -1,4 +1,4 @@ -# # Surface fluxes from prescribed ocean and atmosphere +# Surface fluxes from prescribed ocean and atmosphere # # ClimaOcean uses bulk formulae to estimate the surface exchange of momentum, # heat, and water vapor between the atmosphere and the ocean. @@ -7,7 +7,7 @@ # using ECCO2 data for the ocean and JRA55 data for the atmosphere. # # For this example, we need ClimaOcean with its DataWrangling modules: ECCO2 and JRA55. -# We also need Oceananigans for the ImmersedBoundaryGrid and Field utilities, and CairoMakie to plot. +# # We also need Oceananigans for the ImmersedBoundaryGrid and Field utilities, and CairoMakie to plot. using ClimaOcean using ClimaOcean.ECCO @@ -45,7 +45,8 @@ save("ECCO_continents.png", fig) #hide # January 1st (at 00:00 AM and 03:00 AM). atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory()) -ocean = ocean_simulation(grid) +ocean1 = ocean_simulation(grid) +ocean2 = ocean_simulation(grid) # Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity # to the ECCO2 data by first creating T, S metadata objects, @@ -56,37 +57,66 @@ S_metadata = ECCOMetadata(:salinity) # Note that if a date is not provided to `ECCOMetadata`, then the default Jan 1st, 1992 is used. # To copy the ECCO state into `ocean.model`, we use `set!`, -set!(ocean.model; T=T_metadata, S=S_metadata) +set!(ocean1.model; T=T_metadata, S=S_metadata) +set!(ocean2.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. -coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation=Radiation()) +using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: SkinTemperature -# Now that the surface fluxes are computed, we can extract and visualize them. -# The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`. +similarity_theory = SimilarityTheoryTurbulentFluxes(grid, surface_temperature_type = SkinTemperature(κ = 0.01, δ = 1.0)) +coupled_model1 = OceanSeaIceModel(ocean1; atmosphere, radiation=Radiation()) +coupled_model2 = OceanSeaIceModel(ocean2; atmosphere, radiation=Radiation(), similarity_theory) -fluxes = coupled_model.fluxes.turbulent.fields -λ, φ, z = nodes(fluxes.sensible_heat) +Qnet1 = coupled_model1.ocean.model.tracers.T.boundary_conditions.top.condition; +Qnet2 = coupled_model2.ocean.model.tracers.T.boundary_conditions.top.condition; -fig = Figure(size = (800, 800), fontsize = 15) +Ts1 = coupled_model1.fluxes.turbulent.fields.T_surface +Ts2 = coupled_model2.fluxes.turbulent.fields.T_surface -ax = Axis(fig[1, 1], title = "Sensible heat flux (W m⁻²)", ylabel = "Latitude") -heatmap!(ax, λ, φ, interior(fluxes.sensible_heat, :, :, 1); colormap = :bwr) +fig = Figure() +ax = Axis(fig[1, 1]) +hm = heatmap!((interior(Qnet1, :, :, 1) .- interior(Qnet2, :, :, 1)) .* 1020 .* coupled_model1.fluxes.ocean_heat_capacity, colormap = :balance, colorrange = (-2, 2)) +Colorbar(fig[1, 2], hm) +ax = Axis(fig[1, 3]) +hm = heatmap!((interior(Ts1, :, :, 1) .- interior(Ts2, :, :, 1)), colormap = :balance, colorrange = (-0.2, 0.2)) +Colorbar(fig[1, 4], hm) + +fig2 = Figure() +ax = Axis(fig2[1, 1], title = "iterations for Prescribed temperature") +hm = heatmap!(ax, interior(coupled_model1.fluxes.turbulent.fields.x_momentum, :, :, 1), colorrange = (-1, 1)) #, colormap = :balance) +ax = Axis(fig2[1, 2], title = "iterations for Diagnostic temperature") +hm = heatmap!(ax, interior(coupled_model2.fluxes.turbulent.fields.x_momentum, :, :, 1), colorrange = (-1, 1)) #, colormap = :balance) +Colorbar(fig2[1, 3], hm) +ax = Axis(fig2[1, 4], title = "Difference in iterations (Prescribed - Diagnostic)") +hm = heatmap!(ax, interior(coupled_model1.fluxes.turbulent.fields.x_momentum, :, :, 1) .- interior(coupled_model2.fluxes.turbulent.fields.x_momentum, :, :, 1), colorrange = (-1e-2, 1e-2), colormap = :balance) +Colorbar(fig2[1, 5], hm) + +# # Now that the surface fluxes are computed, we can extract and visualize them. +# # The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`. + +# fluxes = coupled_model.fluxes.turbulent.fields +# λ, φ, z = nodes(fluxes.sensible_heat) + +# fig = Figure(size = (800, 800), fontsize = 15) + +# ax = Axis(fig[1, 1], title = "Sensible heat flux (W m⁻²)", ylabel = "Latitude") +# heatmap!(ax, λ, φ, interior(fluxes.sensible_heat, :, :, 1); colormap = :bwr) -ax = Axis(fig[1, 2], title = "Latent heat flux (W m⁻²)") -heatmap!(ax, λ, φ, interior(fluxes.latent_heat, :, :, 1); colormap = :bwr) +# ax = Axis(fig[1, 2], title = "Latent heat flux (W m⁻²)") +# heatmap!(ax, λ, φ, interior(fluxes.latent_heat, :, :, 1); colormap = :bwr) -ax = Axis(fig[2, 1], title = "Zonal wind stress (N m)", ylabel = "Latitude") -heatmap!(ax, λ, φ, interior(fluxes.x_momentum, :, :, 1); colormap = :bwr) +# ax = Axis(fig[2, 1], title = "Zonal wind stress (N m)", ylabel = "Latitude") +# heatmap!(ax, λ, φ, interior(fluxes.x_momentum, :, :, 1); colormap = :bwr) -ax = Axis(fig[2, 2], title = "Meridional wind stress (N m)", xlabel = "Longitude") -heatmap!(ax, λ, φ, interior(fluxes.y_momentum, :, :, 1); colormap = :bwr) +# ax = Axis(fig[2, 2], title = "Meridional wind stress (N m)", xlabel = "Longitude") +# heatmap!(ax, λ, φ, interior(fluxes.y_momentum, :, :, 1); colormap = :bwr) -ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻² s⁻¹)", xlabel = "Longitude", ylabel = "Latitude") -heatmap!(ax, λ, φ, interior(fluxes.water_vapor, :, :, 1); colormap = :bwr) +# ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻² s⁻¹)", xlabel = "Longitude", ylabel = "Latitude") +# heatmap!(ax, λ, φ, interior(fluxes.water_vapor, :, :, 1); colormap = :bwr) -save("fluxes.png", fig) -# ![](fluxes.png) +# save("fluxes.png", fig) +# # ![](fluxes.png) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 279d5ab4..d61c48e9 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -3,7 +3,7 @@ using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting using OrthogonalSphericalShellGrids using Oceananigans using Oceananigans.Units -using Oceananigans.OceanSeaIceModels.CrossRealmFluxes: DiagnosticSurfaceTemperature, DiffusiveFlux +using Oceananigans.OceanSeaIceModels.CrossRealmFluxes: SkinTemperature, DiffusiveFlux using CFTime using Dates using Printf @@ -88,7 +88,7 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), ##### Atmospheric forcing ##### -surface_temperature_type = DiagnosticSurfaceTemperature(κ=0.01, δ=1.0) +surface_temperature_type = SkinTemperature(κ=0.01, δ=1.0) similarity_theory = SimilarityTheoryTurbulentFluxes(grid, surface_temperature_type) radiation = Radiation(arch) atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) diff --git a/experiments/single_column_experiments/os_papa_surface_temperature.jl b/experiments/single_column_experiments/os_papa_surface_temperature.jl index cf428d82..2bed91d2 100644 --- a/experiments/single_column_experiments/os_papa_surface_temperature.jl +++ b/experiments/single_column_experiments/os_papa_surface_temperature.jl @@ -1,8 +1,8 @@ # This simulation tests the difference between using a -# `PrescribedSurfaceTemperature` and a `DiagnosticSurfaceTemperature` +# `BulkTemperature` and a `SkinTemperature` using ClimaOcean -using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: SimilarityTheoryTurbulentFluxes, DiagnosticSurfaceTemperature +using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: SimilarityTheoryTurbulentFluxes, SkinTemperature using Oceananigans using Oceananigans.Units using Oceananigans.BuoyancyModels: buoyancy_frequency @@ -41,7 +41,7 @@ atmosphere = JRA55_prescribed_atmosphere(1:last_time; radiation = Radiation() coupled_model_prescribed = OceanSeaIceModel(ocean1; atmosphere, radiation) -similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = DiagnosticSurfaceTemperature(; κ = 0.5)) +similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = SkinTemperature(; κ = 0.5)) coupled_model_diagnostic = OceanSeaIceModel(ocean2; atmosphere, radiation, similarity_theory) for (coupled_model, suffix) in zip([coupled_model_diagnostic, coupled_model_prescribed], diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 7f99c9ca..52b6e461 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -114,7 +114,7 @@ struct RelativeVelocity end water_mole_fraction = convert(FT, 0.98), roughness_lengths = default_roughness_lengths(FT), similarity_profile_type = LogarithmicSimilarityProfile(), - surface_temperature_type = PrescribedSurfaceTemperature(), + surface_temperature_type = BulkTemperature(), bulk_velocity = RelativeVelocity(), tolerance = 1e-8, maxiter = 100, @@ -159,7 +159,7 @@ function SimilarityTheoryTurbulentFluxes(FT::DataType = Float64; water_mole_fraction = convert(FT, 0.98), roughness_lengths = default_roughness_lengths(FT), similarity_profile_type = LogarithmicSimilarityProfile(), - surface_temperature_type = PrescribedSurfaceTemperature(), + surface_temperature_type = BulkTemperature(), bulk_velocity = RelativeVelocity(), tolerance = 1e-8, maxiter = 100, @@ -182,7 +182,7 @@ function SimilarityTheoryTurbulentFluxes(FT::DataType = Float64; fields) end -function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature_type = PrescribedSurfaceTemperature(), kw...) +function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature_type = BulkTemperature(), kw...) water_vapor = Field{Center, Center, Nothing}(grid) latent_heat = Field{Center, Center, Nothing}(grid) sensible_heat = Field{Center, Center, Nothing}(grid) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 172930fd..9d1c41e6 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -15,16 +15,16 @@ regularize_surface_temperature_type(surface_temperature_type, grid) = surface_te #### Prescribed surface temperature (the easiest case) #### -struct PrescribedSurfaceTemperature end +struct BulkTemperature end # Do nothing (just copy the temperature) -@inline compute_surface_temperature(::PrescribedSurfaceTemperature, θ₀, args...) = θ₀ +@inline compute_surface_temperature(::BulkTemperature, θ₀, args...) = θ₀ #### #### Diagnostic surface temperature calculated as a flux balance #### -struct DiagnosticSurfaceTemperature{I} +struct SkinTemperature{I} internal_flux :: I end @@ -33,10 +33,10 @@ struct DiffusiveFlux{Z, K} κ :: K # diffusivity in m²s⁻¹ end -# A default constructor for DiagnosticSurfaceTemperature -function DiagnosticSurfaceTemperature(; κ = 0.1, δ = nothing) +# A default constructor for SkinTemperature +function SkinTemperature(; κ = 0.1, δ = nothing) internal_flux = DiffusiveFlux(; κ, δ) - return DiagnosticSurfaceTemperature(internal_flux) + return SkinTemperature(internal_flux) end DiffusiveFlux(; κ = 1e-2, δ = nothing) = DiffusiveFlux(δ, κ) @@ -48,8 +48,8 @@ function DiffusiveFlux(grid; κ = 0.1, δ = nothing) return DiffusiveFlux(δ, κ) end -regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveFlux}, grid) = - DiagnosticSurfaceTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ, δ = T.internal_flux.δ)) +regularize_surface_temperature_type(T::SkinTemperature{<:DiffusiveFlux}, grid) = + SkinTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ, δ = T.internal_flux.δ)) # The flux balance could be solved either # @@ -71,7 +71,7 @@ regularize_surface_temperature_type(T::DiagnosticSurfaceTemperature{<:DiffusiveF # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. -@inline function compute_surface_temperature(st::DiagnosticSurfaceTemperature, θₛ, ℂ, 𝒬₀, +@inline function compute_surface_temperature(st::SkinTemperature, θₛ, ℂ, 𝒬₀, ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, prescribed_heat_fluxes, radiation_properties) diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index c401cce7..32ce6afb 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -200,6 +200,10 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, end end +@testset "SkinTemperature" begin + +end + @testset "Fluxes regression" begin for arch in test_architectures @info "Testing fluxes regression..." @@ -272,4 +276,3 @@ end end end - From ae1b6a23c3bca3989037737183330cbf69ab39b0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 11:29:24 +0100 Subject: [PATCH 30/48] add a comment --- .../CrossRealmFluxes/surface_temperature.jl | 24 +++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 9d1c41e6..054a9e9e 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -12,18 +12,38 @@ import Thermodynamics as AtmosphericThermodynamics regularize_surface_temperature_type(surface_temperature_type, grid) = surface_temperature_type #### -#### Prescribed surface temperature (the easiest case) +#### Bulk surface temperature (the easiest case) #### +""" + struct BulkTemperature end + +A type to represent the surface temperature used in the flux calculation. +The surface temperature is not calculated but provided by either the ocean or the sea ice model. +""" struct BulkTemperature end # Do nothing (just copy the temperature) @inline compute_surface_temperature(::BulkTemperature, θ₀, args...) = θ₀ #### -#### Diagnostic surface temperature calculated as a flux balance +#### Skin surface temperature calculated as a flux balance #### +""" + struct SkinTemperature + internal_flux :: I + end + +A type to represent the surface temperature used in the flux calculation. +The surface temperature is calculated from the flux balance at the surface. +In particular, the surface temperature ``θₛ`` is the root of: + +F(θₛ) - Jᵀ = 0 (all fluxes positive upwards) + +where Jᵀ are the fluxes at the top of the surface (turbulent + radiative), and F is the internal diffusive flux +dependent on the surface temperature itself. +""" struct SkinTemperature{I} internal_flux :: I end From a0140ffb952f0d314c270e4711bebdbb5b2ecca3 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 11:30:58 +0100 Subject: [PATCH 31/48] revert unrelated files --- examples/generate_surface_fluxes.jl | 74 ++++++------------- .../one_degree_simulation.jl | 13 +--- 2 files changed, 25 insertions(+), 62 deletions(-) diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 25aa19e0..946a8be5 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -1,4 +1,4 @@ -# Surface fluxes from prescribed ocean and atmosphere +# # Surface fluxes from prescribed ocean and atmosphere # # ClimaOcean uses bulk formulae to estimate the surface exchange of momentum, # heat, and water vapor between the atmosphere and the ocean. @@ -7,7 +7,7 @@ # using ECCO2 data for the ocean and JRA55 data for the atmosphere. # # For this example, we need ClimaOcean with its DataWrangling modules: ECCO2 and JRA55. -# # We also need Oceananigans for the ImmersedBoundaryGrid and Field utilities, and CairoMakie to plot. +# We also need Oceananigans for the ImmersedBoundaryGrid and Field utilities, and CairoMakie to plot. using ClimaOcean using ClimaOcean.ECCO @@ -45,8 +45,7 @@ save("ECCO_continents.png", fig) #hide # January 1st (at 00:00 AM and 03:00 AM). atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory()) -ocean1 = ocean_simulation(grid) -ocean2 = ocean_simulation(grid) +ocean = ocean_simulation(grid) # Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity # to the ECCO2 data by first creating T, S metadata objects, @@ -57,66 +56,37 @@ S_metadata = ECCOMetadata(:salinity) # Note that if a date is not provided to `ECCOMetadata`, then the default Jan 1st, 1992 is used. # To copy the ECCO state into `ocean.model`, we use `set!`, -set!(ocean1.model; T=T_metadata, S=S_metadata) -set!(ocean2.model; T=T_metadata, S=S_metadata) +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. -using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: SkinTemperature +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation=Radiation()) -similarity_theory = SimilarityTheoryTurbulentFluxes(grid, surface_temperature_type = SkinTemperature(κ = 0.01, δ = 1.0)) -coupled_model1 = OceanSeaIceModel(ocean1; atmosphere, radiation=Radiation()) -coupled_model2 = OceanSeaIceModel(ocean2; atmosphere, radiation=Radiation(), similarity_theory) +# Now that the surface fluxes are computed, we can extract and visualize them. +# The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`. -Qnet1 = coupled_model1.ocean.model.tracers.T.boundary_conditions.top.condition; -Qnet2 = coupled_model2.ocean.model.tracers.T.boundary_conditions.top.condition; +fluxes = coupled_model.fluxes.turbulent.fields +λ, φ, z = nodes(fluxes.sensible_heat) -Ts1 = coupled_model1.fluxes.turbulent.fields.T_surface -Ts2 = coupled_model2.fluxes.turbulent.fields.T_surface +fig = Figure(size = (800, 800), fontsize = 15) -fig = Figure() -ax = Axis(fig[1, 1]) -hm = heatmap!((interior(Qnet1, :, :, 1) .- interior(Qnet2, :, :, 1)) .* 1020 .* coupled_model1.fluxes.ocean_heat_capacity, colormap = :balance, colorrange = (-2, 2)) -Colorbar(fig[1, 2], hm) -ax = Axis(fig[1, 3]) -hm = heatmap!((interior(Ts1, :, :, 1) .- interior(Ts2, :, :, 1)), colormap = :balance, colorrange = (-0.2, 0.2)) -Colorbar(fig[1, 4], hm) - -fig2 = Figure() -ax = Axis(fig2[1, 1], title = "iterations for Prescribed temperature") -hm = heatmap!(ax, interior(coupled_model1.fluxes.turbulent.fields.x_momentum, :, :, 1), colorrange = (-1, 1)) #, colormap = :balance) -ax = Axis(fig2[1, 2], title = "iterations for Diagnostic temperature") -hm = heatmap!(ax, interior(coupled_model2.fluxes.turbulent.fields.x_momentum, :, :, 1), colorrange = (-1, 1)) #, colormap = :balance) -Colorbar(fig2[1, 3], hm) -ax = Axis(fig2[1, 4], title = "Difference in iterations (Prescribed - Diagnostic)") -hm = heatmap!(ax, interior(coupled_model1.fluxes.turbulent.fields.x_momentum, :, :, 1) .- interior(coupled_model2.fluxes.turbulent.fields.x_momentum, :, :, 1), colorrange = (-1e-2, 1e-2), colormap = :balance) -Colorbar(fig2[1, 5], hm) - -# # Now that the surface fluxes are computed, we can extract and visualize them. -# # The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`. - -# fluxes = coupled_model.fluxes.turbulent.fields -# λ, φ, z = nodes(fluxes.sensible_heat) - -# fig = Figure(size = (800, 800), fontsize = 15) - -# ax = Axis(fig[1, 1], title = "Sensible heat flux (W m⁻²)", ylabel = "Latitude") -# heatmap!(ax, λ, φ, interior(fluxes.sensible_heat, :, :, 1); colormap = :bwr) +ax = Axis(fig[1, 1], title = "Sensible heat flux (W m⁻²)", ylabel = "Latitude") +heatmap!(ax, λ, φ, interior(fluxes.sensible_heat, :, :, 1); colormap = :bwr) -# ax = Axis(fig[1, 2], title = "Latent heat flux (W m⁻²)") -# heatmap!(ax, λ, φ, interior(fluxes.latent_heat, :, :, 1); colormap = :bwr) +ax = Axis(fig[1, 2], title = "Latent heat flux (W m⁻²)") +heatmap!(ax, λ, φ, interior(fluxes.latent_heat, :, :, 1); colormap = :bwr) -# ax = Axis(fig[2, 1], title = "Zonal wind stress (N m)", ylabel = "Latitude") -# heatmap!(ax, λ, φ, interior(fluxes.x_momentum, :, :, 1); colormap = :bwr) +ax = Axis(fig[2, 1], title = "Zonal wind stress (N m)", ylabel = "Latitude") +heatmap!(ax, λ, φ, interior(fluxes.x_momentum, :, :, 1); colormap = :bwr) -# ax = Axis(fig[2, 2], title = "Meridional wind stress (N m)", xlabel = "Longitude") -# heatmap!(ax, λ, φ, interior(fluxes.y_momentum, :, :, 1); colormap = :bwr) +ax = Axis(fig[2, 2], title = "Meridional wind stress (N m)", xlabel = "Longitude") +heatmap!(ax, λ, φ, interior(fluxes.y_momentum, :, :, 1); colormap = :bwr) -# ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻² s⁻¹)", xlabel = "Longitude", ylabel = "Latitude") -# heatmap!(ax, λ, φ, interior(fluxes.water_vapor, :, :, 1); colormap = :bwr) +ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻² s⁻¹)", xlabel = "Longitude", ylabel = "Latitude") +heatmap!(ax, λ, φ, interior(fluxes.water_vapor, :, :, 1); colormap = :bwr) -# save("fluxes.png", fig) -# # ![](fluxes.png) +save("fluxes.png", fig) +# ![](fluxes.png) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index d61c48e9..08577fed 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -3,7 +3,6 @@ using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting using OrthogonalSphericalShellGrids using Oceananigans using Oceananigans.Units -using Oceananigans.OceanSeaIceModels.CrossRealmFluxes: SkinTemperature, DiffusiveFlux using CFTime using Dates using Printf @@ -19,7 +18,7 @@ arch = GPU() Nx = 360 Ny = 180 -Nz = 50 +Nz = 100 z_faces = exponential_z_faces(; Nz, depth=5000, h=34) @@ -47,7 +46,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) catke = ClimaOcean.OceanSimulations.default_ocean_closure() -viscous_closure = HorizontalScalarBiharmonicDiffusivity(ν = 1e11) +viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000) closure = (gm, catke, viscous_closure) @@ -88,8 +87,6 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), ##### Atmospheric forcing ##### -surface_temperature_type = SkinTemperature(κ=0.01, δ=1.0) -similarity_theory = SimilarityTheoryTurbulentFluxes(grid, surface_temperature_type) radiation = Radiation(arch) atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) @@ -98,7 +95,7 @@ atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) ##### coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=2minutes, stop_time=30days) +simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days) ##### ##### Run it! @@ -130,7 +127,3 @@ end add_callback!(simulation, progress, IterationInterval(10)) run!(simulation) - -simulation.Δt = 25minutes - -run!(simulation) From 3ddb3ad344ac39a30a3725fa5ed1086946bfef96 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 11:33:06 +0100 Subject: [PATCH 32/48] lets not regularize... --- .../similarity_theory_turbulent_fluxes.jl | 2 -- .../CrossRealmFluxes/surface_temperature.jl | 17 ++--------------- 2 files changed, 2 insertions(+), 17 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 52b6e461..e0887f89 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -192,8 +192,6 @@ function SimilarityTheoryTurbulentFluxes(grid::AbstractGrid; surface_temperature fields = (; latent_heat, sensible_heat, water_vapor, x_momentum, y_momentum, T_surface) - surface_temperature_type = regularize_surface_temperature_type(surface_temperature_type, grid) - return SimilarityTheoryTurbulentFluxes(eltype(grid); surface_temperature_type, kw..., fields) end diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 054a9e9e..553be35a 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -8,9 +8,6 @@ import Thermodynamics as AtmosphericThermodynamics @inline upwelling_radiation(θ₀, ::Nothing) = zero(θ₀) @inline upwelling_radiation(θ₀, r) = r.σ * r.ϵ * θ₀^4 -# For any surface temperture type that does not depend on the grid -regularize_surface_temperature_type(surface_temperature_type, grid) = surface_temperature_type - #### #### Bulk surface temperature (the easiest case) #### @@ -54,22 +51,12 @@ struct DiffusiveFlux{Z, K} end # A default constructor for SkinTemperature -function SkinTemperature(; κ = 0.1, δ = nothing) +function SkinTemperature(; κ = 0.1, δ = 1.0) internal_flux = DiffusiveFlux(; κ, δ) return SkinTemperature(internal_flux) end -DiffusiveFlux(; κ = 1e-2, δ = nothing) = DiffusiveFlux(δ, κ) - -function DiffusiveFlux(grid; κ = 0.1, δ = nothing) - if isnothing(δ) - δ = @allowscalar Δzᶜᶜᶜ(1, 1, grid.Nz, grid) - end - return DiffusiveFlux(δ, κ) -end - -regularize_surface_temperature_type(T::SkinTemperature{<:DiffusiveFlux}, grid) = - SkinTemperature(DiffusiveFlux(grid; κ = T.internal_flux.κ, δ = T.internal_flux.δ)) +DiffusiveFlux(; κ = 1e-2, δ = 1.0) = DiffusiveFlux(δ, κ) # The flux balance could be solved either # From 5625be1eaac1eae6a149f997f794423c3ed97d51 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 12:05:25 +0100 Subject: [PATCH 33/48] move to examples --- examples/one_degree_global_simulation.jl | 133 +++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 examples/one_degree_global_simulation.jl diff --git a/examples/one_degree_global_simulation.jl b/examples/one_degree_global_simulation.jl new file mode 100644 index 00000000..f5f30d86 --- /dev/null +++ b/examples/one_degree_global_simulation.jl @@ -0,0 +1,133 @@ +using ClimaOcean +using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting +using OrthogonalSphericalShellGrids +using Oceananigans +using Oceananigans.Units +using CFTime +using Dates +using Printf +using CUDA: @allowscalar, device! + +using Oceananigans.Grids: znode + +arch = GPU() + +##### +##### Grid and Bathymetry +##### + +Nx = 360 +Ny = 180 +Nz = 100 + +z_faces = exponential_z_faces(; Nz, depth=5000, h=34) + +underlying_grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + first_pole_longitude = 70, + north_poles_latitude = 55) + +bottom_height = regrid_bathymetry(underlying_grid; + minimum_depth = 10, + interpolation_passes = 75, + major_basins = 2) + +# Open Gibraltar strait +# TODO: find a better way to do this +tampered_bottom_height = deepcopy(bottom_height) +view(tampered_bottom_height, 102:103, 124, 1) .= -400 + +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_height)) + +##### +##### Closures +##### + +gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) +catke = ClimaOcean.OceanSimulations.default_ocean_closure() +viscous_closure = HorizonalScalarDiffusivity(κ=1000, ν=2000) + +closure = (gm, catke, viscous_closure) + +##### +##### Restoring +##### + +restoring_rate = 1 / 2days +z_below_surface = @allowscalar znode(1, 1, grid.Nz, grid, Center(), Center(), Face()) + +mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0)) + +dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 11, 1) +temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./") +salinity = ECCOMetadata(:salinity; dates, version=ECCO4Monthly(), dir="./") + +# inpainting = NearestNeighborInpainting(30) should be enough to fill the gaps near bathymetry +FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50)) +FS = ECCORestoring(arch, salinity; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50)) +forcing = (T=FT, S=FS) + +##### +##### Ocean simulation +##### + +momentum_advection = VectorInvariant() +tracer_advection = Centered(order=2) + +# Should we add a side drag since this is at a coarser resolution? +ocean = ocean_simulation(grid; momentum_advection, tracer_advection, + closure, forcing, + tracers = (:T, :S, :e)) + +set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), + S=ECCOMetadata(:salinity; dates=first(dates))) + +##### +##### Atmospheric forcing +##### + +radiation = Radiation(arch) +atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) + +##### +##### Coupled simulation +##### + +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) +simulation = Simulation(coupled_model; Δt=2minutes, stop_time=30days) + +##### +##### Run it! +##### + +wall_time = Ref(time_ns()) + +function progress(sim) + ocean = sim.model.ocean + u, v, w = ocean.model.velocities + T = ocean.model.tracers.T + Tmax = maximum(interior(T)) + Tmin = minimum(interior(T)) + umax = (maximum(abs, interior(u)), + maximum(abs, interior(v)), + maximum(abs, interior(w))) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s \n", + prettytime(sim), iteration(sim), prettytime(sim.Δt), + umax..., Tmax, Tmin, prettytime(step_time)) + + wall_time[] = time_ns() + + return nothing +end + +add_callback!(simulation, progress, IterationInterval(10)) + +run!(simulation) + +simulation.Δt = 25minutes + +run!(simulation) From a3276060f93950caa9efd45f4004992f5c712e51 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 12:07:43 +0100 Subject: [PATCH 34/48] not in this PR --- examples/one_degree_global_simulation.jl | 133 ----------------------- 1 file changed, 133 deletions(-) delete mode 100644 examples/one_degree_global_simulation.jl diff --git a/examples/one_degree_global_simulation.jl b/examples/one_degree_global_simulation.jl deleted file mode 100644 index f5f30d86..00000000 --- a/examples/one_degree_global_simulation.jl +++ /dev/null @@ -1,133 +0,0 @@ -using ClimaOcean -using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting -using OrthogonalSphericalShellGrids -using Oceananigans -using Oceananigans.Units -using CFTime -using Dates -using Printf -using CUDA: @allowscalar, device! - -using Oceananigans.Grids: znode - -arch = GPU() - -##### -##### Grid and Bathymetry -##### - -Nx = 360 -Ny = 180 -Nz = 100 - -z_faces = exponential_z_faces(; Nz, depth=5000, h=34) - -underlying_grid = TripolarGrid(arch; - size = (Nx, Ny, Nz), - z = z_faces, - first_pole_longitude = 70, - north_poles_latitude = 55) - -bottom_height = regrid_bathymetry(underlying_grid; - minimum_depth = 10, - interpolation_passes = 75, - major_basins = 2) - -# Open Gibraltar strait -# TODO: find a better way to do this -tampered_bottom_height = deepcopy(bottom_height) -view(tampered_bottom_height, 102:103, 124, 1) .= -400 - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_height)) - -##### -##### Closures -##### - -gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) -catke = ClimaOcean.OceanSimulations.default_ocean_closure() -viscous_closure = HorizonalScalarDiffusivity(κ=1000, ν=2000) - -closure = (gm, catke, viscous_closure) - -##### -##### Restoring -##### - -restoring_rate = 1 / 2days -z_below_surface = @allowscalar znode(1, 1, grid.Nz, grid, Center(), Center(), Face()) - -mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0)) - -dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 11, 1) -temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./") -salinity = ECCOMetadata(:salinity; dates, version=ECCO4Monthly(), dir="./") - -# inpainting = NearestNeighborInpainting(30) should be enough to fill the gaps near bathymetry -FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50)) -FS = ECCORestoring(arch, salinity; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50)) -forcing = (T=FT, S=FS) - -##### -##### Ocean simulation -##### - -momentum_advection = VectorInvariant() -tracer_advection = Centered(order=2) - -# Should we add a side drag since this is at a coarser resolution? -ocean = ocean_simulation(grid; momentum_advection, tracer_advection, - closure, forcing, - tracers = (:T, :S, :e)) - -set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)), - S=ECCOMetadata(:salinity; dates=first(dates))) - -##### -##### Atmospheric forcing -##### - -radiation = Radiation(arch) -atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) - -##### -##### Coupled simulation -##### - -coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=2minutes, stop_time=30days) - -##### -##### Run it! -##### - -wall_time = Ref(time_ns()) - -function progress(sim) - ocean = sim.model.ocean - u, v, w = ocean.model.velocities - T = ocean.model.tracers.T - Tmax = maximum(interior(T)) - Tmin = minimum(interior(T)) - umax = (maximum(abs, interior(u)), - maximum(abs, interior(v)), - maximum(abs, interior(w))) - - step_time = 1e-9 * (time_ns() - wall_time[]) - - @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s \n", - prettytime(sim), iteration(sim), prettytime(sim.Δt), - umax..., Tmax, Tmin, prettytime(step_time)) - - wall_time[] = time_ns() - - return nothing -end - -add_callback!(simulation, progress, IterationInterval(10)) - -run!(simulation) - -simulation.Δt = 25minutes - -run!(simulation) From 1d169f7c4d65e5fa8445239cc68eeedb59603001 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 12:14:51 +0100 Subject: [PATCH 35/48] convert to FT --- .../CrossRealmFluxes/surface_temperature.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 553be35a..91ac8997 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -51,12 +51,12 @@ struct DiffusiveFlux{Z, K} end # A default constructor for SkinTemperature -function SkinTemperature(; κ = 0.1, δ = 1.0) - internal_flux = DiffusiveFlux(; κ, δ) +function SkinTemperature(FT::DataType=Float64; κ=1e-2, δ=1.0) + internal_flux = DiffusiveFlux(FT; κ, δ) return SkinTemperature(internal_flux) end -DiffusiveFlux(; κ = 1e-2, δ = 1.0) = DiffusiveFlux(δ, κ) +DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert(FT, κ)) # The flux balance could be solved either # From ab7b78b7f4e6757bac3a2bbffc8621833f1af6ff Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:03:41 +0100 Subject: [PATCH 36/48] another bugfix --- .../CrossRealmFluxes/atmosphere_ocean_fluxes.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index efc655e1..db433c72 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -329,9 +329,10 @@ end inactive = inactive_node(i, j, kᴺ, grid, c, c, c) maxiter = ifelse(inactive, 1, similarity_theory.maxiter) + time = Time(clock.time) prescribed_heat_fluxes = net_downwelling_radiation(i, j, grid, time, radiation, Rs, Rℓ) - radiative_properties = local_radiation_properties(i, j, 1, grid, radiation, time) + radiative_properties = local_radiation_properties(i, j, kᴺ, grid, radiation, time) turbulent_fluxes, surface_temperature = compute_similarity_theory_fluxes(similarity_theory, dynamic_ocean_state, From 84ec37b974e0825829261677d95deac7614c51cd Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:05:05 +0100 Subject: [PATCH 37/48] rearrange signature --- .../CrossRealmFluxes/atmosphere_ocean_fluxes.jl | 8 ++++---- src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index db433c72..b80ad76b 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -277,8 +277,9 @@ end atmos_thermodynamics_parameters) i, j = @index(Global, NTuple) - kᴺ = size(grid, 3) # index of the top ocean cell - + kᴺ = size(grid, 3) # index of the top ocean cell + time = Time(clock.time) + @inbounds begin uₐ = surface_atmos_state.u[i, j, 1] vₐ = surface_atmos_state.v[i, j, 1] @@ -329,10 +330,9 @@ end inactive = inactive_node(i, j, kᴺ, grid, c, c, c) maxiter = ifelse(inactive, 1, similarity_theory.maxiter) - time = Time(clock.time) prescribed_heat_fluxes = net_downwelling_radiation(i, j, grid, time, radiation, Rs, Rℓ) - radiative_properties = local_radiation_properties(i, j, kᴺ, grid, radiation, time) + radiative_properties = local_radiation_properties(i, j, kᴺ, grid, time, radiation) turbulent_fluxes, surface_temperature = compute_similarity_theory_fluxes(similarity_theory, dynamic_ocean_state, diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl index 52176fd8..54706edc 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl @@ -91,9 +91,9 @@ function Base.show(io::IO, properties::SurfaceProperties) print(io, "└── sea_ice: ", summary(properties.sea_ice)) end -@inline local_radiation_properties(i, j, k, grid, ::Nothing, time) = nothing +@inline local_radiation_properties(i, j, k, grid, time, ::Nothing) = nothing -@inline function local_radiation_properties(i, j, k, grid, r::Radiation, time) +@inline function local_radiation_properties(i, j, k, grid, time, r::Radiation) σ = r.stefan_boltzmann_constant ϵ = stateindex(r.emission.ocean, i, j, k, grid, time) α = stateindex(r.reflection.ocean, i, j, k, grid, time) From 93133734b6f3ca710b8edd4957110acfaf2a20a0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:15:13 +0100 Subject: [PATCH 38/48] remove tyep instability --- .../CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index e0887f89..4d190cec 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -269,7 +269,8 @@ struct COARELogarithmicSimilarityProfile end # Initialize the solver iteration = ifelse(zero_shear_velocity, maxiter+1, 0) - Σ₀ = ifelse(zero_shear_velocity, SimilarityScales(0, 0, 0), Σ★) + Σₙ = SimilarityScales(zero(grid), zero(grid), zero(grid)) + Σ₀ = ifelse(zero_shear_velocity, Σₙ, Σ★) # Iterate until convergence while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) From 542efa1d3d0845bbcd54fb5c60522d762b02b183 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:18:07 +0100 Subject: [PATCH 39/48] better syntax --- .../similarity_theory_turbulent_fluxes.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 4d190cec..4713b2e7 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -252,8 +252,6 @@ struct COARELogarithmicSimilarityProfile end # Initial guess for the characteristic scales u★, θ★, q★. # Does not really matter if we are sophisticated or not, it converges # in about 10 iterations no matter what... - u★ = convert(FT, 1e-4) - Σ★ = SimilarityScales(u★, u★, u★) Δu, Δv = velocity_differences(atmos_state, surface_state, similarity_theory.bulk_velocity) # The inital velocity scale assumes that the gustiness velocity `Uᴳ` is equal to 0.5 ms⁻¹. @@ -264,13 +262,14 @@ struct COARELogarithmicSimilarityProfile end ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ²) # break the cycle if Δu == Δv == gustiness_parameter == 0 since that would mean - # that u★ == 0 so there is no turbulent transfer and the solver will not converge, leading to NaNs. + # that u★ == 0 so there is no turbulent transfer and the solver will not converge, leading to NaNs. zero_shear_velocity = (Δu == 0) && (Δv == 0) && (similarity_theory.gustiness_parameter == 0) # Initialize the solver iteration = ifelse(zero_shear_velocity, maxiter+1, 0) - Σₙ = SimilarityScales(zero(grid), zero(grid), zero(grid)) - Σ₀ = ifelse(zero_shear_velocity, Σₙ, Σ★) + u★ = ifelse(zero_shear_velocity, zero(grid), convert(FT, 1e-4)) + Σ★ = SimilarityScales(u★, u★, u★) + Σ₀ = Σ★ # Iterate until convergence while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) @@ -324,10 +323,10 @@ end # Iterating condition for the characteristic scales solvers @inline function iterating(Σ★, iteration, maxiter, solver) - havent_started = iteration == 0 + hasnt_started = iteration == 0 converged = norm(Σ★) < solver.tolerance reached_maxiter = iteration ≥ maxiter - return !(converged | reached_maxiter) | havent_started + return !(converged | reached_maxiter) | hasnt_started end """ From 0030b9d4590f89fc0d485f6830ec2e4696ff2de8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:20:48 +0100 Subject: [PATCH 40/48] && -> & --- .../CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 4713b2e7..21b25620 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -263,7 +263,7 @@ struct COARELogarithmicSimilarityProfile end # break the cycle if Δu == Δv == gustiness_parameter == 0 since that would mean # that u★ == 0 so there is no turbulent transfer and the solver will not converge, leading to NaNs. - zero_shear_velocity = (Δu == 0) && (Δv == 0) && (similarity_theory.gustiness_parameter == 0) + zero_shear_velocity = (Δu == 0) & (Δv == 0) & (similarity_theory.gustiness_parameter == 0) # Initialize the solver iteration = ifelse(zero_shear_velocity, maxiter+1, 0) From 60ad70458589a83207bf314440c3383746217903 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:26:44 +0100 Subject: [PATCH 41/48] another bugfix --- .../CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 21b25620..a2a952a9 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -267,7 +267,7 @@ struct COARELogarithmicSimilarityProfile end # Initialize the solver iteration = ifelse(zero_shear_velocity, maxiter+1, 0) - u★ = ifelse(zero_shear_velocity, zero(grid), convert(FT, 1e-4)) + u★ = ifelse(zero_shear_velocity, zero(FT), convert(FT, 1e-4)) Σ★ = SimilarityScales(u★, u★, u★) Σ₀ = Σ★ From c8db73971247dead2a9cd211eda6488eee64f0b9 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:33:19 +0100 Subject: [PATCH 42/48] do not change OceanSimulations.jl here --- src/OceanSimulations/OceanSimulations.jl | 29 ++++-------------------- 1 file changed, 4 insertions(+), 25 deletions(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 41c8b2b2..8c95965b 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -16,7 +16,6 @@ using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEEquation using SeawaterPolynomials.TEOS10: TEOS10EquationOfState -using Statistics: mean using Oceananigans.BuoyancyModels: g_Earth using Oceananigans.Coriolis: Ω_Earth @@ -43,30 +42,10 @@ default_or_override(override, alternative_default=nothing) = override # Some defaults default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) -function infer_maximum_Δt(grid) - Δx = mean(xspacings(grid)) - Δy = mean(yspacings(grid)) - Δθ = rad2deg(mean([Δx, Δy])) / grid.radius - - # The maximum Δt is roughly 40minutes / Δθ, giving: - # - 40 minutes for a 1 degree ocean - # - 20 minutes for a 1/4 degree ocean - # - 10 minutes for a 1/8 degree ocean - # - 5 minutes for a 1/16 degree ocean - # - 2.5 minutes for a 1/32 degree ocean - - return 40minutes / Δθ -end - +# 70 substeps is a safe rule of thumb for an ocean at 1/4 - 1/10th of a degree +# TODO: pass the cfl and a given Δt to calculate the number of substeps? const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} - -function default_free_surface(grid::TripolarOfSomeKind; - fixed_Δt = compute_maximum_Δt(grid), - cfl = 0.7) - free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) - @info "Using a $(free_surface)" - return free_surface -end +default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70) function default_ocean_closure() mixing_length = CATKEMixingLength(Cᵇ=0.01) @@ -96,7 +75,7 @@ default_tracer_advection() = FluxFormAdvection(WENO(order=7), # TODO: Specify the grid to a grid on the sphere; otherwise we can provide a different # function that requires latitude and longitude etc for computing coriolis=FPlane... function ocean_simulation(grid; - Δt = infer_maximum_Δt(grid), + Δt = 5minutes, closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), From 8616febd6a3967d15fc5dce0dd5db028800f3577 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:35:52 +0100 Subject: [PATCH 43/48] added a test for SkinTemperature zero fluxes --- test/test_surface_fluxes.jl | 59 +++++++++++++++++++------------------ 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index 32ce6afb..17e789ce 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -62,31 +62,36 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, water_mole_fraction = 1 # turbulent fluxes that force a specific humidity at the ocean's surface - similarity_theory = SimilarityTheoryTurbulentFluxes(grid; water_vapor_saturation, water_mole_fraction) - - # Thermodynamic parameters of the atmosphere - g = similarity_theory.gravitational_acceleration - 𝒬ₐ = Thermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₐ, qₐ) - cp = Thermodynamics.cp_m(ℂₐ, 𝒬ₐ) - ρₐ = Thermodynamics.air_density(ℂₐ, 𝒬ₐ) - ℰv = Thermodynamics.latent_heat_vapor(ℂₐ, 𝒬ₐ) - - # Ensure that the ΔT between atmosphere and ocean is zero - # Note that the Δθ accounts for the "lapse rate" at height h - Tₒ = Tₐ - celsius_to_kelvin + h / cp * g - - set!(ocean.model, u = uₐ, v = vₐ, T = Tₒ) - - # Compute the turbulent fluxes (neglecting radiation) - coupled_model = OceanSeaIceModel(ocean; atmosphere, similarity_theory) - turbulent_fluxes = coupled_model.fluxes.turbulent.fields - - # Make sure all fluxes are (almost) zero! - @test turbulent_fluxes.x_momentum[1, 1, 1] < eps(eltype(grid)) - @test turbulent_fluxes.y_momentum[1, 1, 1] < eps(eltype(grid)) - @test turbulent_fluxes.sensible_heat[1, 1, 1] < eps(eltype(grid)) - @test turbulent_fluxes.latent_heat[1, 1, 1] < eps(eltype(grid)) - @test turbulent_fluxes.water_vapor[1, 1, 1] < eps(eltype(grid)) + for Tmode in (BulkTemperature, SkinTemperature) + similarity_theory = SimilarityTheoryTurbulentFluxes(grid; + water_vapor_saturation, + water_mole_fraction, + surface_temperature_type = Tmode()) + + # Thermodynamic parameters of the atmosphere + g = similarity_theory.gravitational_acceleration + 𝒬ₐ = Thermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₐ, qₐ) + cp = Thermodynamics.cp_m(ℂₐ, 𝒬ₐ) + ρₐ = Thermodynamics.air_density(ℂₐ, 𝒬ₐ) + ℰv = Thermodynamics.latent_heat_vapor(ℂₐ, 𝒬ₐ) + + # Ensure that the ΔT between atmosphere and ocean is zero + # Note that the Δθ accounts for the "lapse rate" at height h + Tₒ = Tₐ - celsius_to_kelvin + h / cp * g + + set!(ocean.model, u = uₐ, v = vₐ, T = Tₒ) + + # Compute the turbulent fluxes (neglecting radiation) + coupled_model = OceanSeaIceModel(ocean; atmosphere, similarity_theory) + turbulent_fluxes = coupled_model.fluxes.turbulent.fields + + # Make sure all fluxes are (almost) zero! + @test turbulent_fluxes.x_momentum[1, 1, 1] < eps(eltype(grid)) + @test turbulent_fluxes.y_momentum[1, 1, 1] < eps(eltype(grid)) + @test turbulent_fluxes.sensible_heat[1, 1, 1] < eps(eltype(grid)) + @test turbulent_fluxes.latent_heat[1, 1, 1] < eps(eltype(grid)) + @test turbulent_fluxes.water_vapor[1, 1, 1] < eps(eltype(grid)) + end @info " Testing neutral fluxes..." @@ -200,10 +205,6 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, end end -@testset "SkinTemperature" begin - -end - @testset "Fluxes regression" begin for arch in test_architectures @info "Testing fluxes regression..." From 54dd6c1f08f1fb47c6faf43fe8f9c19780915422 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:36:22 +0100 Subject: [PATCH 44/48] import types --- test/test_surface_fluxes.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index 17e789ce..8babcccb 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -5,7 +5,9 @@ using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: convert_to_kelvin, SimilarityScales, seawater_saturation_specific_humidity, - surface_flux + surface_flux, + SkinTemperature, + BulkTemperature using Thermodynamics using CUDA From fc8ddc0d9be44e54633b44509e19d18c13426fdb Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 15:41:28 +0100 Subject: [PATCH 45/48] fix tests --- test/test_surface_fluxes.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index 8babcccb..6aba08c7 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -30,7 +30,6 @@ import Oceananigans.Fields: _fractional_indices _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, nothing, nothing) @testset "Test surface fluxes" begin - @info " Testing zero fluxes..." for arch in test_architectures grid = LatitudeLongitudeGrid(arch; size = 1, @@ -63,19 +62,22 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, water_vapor_saturation = FixedSpecificHumidity(qₐ) water_mole_fraction = 1 + # Thermodynamic parameters of the atmosphere + 𝒬ₐ = Thermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₐ, qₐ) + cp = Thermodynamics.cp_m(ℂₐ, 𝒬ₐ) + ρₐ = Thermodynamics.air_density(ℂₐ, 𝒬ₐ) + ℰv = Thermodynamics.latent_heat_vapor(ℂₐ, 𝒬ₐ) + # turbulent fluxes that force a specific humidity at the ocean's surface for Tmode in (BulkTemperature, SkinTemperature) + @info " Testing zero fluxes with $(Tmode)..." + similarity_theory = SimilarityTheoryTurbulentFluxes(grid; water_vapor_saturation, water_mole_fraction, surface_temperature_type = Tmode()) - # Thermodynamic parameters of the atmosphere g = similarity_theory.gravitational_acceleration - 𝒬ₐ = Thermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₐ, qₐ) - cp = Thermodynamics.cp_m(ℂₐ, 𝒬ₐ) - ρₐ = Thermodynamics.air_density(ℂₐ, 𝒬ₐ) - ℰv = Thermodynamics.latent_heat_vapor(ℂₐ, 𝒬ₐ) # Ensure that the ΔT between atmosphere and ocean is zero # Note that the Δθ accounts for the "lapse rate" at height h @@ -130,6 +132,7 @@ _fractional_indices(at_node, grid, ::Nothing, ::Nothing, ::Nothing) = (nothing, 𝒬ₒ = Thermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₒ, qₒ) qₒ = Thermodynamics.vapor_specific_humidity(ℂₐ, 𝒬ₒ) + g = similarity_theory.gravitational_acceleration # Differences! Δu = uₐ From 2ea71e349afda36e6c73a241eecc1100c577e55b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 18:26:36 +0100 Subject: [PATCH 46/48] Update src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl Co-authored-by: Gregory L. Wagner --- src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 91ac8997..10e52a5d 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -15,7 +15,7 @@ import Thermodynamics as AtmosphericThermodynamics """ struct BulkTemperature end -A type to represent the surface temperature used in the flux calculation. +Represents the surface temperature used in fixed-point iteration for surface fluxes following similarity theory. The surface temperature is not calculated but provided by either the ocean or the sea ice model. """ struct BulkTemperature end From f915ed34abbafb255ff1efd82324f96b6b489aaa Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 5 Dec 2024 15:41:55 +0100 Subject: [PATCH 47/48] change comment + only one papa simulation --- examples/single_column_os_papa_simulation.jl | 3 + .../os_papa_surface_temperature.jl | 298 ------------------ src/ClimaOcean.jl | 2 + .../CrossRealmFluxes/CrossRealmFluxes.jl | 4 +- .../CrossRealmFluxes/surface_temperature.jl | 2 +- src/OceanSeaIceModels/OceanSeaIceModels.jl | 1 + 6 files changed, 10 insertions(+), 300 deletions(-) delete mode 100644 experiments/single_column_experiments/os_papa_surface_temperature.jl diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl index 2f12a0d1..11e5b585 100644 --- a/examples/single_column_os_papa_simulation.jl +++ b/examples/single_column_os_papa_simulation.jl @@ -97,8 +97,11 @@ lines!(axq, t_days, qa) current_figure() # We continue constructing a simulation. +# For the fluxes computation we use a `SkinTemperature` formulation that computes +# the skin temperature from a balance between internal and external heat fluxes. radiation = Radiation() +similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type=SkinTemperature()) coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) diff --git a/experiments/single_column_experiments/os_papa_surface_temperature.jl b/experiments/single_column_experiments/os_papa_surface_temperature.jl deleted file mode 100644 index 2bed91d2..00000000 --- a/experiments/single_column_experiments/os_papa_surface_temperature.jl +++ /dev/null @@ -1,298 +0,0 @@ -# This simulation tests the difference between using a -# `BulkTemperature` and a `SkinTemperature` - -using ClimaOcean -using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: SimilarityTheoryTurbulentFluxes, SkinTemperature -using Oceananigans -using Oceananigans.Units -using Oceananigans.BuoyancyModels: buoyancy_frequency -using Oceananigans.Units: Time -using Printf - -# Ocean station papa location -location_name = "ocean_station_papa" -λ★, φ★ = 35.1, 50.1 - -grid = RectilinearGrid(size = 200, - x = λ★, - y = φ★, - z = (-400, 0), - topology = (Flat, Flat, Bounded)) - -ocean1 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) -ocean2 = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) - -# We set initial conditions from ECCO: -set!(ocean1.model.tracers.T, ECCOMetadata(:temperature)) -set!(ocean1.model.tracers.S, ECCOMetadata(:salinity)) -set!(ocean2.model.tracers.T, ECCOMetadata(:temperature)) -set!(ocean2.model.tracers.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()) - -# We continue constructing a simulation. - -radiation = Radiation() -coupled_model_prescribed = OceanSeaIceModel(ocean1; atmosphere, radiation) - -similarity_theory = SimilarityTheoryTurbulentFluxes(grid; surface_temperature_type = SkinTemperature(; κ = 0.5)) -coupled_model_diagnostic = OceanSeaIceModel(ocean2; atmosphere, radiation, similarity_theory) - -for (coupled_model, suffix) in zip([coupled_model_diagnostic, coupled_model_prescribed], - ["diagnostic", "prescribed"]) - - simulation = Simulation(coupled_model, Δt=ocean1.Δt, stop_time=10days) - - wall_clock = Ref(time_ns()) - - function progress(sim) - msg = "Ocean Station Papa" - msg *= string(", iter: ", iteration(sim), ", time: ", prettytime(sim)) - - elapsed = 1e-9 * (time_ns() - wall_clock[]) - msg *= string(", wall time: ", prettytime(elapsed)) - wall_clock[] = time_ns() - - u, v, w = sim.model.ocean.model.velocities - msg *= @sprintf(", max|u|: (%.2e, %.2e)", maximum(abs, u), maximum(abs, v)) - - T = sim.model.ocean.model.tracers.T - S = sim.model.ocean.model.tracers.S - e = sim.model.ocean.model.tracers.e - - τx = first(sim.model.fluxes.total.ocean.momentum.u) - τy = first(sim.model.fluxes.total.ocean.momentum.v) - Q = first(sim.model.fluxes.total.ocean.heat) - - u★ = sqrt(sqrt(τx^2 + τy^2)) - Ts = first(interior(sim.model.fluxes.turbulent.fields.T_surface, 1, 1, 1)) - - Nz = size(T, 3) - msg *= @sprintf(", u★: %.2f m s⁻¹", u★) - msg *= @sprintf(", Q: %.2f W m⁻²", Q) - msg *= @sprintf(", T₀: %.2f ᵒC", Ts) - msg *= @sprintf(", S₀: %.2f g/kg", first(interior(S, 1, 1, Nz))) - msg *= @sprintf(", e₀: %.2e m² s⁻²", first(interior(e, 1, 1, Nz))) - - @info msg - end - - simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) - - # Build flux outputs - τx = coupled_model.fluxes.total.ocean.momentum.u - τy = coupled_model.fluxes.total.ocean.momentum.v - JT = coupled_model.fluxes.total.ocean.tracers.T - Js = coupled_model.fluxes.total.ocean.tracers.S - E = coupled_model.fluxes.turbulent.fields.water_vapor - Qc = coupled_model.fluxes.turbulent.fields.sensible_heat - Qv = coupled_model.fluxes.turbulent.fields.latent_heat - Ts = coupled_model.fluxes.turbulent.fields.T_surface - ρₒ = coupled_model.fluxes.ocean_reference_density - cₚ = coupled_model.fluxes.ocean_heat_capacity - - Q = ρₒ * cₚ * JT - ρτx = ρₒ * τx - ρτy = ρₒ * τy - N² = buoyancy_frequency(coupled_model.ocean.model) - κc = coupled_model.ocean.model.diffusivity_fields.κc - - fluxes = (; ρτx, ρτy, E, Js, Qv, Qc, Ts) - auxiliary_fields = (; N², κc) - fields = merge(coupled_model.ocean.model.velocities, coupled_model.ocean.model.tracers, auxiliary_fields) - - # Slice fields at the surface - outputs = merge(fields, fluxes) - - filename = "single_column_omip_$(location_name)_$(suffix)" - - simulation.output_writers[:jld2] = JLD2OutputWriter(coupled_model.ocean.model, outputs; filename, - schedule = TimeInterval(3hours), - overwrite_existing = true) - - run!(simulation) -end - -##### -##### Visualization -##### - -filename_prescribed = "single_column_omip_$(location_name)_prescribed.jld2" -filename_diagnostic = "single_column_omip_$(location_name)_diagnostic.jld2" - -# Diagnosed -ud = FieldTimeSeries(filename_diagnostic, "u") -vd = FieldTimeSeries(filename_diagnostic, "v") -Td = FieldTimeSeries(filename_diagnostic, "v") -ed = FieldTimeSeries(filename_diagnostic, "T") -Sd = FieldTimeSeries(filename_diagnostic, "S") -Nd² = FieldTimeSeries(filename_diagnostic, "N²") -κd = FieldTimeSeries(filename_diagnostic, "κc") - -Tsd = FieldTimeSeries(filename_diagnostic, "Ts") -Qvd = FieldTimeSeries(filename_diagnostic, "Qv") -Qcd = FieldTimeSeries(filename_diagnostic, "Qc") -Jsd = FieldTimeSeries(filename_diagnostic, "Js") -Evd = FieldTimeSeries(filename_diagnostic, "E") -ρτxd = FieldTimeSeries(filename_diagnostic, "ρτx") -ρτyd = FieldTimeSeries(filename_diagnostic, "ρτy") - -# Prescribed -up = FieldTimeSeries(filename_prescribed, "u") -vp = FieldTimeSeries(filename_prescribed, "v") -Tp = FieldTimeSeries(filename_prescribed, "T") -Sp = FieldTimeSeries(filename_prescribed, "S") -ep = FieldTimeSeries(filename_prescribed, "e") -Np² = FieldTimeSeries(filename_prescribed, "N²") -κp = FieldTimeSeries(filename_prescribed, "κc") - -Tsp = FieldTimeSeries(filename_prescribed, "Ts") -Qvp = FieldTimeSeries(filename_prescribed, "Qv") -Qcp = FieldTimeSeries(filename_prescribed, "Qc") -Jsp = FieldTimeSeries(filename_prescribed, "Js") -Evp = FieldTimeSeries(filename_prescribed, "E") -ρτxp = FieldTimeSeries(filename_prescribed, "ρτx") -ρτyp = FieldTimeSeries(filename_prescribed, "ρτy") - -Nz = size(ud, 3) -times = Qcd.times - -fig = Figure(size=(1800, 1800)) - -axτ = Axis(fig[1, 1:3], xlabel="Days since Oct 1 1992", ylabel="Wind stress (N m⁻²)") -axQ = Axis(fig[1, 4:6], xlabel="Days since Oct 1 1992", ylabel="Heat flux (W m⁻²)") -axu = Axis(fig[2, 1:3], xlabel="Days since Oct 1 1992", ylabel="Velocities (m s⁻¹)") -axT = Axis(fig[2, 4:6], xlabel="Days since Oct 1 1992", ylabel="Surface temperature (ᵒC)") -axF = Axis(fig[3, 1:3], xlabel="Days since Oct 1 1992", ylabel="Freshwater volume flux (m s⁻¹)") -axS = Axis(fig[3, 4:6], xlabel="Days since Oct 1 1992", ylabel="Surface salinity (g kg⁻¹)") - -axuz = Axis(fig[4:5, 1:2], xlabel="Velocities (m s⁻¹)", ylabel="z (m)") -axTz = Axis(fig[4:5, 3:4], xlabel="Temperature (ᵒC)", ylabel="z (m)") -axSz = Axis(fig[4:5, 5:6], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)") -axNz = Axis(fig[6:7, 1:2], xlabel="Buoyancy frequency (s⁻²)", ylabel="z (m)") -axκz = Axis(fig[6:7, 3:4], xlabel="Eddy diffusivity (m² s⁻¹)", ylabel="z (m)", xscale=log10) -axez = Axis(fig[6:7, 5:6], xlabel="Turbulent kinetic energy (m² s⁻²)", ylabel="z (m)", xscale=log10) - -title = @sprintf("Single-column simulation at %.2f, %.2f", φ★, λ★) -Label(fig[0, 1:6], title) - -n = Observable(1) - -times = (times .- times[1]) ./days -Nt = length(times) -tn = @lift times[$n] - -colors = Makie.wong_colors() - -ρₒ = coupled_model_prescribed.fluxes.ocean_reference_density -τxp = interior(ρτxp, 1, 1, 1, :) ./ ρₒ -τyp = interior(ρτyp, 1, 1, 1, :) ./ ρₒ -up★ = @. (τxp^2 + τyp^2)^(1/4) -τxd = interior(ρτxd, 1, 1, 1, :) ./ ρₒ -τyd = interior(ρτyd, 1, 1, 1, :) ./ ρₒ -ud★ = @. (τxd^2 + τyd^2)^(1/4) - -lines!(axu, times, interior(ud, 1, 1, Nz, :), color=colors[1], label="Zonal") -lines!(axu, times, interior(vd, 1, 1, Nz, :), color=colors[2], label="Meridional") -lines!(axu, times, interior(up, 1, 1, Nz, :), color=colors[1], linestyle = :dash, label="Zonal") -lines!(axu, times, interior(vp, 1, 1, Nz, :), color=colors[2], linestyle = :dash, label="Meridional") -lines!(axu, times, ud★, color=colors[3], label="Ocean-side u★") -lines!(axu, times, up★, color=colors[3], label="Ocean-side u★", linestyle = :dash) -vlines!(axu, tn, linewidth=4, color=(:black, 0.5)) -axislegend(axu) - -lines!(axτ, times, interior(ρτxd, 1, 1, 1, :), label="Zonal") -lines!(axτ, times, interior(ρτyd, 1, 1, 1, :), label="Meridional") -lines!(axτ, times, interior(ρτxp, 1, 1, 1, :), linestyle=:dash, label="Zonal") -lines!(axτ, times, interior(ρτyp, 1, 1, 1, :), linestyle=:dash, label="Meridional") -vlines!(axτ, tn, linewidth=4, color=(:black, 0.5)) -axislegend(axτ) - -lines!(axT, times, interior(Tsd, 1, 1, 1, :), color=colors[2], linewidth=4, label="Ocean surface temperature") -lines!(axT, times, interior(Tsp, 1, 1, 1, :), color=colors[2], linewidth=4, linestyle = :dash, label="Ocean surface temperature") -vlines!(axT, tn, linewidth=4, color=(:black, 0.5)) -axislegend(axT) - -lines!(axQ, times, interior(Qvd, 1, 1, 1, 1:Nt), color=colors[2], label="Sensible", linewidth=2) -lines!(axQ, times, interior(Qcd, 1, 1, 1, 1:Nt), color=colors[3], label="Latent", linewidth=2) -lines!(axQ, times, interior(Qvp, 1, 1, 1, 1:Nt), color=colors[2], linestyle=:dash, label="Sensible", linewidth=2) -lines!(axQ, times, interior(Qcp, 1, 1, 1, 1:Nt), color=colors[3], linestyle=:dash, label="Latent", linewidth=2) -vlines!(axQ, tn, linewidth=4, color=(:black, 0.5)) -axislegend(axQ) - -lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), label="Evaporation") -lines!(axF, times, - interior(Evd, 1, 1, 1, 1:Nt), linestyle=:dash, label="Evaporation") -vlines!(axF, tn, linewidth=4, color=(:black, 0.5)) -axislegend(axF) - -lines!(axS, times, interior(Sd, 1, 1, Nz, :)) -lines!(axS, times, interior(Sp, 1, 1, Nz, :), linestyle=:dash) -vlines!(axS, tn, linewidth=4, color=(:black, 0.5)) - -zc = znodes(Tp) -zf = znodes(κp) -upn = @lift interior(up[$n], 1, 1, :) -vpn = @lift interior(vp[$n], 1, 1, :) -Tpn = @lift interior(Tp[$n], 1, 1, :) -Spn = @lift interior(Sp[$n], 1, 1, :) -κpn = @lift interior(κp[$n], 1, 1, :) -epn = @lift interior(ep[$n], 1, 1, :) -Np²n = @lift interior(Np²[$n], 1, 1, :) - -udn = @lift interior(ud[$n], 1, 1, :) -vdn = @lift interior(vd[$n], 1, 1, :) -Tdn = @lift interior(Td[$n], 1, 1, :) -Sdn = @lift interior(Sd[$n], 1, 1, :) -κdn = @lift interior(κd[$n], 1, 1, :) -edn = @lift interior(ed[$n], 1, 1, :) -Nd²n = @lift interior(Nd²[$n], 1, 1, :) - -scatterlines!(axuz, udn, zc, label="u") -scatterlines!(axuz, vdn, zc, label="v") -scatterlines!(axTz, Tdn, zc) -scatterlines!(axSz, Sdn, zc) -scatterlines!(axez, edn, zc) -scatterlines!(axNz, Nd²n, zf) -scatterlines!(axκz, κdn, zf) -scatterlines!(axuz, upn, zc, linestyle=:dash, label="u") -scatterlines!(axuz, vpn, zc, linestyle=:dash, label="v") -scatterlines!(axTz, Tpn, zc, linestyle=:dash) -scatterlines!(axSz, Spn, zc, linestyle=:dash) -scatterlines!(axez, epn, zc, linestyle=:dash) -scatterlines!(axNz, Np²n, zf, linestyle=:dash) -scatterlines!(axκz, κpn, zf, linestyle=:dash) - -axislegend(axuz) - -ulim = max(maximum(abs, up), maximum(abs, vp)) -xlims!(axuz, -ulim, ulim) - -Tmax = maximum(interior(Tp)) -Tmin = minimum(interior(Tp)) -xlims!(axTz, Tmin - 0.1, Tmax + 0.1) - -Nmax = maximum(interior(Np²)) -xlims!(axNz, -Nmax/10, Nmax * 1.05) - -κmax = maximum(interior(κp)) -xlims!(axκz, 1e-9, κmax * 1.1) - -emax = maximum(interior(ep)) -xlims!(axez, 1e-11, emax * 1.1) - -Smax = maximum(interior(Sp)) -Smin = minimum(interior(Sp)) -xlims!(axSz, Smin - 0.2, Smax + 0.2) - -record(fig, "single_column_profiles.mp4", 1:Nt, framerate=24) do nn - @info "Drawing frame $nn of $Nt..." - n[] = nn -end -nothing #hide -# ![](single_column_profiles.mp4) diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index ae1bb084..b2c4f261 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -13,6 +13,8 @@ export Radiation, LatitudeDependentAlbedo, SimilarityTheoryTurbulentFluxes, + SkinTemperature, + BulkTemperature, JRA55_prescribed_atmosphere, JRA55NetCDFBackend, regrid_bathymetry, diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl index 9cac3060..e8bb9b4c 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl @@ -6,7 +6,9 @@ using Adapt export Radiation, OceanSeaIceSurfaceFluxes, LatitudeDependentAlbedo, - SimilarityTheoryTurbulentFluxes + SimilarityTheoryTurbulentFluxes, + SkinTemperature, + BulkTemperature using ..OceanSeaIceModels: default_gravitational_acceleration diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 10e52a5d..d732da2c 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -73,7 +73,7 @@ DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert( # # Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) # plus the (semi-implicit) outgoing longwave radiation and the RHS are the remaining atmospheric and radiative fluxes -# provided explicitly. +# provided explicitly. Here we implement the fully explicit version. @inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = θₒ - Jᵀ / F.κ * F.δ # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 26e16348..85e69d94 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -2,6 +2,7 @@ module OceanSeaIceModels export OceanSeaIceModel, SimilarityTheoryTurbulentFluxes, FreezingLimitedOceanTemperature export Radiation, LatitudeDependentAlbedo +export SkinTemperature, BulkTemperature using Oceananigans using SeawaterPolynomials From 08be1bbb1ea68d64c318f3ef5a56043a6b8868ed Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 5 Dec 2024 15:55:03 +0100 Subject: [PATCH 48/48] add comment --- src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index d732da2c..3da5b500 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -73,7 +73,8 @@ DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert( # # Where the LHS is the internal diffusive flux inside the ocean (within the boundary layer of thickness δ) # plus the (semi-implicit) outgoing longwave radiation and the RHS are the remaining atmospheric and radiative fluxes -# provided explicitly. Here we implement the fully explicit version. +# provided explicitly. Here we implement the fully explicit version, the linearized version is an optimization +# that can be explored in the future. @inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = θₒ - Jᵀ / F.κ * F.δ # he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the