A first stab at a diagnostic surface temperature implementation

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.3 - 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"

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.
 

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)
 

src/ClimaOcean.jl

@@ -13,6 +13,8 @@ export
     Radiation,
     LatitudeDependentAlbedo,
     SimilarityTheoryTurbulentFluxes,
+    SkinTemperature,
+    BulkTemperature,
     JRA55_prescribed_atmosphere,
     JRA55NetCDFBackend,
     regrid_bathymetry,

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)

src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl

@@ -6,7 +6,9 @@ using Adapt
 export Radiation,
        OceanSeaIceSurfaceFluxes,
        LatitudeDependentAlbedo,
-       SimilarityTheoryTurbulentFluxes
+       SimilarityTheoryTurbulentFluxes,
+       SkinTemperature, 
+       BulkTemperature
 
 using ..OceanSeaIceModels: default_gravitational_acceleration
 
@@ -31,6 +33,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")

src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl

@@ -141,8 +141,11 @@ 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,
             atmosphere.reference_height, # height at which the state is known
             atmosphere.boundary_layer_height,
             atmosphere.thermodynamics_parameters)   
@@ -264,21 +267,27 @@ end
                                                                      grid,
                                                                      clock,
                                                                      ocean_state,
+                                                                     ocean_density,
+                                                                     ocean_heat_capacity,
                                                                      ocean_temperature_units,
                                                                      surface_atmos_state,
+                                                                     radiation,
                                                                      atmosphere_reference_height,
                                                                      atmosphere_boundary_layer_height,
                                                                      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]
         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,18 +331,27 @@ 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ℓ) 
+    radiative_properties = local_radiation_properties(i, j, kᴺ, grid, time, radiation)
+
+    turbulent_fluxes, surface_temperature = compute_similarity_theory_fluxes(similarity_theory,
+                                                                             dynamic_ocean_state, 
+                                                                             dynamic_atmos_state, 
+                                                                             prescribed_heat_fluxes,
+                                                                             radiative_properties,
+                                                                             Sₒ,
+                                                                             ocean_density,
+                                                                             ocean_heat_capacity,
+                                                                             atmosphere_boundary_layer_height,
+                                                                             ℂₐ, g, ϰ, maxiter)
 
     # Store fluxes
     Qv = similarity_theory.fields.latent_heat
     Qc = similarity_theory.fields.sensible_heat
     Fv = similarity_theory.fields.water_vapor
     ρτx = similarity_theory.fields.x_momentum
     ρτy = similarity_theory.fields.y_momentum
+    Ts  = similarity_theory.fields.T_surface
 
     @inbounds begin
         # +0: cooling, -0: heating
@@ -342,6 +360,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
 

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, time, ::Nothing) = nothing
+
+@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)
+
+    return (; ϵ, α, σ)
+end