Make surface flux fully explicit

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

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

config/model_configs/baroclinic_wave.yml

@@ -3,6 +3,7 @@ initial_condition: "DryBaroclinicWave"
 dt: "400secs"
 t_end: "10days"
 deep_atmosphere: false
+surface_setup: "NoSurface"
 diagnostics:
   - short_name: [pfull, ua, wa, va, rv, ta, ke]
     period: 1days

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]

config/model_configs/baroclinic_wave_equil.yml

@@ -5,6 +5,7 @@ initial_condition: "MoistBaroclinicWave"
 dt: "450secs"
 t_end: "10days"
 moist: "equil"
+surface_setup: "NoSurface"
 diagnostics:
   - short_name: [pfull, wa, va, rv, hus, ke]
     period: 1days

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:

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:

config/model_configs/baroclinic_wave_hughes2023.yml

@@ -3,6 +3,7 @@ initial_condition: "MoistBaroclinicWave"
 topography: "Hughes2023"
 topo_smoothing: false
 moist: equil
+surface_setup: "NoSurface"
 dt: "200secs"
 t_end: "10days"
 diagnostics:

config/model_configs/baroclinic_wave_topography_dcmip_rs.yml

@@ -6,6 +6,7 @@ initial_condition: "DryBaroclinicWave"
 t_end: "6days"
 dt: "200secs"
 hyperdiff: "ClimaHyperdiffusion"
+surface_setup: "NoSurface"
 netcdf_output_at_levels: false
 diagnostics:
   - short_name: [pfull, wa, va, rv]

config/model_configs/box_cosine_hills_float64_test.yml

@@ -12,6 +12,7 @@ dz_bottom: 10
 dt: "0.5secs" # CFL is slightly lower than for 2D case
 t_end: "30mins"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/box_density_current_test.yml

@@ -14,6 +14,7 @@ z_stretch: false
 dt: "0.25secs"
 t_end: "20secs"
 dt_save_state_to_disk: "5secs"
+surface_setup: "NoSurface"
 netcdf_interpolation_num_points: [40, 40, 80]
 diagnostics:
   - short_name: [wa, ua, va, ta, thetaa, ha]

config/model_configs/box_hydrostatic_balance.yml

@@ -3,6 +3,7 @@ config: "box"
 hyperdiff: false
 dt: "20secs"
 perturb_initstate: false
+surface_setup: "NoSurface"
 diagnostics:
   - short_name: [wa, ua]
     period: 1days

config/model_configs/extruded_plane_cosine_hills_float64_test.yml

@@ -12,6 +12,7 @@ dz_bottom: 10
 dt: "0.5secs" # CFL is slightly lower than for 2D case
 t_end: "9hours"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/hydrostatic_balance_ft64.yml

@@ -4,4 +4,5 @@ t_end: "8days"
 discrete_hydrostatic_balance: true
 hyperdiff: false
 perturb_initstate: false
+surface_setup: "NoSurface"
 FLOAT_TYPE: "Float64"

config/model_configs/plane_agnesi_mountain_float64_test.yml

@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/plane_cosine_hills_float64_test.yml

@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/plane_density_current_test.yml

@@ -9,6 +9,7 @@ z_stretch: false
 x_elem: 80
 config: "plane"
 z_max: 6400.0
+surface_setup: "NoSurface"
 output_default_diagnostics: false
 diagnostics:
   - short_name: thetaa

config/model_configs/plane_no_topography_float64_test.yml

@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/plane_schar_mountain_float32_test.yml

@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/plane_schar_mountain_float64_test.yml

@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+surface_setup: "NoSurface"
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false

config/model_configs/single_column_hydrostatic_balance_ft64.yml

@@ -5,6 +5,7 @@ dt_save_state_to_disk: "5days"
 initial_condition: "IsothermalProfile"
 config: "column"
 hyperdiff: false
+surface_setup: "NoSurface"
 dt: "3hours"
 FLOAT_TYPE: "Float64"
 reproducibility_test: true

config/model_configs/single_column_nonorographic_gravity_wave.yml

@@ -4,4 +4,5 @@ t_end: "1500secs"
 config: "column"
 hyperdiff: false
 dt: "400secs"
+surface_setup: "NoSurface"
 non_orographic_gravity_wave: true

config/model_configs/single_column_radiative_equilibrium_allsky_idealized_clouds.yml

@@ -12,6 +12,7 @@ idealized_clouds: true
 config: "column"
 insolation: "timevarying"
 dt_save_to_sol: "30hours"
+surface_setup: "NoSurface"
 rad: "allskywithclear"
 diagnostics:
   - short_name: rhoa

config/model_configs/single_column_radiative_equilibrium_clearsky.yml

@@ -10,4 +10,5 @@ idealized_h2o: true
 t_end: "654days"
 config: "column"
 dt_save_to_sol: "30hours"
+surface_setup: "NoSurface"
 rad: "clearsky"

config/model_configs/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.yml

@@ -11,5 +11,6 @@ dt: "3hours"
 dt_rad: "3hours"
 dt_save_to_sol: "30hours"
 dt_save_state_to_disk: "100days"
+surface_setup: "NoSurface"
 prognostic_surface: true
 toml: [toml/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.toml]

config/model_configs/single_column_radiative_equilibrium_gray.yml

@@ -10,6 +10,7 @@ t_end: "654days"
 dz_bottom: 30.0
 config: "column"
 dt_save_to_sol: "30hours"
+surface_setup: "NoSurface"
 rad: "gray"
 # [2, 2, 80] instead of [1, 1, 80] because Julia ranges are inclusive of the
 # extrema. Given that our columns are 3D, we cannot map the horizontal dimension

config/model_configs/sphere_nonorographic_gravity_wave.yml

@@ -1,3 +1,4 @@
 t_end: "1500secs"
 dt: "400secs"
+surface_setup: "NoSurface"
 non_orographic_gravity_wave: true

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(

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)

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

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()

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