Merge branch 'main' into ss-nc/fix-climatoloy-experiment

Author: simone-silvestri

src/ClimaOcean.jl

@@ -53,21 +53,20 @@ const SKOFTS = SomeKindOfFieldTimeSeries
 @inline stateindex(a::AbstractArray, i, j, k, args...) = @inbounds a[i, j, k]
 @inline stateindex(a::SKOFTS, i, j, k, grid, time, args...) = @inbounds a[i, j, k, time]
 
-@inline function stateindex(a::Function, i, j, k, grid, time, loc)
-    LX, LY, LZ = loc
+@inline function stateindex(a::Function, i, j, k, grid, time, (LX, LY, LZ), args...)
     λ, φ, z = node(i, j, k, grid, LX(), LY(), LZ())
     return a(λ, φ, z, time)
 end
 
-@inline function stateindex(a::Tuple, i, j, k, grid, time)
+@inline function stateindex(a::Tuple, i, j, k, grid, time, args...)
     N = length(a)
     ntuple(Val(N)) do n
-        stateindex(a[n], i, j, k, grid, time)
+        stateindex(a[n], i, j, k, grid, time, args...)
     end
 end
 
-@inline function stateindex(a::NamedTuple, i, j, k, grid, time)
-    vals = stateindex(values(a), i, j, k, grid, time)
+@inline function stateindex(a::NamedTuple, i, j, k, grid, time, args...)
+    vals = stateindex(values(a), i, j, k, grid, time, args...)
     names = keys(a)
     return NamedTuple{names}(vals)
 end

src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl

@@ -100,14 +100,15 @@ end
 
     # Compute radiation fluxes
     σ = atmos_ocean_properties.radiation.σ
-    α = stateindex(atmos_ocean_properties.radiation.α, i, j, kᴺ, grid, time)
-    ϵ = stateindex(atmos_ocean_properties.radiation.ϵ, i, j, kᴺ, grid, time)
-    Qu = upwelling_radiation(Tₛ, σ, ϵ)
-    Qdℓ = downwelling_longwave_radiation(Qℓ, ϵ)
-    Qds = downwelling_shortwave_radiation(Qs, α)
+    α = atmos_ocean_properties.radiation.α
+    ϵ = atmos_ocean_properties.radiation.ϵ
+    Qu  = upwelling_radiation(i, j, kᴺ, grid, time, Tₛ, σ, ϵ) 
+    Qdℓ = downwelling_longwave_radiation(i, j, kᴺ, grid, time, ϵ, Qℓ)
+    Qds = downwelling_shortwave_radiation(i, j, kᴺ, grid, time, α, Qs)
     ΣQao = Qu + Qc + Qv + Qdℓ + Qds
 
     @inbounds begin
+        # Write radiative components of the heat flux for diagnostic purposes
         atmos_ocean_fluxes.upwelling_longwave[i, j, 1] = Qu
         atmos_ocean_fluxes.downwelling_longwave[i, j, 1] = Qdℓ
         atmos_ocean_fluxes.downwelling_shortwave[i, j, 1] = Qds
@@ -244,10 +245,10 @@ end
 
     # Compute radiation fluxes
     σ = atmos_sea_ice_properties.radiation.σ
-    α = stateindex(atmos_sea_ice_properties.radiation.α, i, j, kᴺ, grid, time)
-    ϵ = stateindex(atmos_sea_ice_properties.radiation.ϵ, i, j, kᴺ, grid, time)
-    Qu = upwelling_radiation(Ts, σ, ϵ)
-    Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ)
+    α = atmos_sea_ice_properties.radiation.α
+    ϵ = atmos_sea_ice_properties.radiation.ϵ
+    Qu = upwelling_radiation(i, j, kᴺ, grid, time, Ts, σ, ϵ) 
+    Qd = net_downwelling_radiation(i, j, kᴺ, grid, time, α, ϵ, Qs, Qℓ)
 
     ΣQt = (Qd + Qu + Qc + Qv) * ℵi # If ℵi == 0 there is no heat flux from the top!
     ΣQb = Qf + Qi

src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl

