explicit surface flux

CliMA · Mar 2, 2025 · 84537d0 · 84537d0
1 parent 48b6e9f
commit 84537d0
Show file tree

Hide file tree

Showing 15 changed files with 104 additions and 44 deletions.
diff --git a/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean.yml b/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean.yml
@@ -11,4 +11,4 @@ rad: "clearsky"
 insolation: "timevarying"
 dt_rad: "1hours"
 prognostic_surface: "PrognosticSurfaceTemperature"
-#check_conservation: true
+check_conservation: true
diff --git a/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean_ft64.yml b/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean_ft64.yml
@@ -11,4 +11,4 @@ rad: "clearsky"
 insolation: "timevarying"
 dt_rad: "1hours"
 prognostic_surface: "PrognosticSurfaceTemperature"
-#check_conservation: true
+check_conservation: true
diff --git a/config/model_configs/baroclinic_wave_conservation.yml b/config/model_configs/baroclinic_wave_conservation.yml
@@ -2,6 +2,7 @@ initial_condition: "DryBaroclinicWave"
 dt: "400secs"
 t_end: "5days"
 dt_save_state_to_disk: "5days"
+surface_setup: "NoSurface"
 check_conservation: true
 diagnostics:
   - short_name: [massa, energya]

diff --git a/config/model_configs/baroclinic_wave_equil_conservation.yml b/config/model_configs/baroclinic_wave_equil_conservation.yml
@@ -2,6 +2,7 @@ initial_condition: "MoistBaroclinicWave"
 dt: "400secs"
 t_end: "5days"
 moist: "equil"
+surface_setup: "NoSurface"
 dt_save_state_to_disk: "5days"
 check_conservation: true
 diagnostics:

diff --git a/config/model_configs/baroclinic_wave_equil_conservation_ft64.yml b/config/model_configs/baroclinic_wave_equil_conservation_ft64.yml
@@ -3,6 +3,7 @@ initial_condition: "MoistBaroclinicWave"
 dt: "400secs"
 t_end: "5days"
 moist: "equil"
+surface_setup: "NoSurface"
 dt_save_state_to_disk: "5days"
 check_conservation: true
 diagnostics:

diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
@@ -116,6 +116,7 @@ include(
 include(
     joinpath("prognostic_equations", "vertical_diffusion_boundary_layer.jl"),
 )
+include(joinpath("prognostic_equations", "surface_flux.jl"))
 include(joinpath("parameterized_tendencies", "sponge", "rayleigh_sponge.jl"))
 include(joinpath("parameterized_tendencies", "sponge", "viscous_sponge.jl"))
 include(

diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
@@ -36,7 +36,7 @@ function ClimaAtmosParameters(toml_dict::TD) where {TD <: CP.AbstractTOMLDict}
     IP = typeof(insolation_params)
 
     surface_fluxes_params =
-        SF.Parameters.SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
+        SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
     SFP = typeof(surface_fluxes_params)
 
     surface_temp_params = SurfaceTemperatureParameters(toml_dict)

diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -178,7 +178,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; sfc_conditions) = p.precomputed
     (; ᶜρa⁰, ᶜu⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰, ᶜtke⁰, ᶜmixing_length) = p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
@@ -192,7 +191,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
         # energy
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+            bottom = Operators.SetValue(C3(FT(0))),
         )
         @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜmse⁰ + ᶜK⁰)))
         if use_prognostic_tke(turbconv_model)
@@ -216,7 +215,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
             ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
             ᶜdivᵥ_ρq_tot = Operators.DivergenceF2C(
                 top = Operators.SetValue(C3(FT(0))),
-                bottom = Operators.SetValue(sfc_conditions.ρ_flux_q_tot),
+                bottom = Operators.SetValue(C3(FT(0))),
             )
             @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρq_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜq_tot⁰)))
             @. Yₜ.c.ρq_tot -= ᶜρχₜ_diffusion
@@ -228,12 +227,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
         bc_strain_rate = compute_strain_rate_face(ᶜu⁰)
         @. ᶠstrain_rate = bc_strain_rate
         @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(-(2 * ᶠρaK_u * ᶠstrain_rate)) / Y.c.ρ)
-        # apply boundary condition for momentum flux
-        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
-            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
-        )
-        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
     end
 
     # TODO: Add tracer flux
