Try #2174:

.buildkite/pipeline.yml

@@ -762,6 +762,15 @@ steps:
       queue: central
       slurm_ntasks: 1
 
+  - label: "cpu_land_overland_flow_vcatchment"
+    key: "cpu_land_overland_flow_vcatchment"
+    command:
+      - "mpiexec julia --color=yes --project test/Land/Model/test_overland_flow_vcatchment.jl"
+    agents:
+      config: cpu
+      queue: central
+      slurm_ntasks: 1
+
   - label: "gpu_thermodynamics"
     key: "gpu_thermodynamics"
     command:
@@ -1992,6 +2001,16 @@ steps:
       slurm_ntasks: 1
       slurm_gres: "gpu:1"
 
+  - label: "gpu_land_overland_flow_vcatchment"
+    key: "gpu_land_overland_flow_vcatchment"
+    command:
+      - "mpiexec julia --color=yes --project test/Land/Model/test_overland_flow_vcatchment.jl"
+    agents:
+      config: gpu
+      queue: central
+      slurm_ntasks: 1
+      slurm_gres: "gpu:1"
+
   - label: "gpu_unittests"
     key: "gpu_unittests"
     command:

.gitignore

@@ -5,6 +5,7 @@
 *.jl.*.cov
 *.jl.mem
 *.DS_Store
+*.vscode
 *~
 TAGS
 

docs/list_of_apis.jl