@@ -122,6 +122,16 @@ end
     needs_to_converge = stop_criteria isa ConvergenceStopCriteria
     not_water = inactive_node(i, j, kᴺ, grid, Center(), Center(), Center())
 
+    # Compute local radiative properties and rebuild the interface properties
+    α = stateindex(interface_properties.radiation.α, i, j, kᴺ, grid, time, (Center, Center, Center), Qs)
+    ϵ = stateindex(interface_properties.radiation.ϵ, i, j, kᴺ, grid, time, (Center, Center, Center))
+    σ = interface_properties.radiation.σ
+
+    interface_properties = InterfaceProperties((; α, ϵ, σ),
+                                               interface_properties.specific_humidity_formulation,
+                                               interface_properties.temperature_formulation,
+                                               interface_properties.velocity_formulation)
+
     if needs_to_converge && not_water
         interface_state = zero_interface_state(FT)
     else

src/OceanSeaIceModels/InterfaceComputations/interface_states.jl

@@ -299,8 +299,10 @@ end
     ϵ = interface_properties.radiation.ϵ
     α = interface_properties.radiation.α
 
+    Qs = downwelling_radiation.Qs
+    Qℓ = downwelling_radiation.Qℓ
     Qu = upwelling_radiation(Tₛ⁻, σ, ϵ)
-    Qd = net_downwelling_radiation(downwelling_radiation, α, ϵ)
+    Qd = net_downwelling_radiation(Qs, Qℓ, α, ϵ)
 
     u★ = interface_state.u★
     θ★ = interface_state.θ★

src/OceanSeaIceModels/InterfaceComputations/latitude_dependent_albedo.jl

@@ -1,4 +1,5 @@
 using Oceananigans.Grids: ηnode
+using Oceananigans.Fields: instantiate
 
 struct LatitudeDependentAlbedo{FT}
     direct :: FT
@@ -45,8 +46,8 @@ end
 
 Base.show(io::IO, α::LatitudeDependentAlbedo) = print(io, summary(α))
 
-@inline function stateindex(α::LatitudeDependentAlbedo, i, j, k, grid, time)
-    φ = ηnode(i, j, k, grid, Center(), Center(), Center())
+@inline function stateindex(α::LatitudeDependentAlbedo, i, j, k, grid, time, (LX, LY, LZ), args...)
+    φ = ηnode(i, j, k, grid, LX(), LY(), LZ())
     α₀ = α.diffuse
     α₁ = α.direct
     return α₀ - α₁ * hack_cosd(2φ)

src/OceanSeaIceModels/InterfaceComputations/radiation.jl