@@ -253,7 +246,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; sfc_conditions) = p.precomputed
     (; ᶜu, ᶜh_tot, ᶜspecific, ᶜtke⁰, ᶜmixing_length) = p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
@@ -267,7 +259,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
         # energy
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+            bottom = Operators.SetValue(C3(FT(0))),
         )
         @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜh_tot)))
 
@@ -288,7 +280,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
             ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
             ᶜdivᵥ_ρq_tot = Operators.DivergenceF2C(
                 top = Operators.SetValue(C3(FT(0))),
-                bottom = Operators.SetValue(sfc_conditions.ρ_flux_q_tot),
+                bottom = Operators.SetValue(C3(FT(0))),
             )
             @. ᶜρχₜ_diffusion =
                 ᶜdivᵥ_ρq_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜspecific.q_tot)))
@@ -301,12 +293,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
         bc_strain_rate = compute_strain_rate_face(ᶜu)
         @. ᶠstrain_rate = bc_strain_rate
         @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(-(2 * ᶠρaK_u * ᶠstrain_rate)) / Y.c.ρ)
-        # apply boundary condition for momentum flux
-        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
-            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
-        )
-        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
     end
 
     # TODO: Add tracer flux

diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl
@@ -109,6 +109,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
         edmfx_sgs_diffusive_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end
 
+    surface_flux_tendency!(Yₜ, Y, p, t, p.atmos.surface_flux_model)
+
     radiation_tendency!(Yₜ, Y, p, t, p.atmos.radiation_mode)
     edmfx_entr_detr_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     if p.atmos.sgs_mf_mode == Explicit()

diff --git a/src/prognostic_equations/surface_flux.jl b/src/prognostic_equations/surface_flux.jl
@@ -0,0 +1,51 @@
+#####
+##### Surface flux tendencies
+#####
+
+import ClimaCore.Geometry: ⊗
+import ClimaCore.Operators as Operators
+
+surface_flux_tendency!(Yₜ, Y, p, t) =
+    surface_flux_tendency!(Yₜ, Y, p, t, p.atmos.surface_flux_model)
+
+surface_flux_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing
+
+function surface_flux_tendency!(Yₜ, Y, p, t, _)
+    FT = eltype(Y)
+    (; ᶜh_tot, ᶜspecific, sfc_conditions) = p.precomputed
+    ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
+
+    if diffuse_momentum(p.atmos.vert_diff) || !isnothing(p.atmos.turbconv_model)
+        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
+            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
+        )
+        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
+    end
+
+    ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
+        top = Operators.SetValue(C3(FT(0))),
+        bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+    )
+    @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(FT(0) * ᶠgradᵥ(ᶜh_tot)))
+
+    ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
+    ρ_flux_χ = p.scratch.sfc_temp_C3
+    for (ᶜρχₜ, ᶜχ, χ_name) in matching_subfields(Yₜ.c, ᶜspecific)
+        χ_name == :e_tot && continue
+        if χ_name == :q_tot
+            @. ρ_flux_χ = sfc_conditions.ρ_flux_q_tot
+        else
+            @. ρ_flux_χ = C3(FT(0))
+        end
+        ᶜdivᵥ_ρχ = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0))),
+            bottom = Operators.SetValue(ρ_flux_χ),
+        )
+        @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρχ(-(FT(0) * ᶠgradᵥ(ᶜχ)))
+        @. ᶜρχₜ -= ᶜρχₜ_diffusion
+        if !(χ_name in (:q_rai, :q_sno))
+            @. Yₜ.c.ρ -= ᶜρχₜ_diffusion
+        end
+    end
+end
diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
@@ -2,15 +2,7 @@
 ##### Vertical diffusion boundary layer parameterization
 #####
 
-import StaticArrays
 import ClimaCore.Geometry: ⊗
-import ClimaCore.Utilities: half
-import LinearAlgebra: norm
-import Thermodynamics as TD
-import SurfaceFluxes as SF
-import ClimaCore.Spaces as Spaces
-import ClimaCore.Geometry as Geometry
-import ClimaCore.Fields as Fields
 import ClimaCore.Operators as Operators
 
 vertical_diffusion_boundary_layer_tendency!(Yₜ, Y, p, t) =
