Compute reference density with derivative

src/Thermodynamics/reference_states.jl

@@ -1,4 +1,4 @@
-using Oceananigans: Oceananigans, Center, Field, set!, fill_halo_regions!
+using Oceananigans: Oceananigans, Center, Field, set!, fill_halo_regions!, ∂z
 using Oceananigans.BoundaryConditions: FieldBoundaryConditions, ValueBoundaryCondition
 
 using Adapt: Adapt, adapt
@@ -100,17 +100,18 @@ function ReferenceState(grid, constants=ThermodynamicConstants(eltype(grid));
     θ₀ = convert(FT, potential_temperature)
     pˢᵗ = convert(FT, standard_pressure)
     loc = (nothing, nothing, Center())
-
-    ρ₀ = surface_density(p₀, θ₀, constants)
-    ρ_bcs = FieldBoundaryConditions(grid, loc, bottom=ValueBoundaryCondition(ρ₀))
-    ρᵣ = Field{Nothing, Nothing, Center}(grid, boundary_conditions=ρ_bcs)
-    set!(ρᵣ, z -> adiabatic_hydrostatic_density(z, p₀, θ₀, constants))
-    fill_halo_regions!(ρᵣ)
-
+
     p_bcs = FieldBoundaryConditions(grid, loc, bottom=ValueBoundaryCondition(p₀))
     pᵣ = Field{Nothing, Nothing, Center}(grid, boundary_conditions=p_bcs)
     set!(pᵣ, z -> adiabatic_hydrostatic_pressure(z, p₀, θ₀, constants))
     fill_halo_regions!(pᵣ)
+
+    ρ₀ = surface_density(p₀, θ₀, constants)
+    ρ_bcs = FieldBoundaryConditions(grid, loc, bottom=ValueBoundaryCondition(ρ₀))
+    ρᵣ = Field{Nothing, Nothing, Center}(grid, boundary_conditions=ρ_bcs)
+    g = thermo.gravitational_acceleration
+    set!(ρᵣ, - ∂z(pᵣ) / g)
+    fill_halo_regions!(ρᵣ)
 
     return ReferenceState(p₀, θ₀, pˢᵗ, pᵣ, ρᵣ)
 end