@@ -22,6 +22,9 @@ apis = [
             "APIs/Land/SoilWaterParameterizations.md",
         "Soil Heat Parameterizations" =>
             "APIs/Land/SoilHeatParameterizations.md",
+        "Surface Flow" => "APIs/Land/SurfaceFlow.md",
+        "Radiative Boundary Conditions" =>
+            "APIs/Land/RadiativeEnergyFlux.md",
     ],
     "Common" => [
         "Orientations" => "APIs/Common/Orientations.md",

docs/src/APIs/Land/SurfaceFlow.md

@@ -0,0 +1,17 @@
+# SurfaceFlow
+
+```@meta
+CurrentModule = ClimateMachine.Land.SurfaceFlow
+```
+
+## SurfaceFlow types and functions
+```@docs
+OverlandFlowModel
+NoSurfaceFlowModel
+surface_boundary_flux!
+surface_boundary_state!
+calculate_velocity
+Precip
+VolumeAdvection
+SurfaceWaterArea
+```

src/Driver/driver_configs.jl

@@ -686,7 +686,7 @@ function MultiColumnLandModel(
     zmin = zero(FT),
     array_type = ClimateMachine.array_type(),
     mpicomm = MPI.COMM_WORLD,
-    boundary = ((3, 3), (3, 3), (1, 2)),
+    boundary = ((3, 4), (5, 6), (1, 2)),
     solver_type = ExplicitSolverType(),
     periodicity = (false, false, false),
     meshwarp = (x...) -> identity(x),

src/Land/Model/LandModel.jl

@@ -30,7 +30,8 @@ import ..BalanceLaws:
     compute_gradient_flux!,
     nodal_init_state_auxiliary!,
     init_state_prognostic!,
-    nodal_update_auxiliary_state!
+    nodal_update_auxiliary_state!,
+    wavespeed
 using ..DGMethods: LocalGeometry, DGModel
 export LandModel
 
@@ -54,11 +55,13 @@ Users may over-ride prescribed default values for each field.
 # Fields
 $(DocStringExtensions.FIELDS)
 """
-struct LandModel{PS, S, LBC, SRC, SRCDT, IS} <: BalanceLaw
+struct LandModel{PS, S, SF, LBC, SRC, SRCDT, IS} <: BalanceLaw
     "Parameter set"
     param_set::PS
     "Soil model"
     soil::S
+    "Surface Flow model"
+    surface::SF
     "struct of boundary conditions"
     boundary_conditions::LBC
     "Source Terms (Problem specific source terms)"
@@ -75,6 +78,7 @@ parameter_set(m::LandModel) = m.param_set
     LandModel(
         param_set::AbstractParameterSet,
         soil::BalanceLaw;
+        surface::BalanceLaw = NoSurfaceFlowModel(),
         boundary_conditions::LBC = (),
         source::SRC = (),
         init_state_prognostic::IS = nothing
@@ -85,6 +89,7 @@ Constructor for the LandModel structure.
 function LandModel(
     param_set::AbstractParameterSet,
     soil::BalanceLaw;
+    surface::BalanceLaw = NoSurfaceFlowModel(),
     boundary_conditions::LBC = LandDomainBC(),
     source::SRC = (),
     init_state_prognostic::IS = nothing,
@@ -94,6 +99,7 @@ function LandModel(
     land = (
         param_set,
         soil,
+        surface,
         boundary_conditions,
         source,
         source_dt,
@@ -106,26 +112,32 @@ end
 function vars_state(land::LandModel, st::Prognostic, FT)
     @vars begin
         soil::vars_state(land.soil, st, FT)
+        surface::vars_state(land.surface, st, FT)
     end
 end
 
 
 function vars_state(land::LandModel, st::Auxiliary, FT)
     @vars begin
+        x::FT
+        y::FT
         z::FT
         soil::vars_state(land.soil, st, FT)
+        surface::vars_state(land.surface, st, FT)
     end
 end
 
 function vars_state(land::LandModel, st::Gradient, FT)
     @vars begin
         soil::vars_state(land.soil, st, FT)
+        surface::vars_state(land.surface, st, FT)
     end
 end
 
 function vars_state(land::LandModel, st::GradientFlux, FT)
     @vars begin
         soil::vars_state(land.soil, st, FT)
+        surface::vars_state(land.surface, st, FT)
     end
 end
 
@@ -135,8 +147,11 @@ function nodal_init_state_auxiliary!(
     tmp::Vars,
     geom::LocalGeometry,
 )
+    aux.x = geom.coord[1]
+    aux.y = geom.coord[2]
     aux.z = geom.coord[3]
     land_init_aux!(land, land.soil, aux, geom)
+    land_init_aux!(land, land.surface, aux, geom)
 end
 
 function flux_first_order!(
@@ -145,8 +160,19 @@ function flux_first_order!(
     state::Vars,
     aux::Vars,
     t::Real,
-    directions,
-) end
+    direction,
+)
+    tend = Flux{FirstOrder}()
+    _args = (; state, aux, t, direction)
+    args = merge(_args, (precomputed = precompute(land, _args, tend),))
+
+    map(prognostic_vars(land)) do prog
+        var, name = get_prog_state(flux, prog)
+        val = Σfluxes(prog, eq_tends(prog, land, tend), land, args)
+        setproperty!(var, name, val)
+    end
+    nothing
+end
 
 
 function compute_gradient_argument!(
@@ -210,6 +236,7 @@ function nodal_update_auxiliary_state!(
     t::Real,
 )
     land_nodal_update_auxiliary_state!(land, land.soil, state, aux, t)
+    land_nodal_update_auxiliary_state!(land, land.surface, state, aux, t)
 end
 
 
@@ -246,8 +273,8 @@ function init_state_prognostic!(
     land.init_state_prognostic(land, state, aux, coords, t, args...)
 end
 
-include("prog_types.jl")
 
+include("prog_types.jl")
 include("RadiativeEnergyFlux.jl")
 using .RadiativeEnergyFlux
 include("SoilWaterParameterizations.jl")
@@ -260,11 +287,22 @@ include("soil_heat.jl")
 include("Runoff.jl")
 using .Runoff
 include("land_bc.jl")
+include("SurfaceFlow.jl")
+using .SurfaceFlow
 include("soil_bc.jl")
-
 include("prognostic_vars.jl")
-
-include("source.jl")
 include("land_tendencies.jl")
+include("source.jl")
+
+function wavespeed(m::LandModel, n⁻, state::Vars, aux::Vars, t::Real, direction)
+    FT = eltype(state)
+    g = FT(9.8)
+    width = m.surface.width(aux.x, aux.y)
+    area = max(eltype(state)(0.0), state.surface.area)
+    height = area ./ width
+    v = calculate_velocity(m.surface, aux.x, aux.y, height)
+    speed = abs(norm(v))
+    return speed
+end
 
 end # Module