@@ -36,34 +28,21 @@ function vertical_diffusion_boundary_layer_tendency!(
         @. Yₜ.c.uₕ -= C12(
             ᶜdivᵥ(-2 * ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_u) * ᶠstrain_rate) / Y.c.ρ,
         )
-
-        # apply boundary condition for momentum flux
-        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
-            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
-        )
-        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
     end
 
     ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
         top = Operators.SetValue(C3(FT(0))),
-        bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+        bottom = Operators.SetValue(C3(FT(0))),
     )
     @. Yₜ.c.ρe_tot -=
         ᶜdivᵥ_ρe_tot(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜh_tot)))
 
     ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
-    ρ_flux_χ = p.scratch.sfc_temp_C3
     for (ᶜρχₜ, ᶜχ, χ_name) in matching_subfields(Yₜ.c, ᶜspecific)
         χ_name == :e_tot && continue
-        if χ_name == :q_tot
-            @. ρ_flux_χ = sfc_conditions.ρ_flux_q_tot
-        else
-            @. ρ_flux_χ = C3(FT(0))
-        end
         ᶜdivᵥ_ρχ = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
-            bottom = Operators.SetValue(ρ_flux_χ),
+            bottom = Operators.SetValue(C3(FT(0))),
         )
         @. ᶜρχₜ_diffusion =
             ᶜdivᵥ_ρχ(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜχ)))

diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
@@ -107,6 +107,31 @@ function get_surface_model(parsed_args)
     end
 end
 
+function get_surface_flux_model(parsed_args)
+    surface_setup_name = parsed_args["surface_setup"]
+    return if surface_setup_name == "NoSurface"
+        nothing
+    elseif surface_setup_name in (
+        "PrescribedSurface",
+        "DefaultExchangeCoefficients",
+        "DefaultMoninObukhov",
+        "GABLS",
+        "Soares",
+        "Bomex",
+        "DYCOMS_RF01",
+        "DYCOMS_RF02",
+        "Rico",
+        "TRMM_LBA",
+        "SimplePlume",
+        "ISDAC",
+        "GCM",
+    )
+        SurfaceFlux()
+    else
+        error("Uncaught surface flux model `$surface_setup`.")
+    end
+end
+
 function get_surface_albedo_model(parsed_args, params, ::Type{FT}) where {FT}
     albedo_name = parsed_args["albedo_model"]
     return if albedo_name in ("ConstantAlbedo",)

diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
@@ -106,6 +106,7 @@ function get_atmos(config::AtmosConfig, params)
         sfc_temperature = get_sfc_temperature_form(parsed_args),
         insolation = get_insolation_form(parsed_args),
         surface_model = get_surface_model(parsed_args),
+        surface_flux_model = get_surface_flux_model(parsed_args),
         surface_albedo = get_surface_albedo_model(parsed_args, params, FT),
         numerics = get_numerics(parsed_args),
     )

diff --git a/src/solver/types.jl b/src/solver/types.jl
@@ -219,6 +219,8 @@ end
 diffuse_momentum(::DecayWithHeightDiffusion{DM}) where {DM} = DM
 diffuse_momentum(::Nothing) = false
 
+struct SurfaceFlux end
+
 abstract type AbstractSponge end
 Base.Broadcast.broadcastable(x::AbstractSponge) = tuple(x)
 Base.@kwdef struct ViscousSponge{FT} <: AbstractSponge
@@ -523,6 +525,7 @@ Base.@kwdef struct AtmosModel{
     ST,
     IN,
     SM,
+    SFM,
     SA,
     NUM,
 }
@@ -561,6 +564,7 @@ Base.@kwdef struct AtmosModel{
     sfc_temperature::ST = nothing
     insolation::IN = nothing
     surface_model::SM = nothing
+    surface_flux_model::SFM = nothing
     surface_albedo::SA = nothing
     numerics::NUM = nothing
 end

diff --git a/src/surface_conditions/surface_setups.jl b/src/surface_conditions/surface_setups.jl
@@ -1,5 +1,13 @@
 abstract type SurfaceSetup end
 
+"""
+    NoSurface()
+
+Used to indicate that there is no surface flux.
+"""
+struct NoSurface <: SurfaceSetup end
+(surface_setup::NoSurface)(params) = nothing
+
 """
     PrescribedSurface()