@@ -17,7 +17,7 @@ Adapt.adapt_structure(to, r :: Radiation) =  Radiation(Adapt.adapt(to, r.emissio
     Radiation([arch = CPU(), FT=Float64];
               ocean_emissivity = 0.97,
               sea_ice_emissivity = 1.0,
-              ocean_albedo = LatitudeDependentAlbedo(FT),
+              ocean_albedo = 0.05,
               sea_ice_albedo = 0.7,
               stefan_boltzmann_constant = 5.67e-8)
 
@@ -34,14 +34,14 @@ Keyword Arguments
 
 - `ocean_emissivity`: The emissivity of the ocean surface. Default: `0.97`.
 - `sea_ice_emissivity`: The emissivity of the sea ice surface. Default: `1.0`.
-- `ocean_albedo`: The albedo of the ocean surface. Default: `LatitudeDependentAlbedo(FT)`.
+- `ocean_albedo`: The albedo of the ocean surface. Default: `0.05`.
 - `sea_ice_albedo`: The albedo of the sea ice surface. Default: `0.7`.
 - `stefan_boltzmann_constant`: The Stefan-Boltzmann constant. Default: `5.67e-8`.
 """
 function Radiation(arch = CPU(), FT=Oceananigans.defaults.FloatType;
                    ocean_emissivity = 0.97,
                    sea_ice_emissivity = 1.0,
-                   ocean_albedo = LatitudeDependentAlbedo(FT),
+                   ocean_albedo = 0.05,
                    sea_ice_albedo = 0.7,
                    stefan_boltzmann_constant = 5.67e-8)
 
@@ -90,10 +90,30 @@ function Base.show(io::IO, properties::SurfaceProperties)
     print(io, "└── sea_ice: ", summary(properties.sea_ice))
 end
 
-@inline upwelling_radiation(T, σ, ϵ) = σ * ϵ * T^4
-@inline net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ) = - (1 - α) * Qs - ϵ * Qℓ
-@inline net_downwelling_radiation(r, α, ϵ) = - (1 - α) * r.Qs - ϵ * r.Qℓ
+const CCC = (Center, Center, Center)
+
+@inline function upwelling_radiation(i, j, k, grid, time, T, σ, ϵ) 
+    ϵi = stateindex(ϵ, i, j, k, grid, time, CCC)
+    return σ * ϵi * T^4
+end
 
 # Split the individual bands
-@inline downwelling_longwave_radiation(Qℓ, ϵ) = - ϵ * Qℓ
-@inline downwelling_shortwave_radiation(Qs, α) = - (1 - α) * Qs 
+@inline function downwelling_longwave_radiation(i, j, k, grid, time, ϵ, Qℓ)  
+    ϵi = stateindex(ϵ, i, j, k, grid, time, CCC)
+    return - ϵi * Qℓ
+end
+
+@inline function downwelling_shortwave_radiation(i, j, k, grid, time, α, Qs)  
+    αi = stateindex(α, i, j, k, grid, time, CCC, Qs) 
+    return - (1 - αi) * Qs 
+end
+
+@inline net_downwelling_radiation(i, j, k, grid, time, α, ϵ, Qs, Qℓ) = 
+    downwelling_shortwave_radiation(i, j, k, grid, time, α, Qs) + 
+    downwelling_longwave_radiation(i, j, k, grid, time, ϵ, Qℓ)
+
+# Inside the solver we lose both spatial and temporal information, but the
+# radiative properties have already been computed correctly
+@inline net_downwelling_radiation(Qs, Qℓ, α, ϵ) = - (1 - α) * Qs - ϵ * Qℓ
+
+@inline upwelling_radiation(T, σ, ϵ) = σ * ϵ * T^4

src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl

@@ -277,7 +277,7 @@ These stability functions are obtained by regression to experimental data.
 The stability parameter for stable atmospheric conditions is defined as
 ```math
 dζ = min(ζmax, Aˢζ)
-ψₛ = - (Bˢ ζ + Cˢ ( ζ - Dˢ ) ) exp( - dζ) - Cˢ Dˢ
+ψₛ = - Bˢ * ζ⁺ - Cˢ * (ζ⁺ - Dˢ) * exp(- dζ) - Cˢ * Dˢ
 ```
 
 While the stability parameter for unstable atmospheric conditions is calculated

src/OceanSeaIceModels/InterfaceComputations/tabulated_albedo.jl

@@ -115,10 +115,9 @@ Base.show(io::IO, α::TabulatedAlbedo) = print(io, summary(α))
 @inline simulation_day(time::Time{<:Number})      = time.time ÷ 86400
 @inline seconds_in_day(time::Time{<:Number}, day) = time.time - day * 86400
 
-@inline function net_downwelling_radiation(i, j, grid, time, radiation::Radiation{<:Any, <:Any, <:SurfaceProperties{<:TabulatedAlbedo}}, Qs, Qℓ)
-    α = radiation.reflection.ocean
+@inline function stateindex(α::TabulatedAlbedo, i, j, k, grid, time, loc, Qs)
     FT = eltype(α)
-    λ, φ, z = _node(i, j, 1, grid, Center(), Center(), Center())
+    λ, φ, z = _node(i, j, k, grid, Center(), Center(), Center())
 
     φ = deg2rad(φ)
     λ = deg2rad(λ)
@@ -164,8 +163,6 @@ Base.show(io::IO, α::TabulatedAlbedo) = print(io, summary(α))
                    ϕ₃(ξ, η) * getindex(α.α_table, i⁺, j⁻) +
                    ϕ₄(ξ, η) * getindex(α.α_table, i⁺, j⁺)
 
-    ϵ = stateindex(radiation.emission.ocean, i, j, 1, grid, time)
-
-    return - (1 - α) * Qs - ϵ * Qℓ
+    return α
 end