export LandDomainBC,
    LandComponentBC,
    Dirichlet,
    Neumann,
    NoBC,
    SurfaceDrivenWaterBoundaryConditions,
    SurfaceDrivenHeatBoundaryConditions

"""
   AbstractBoundaryConditions 
"""
abstract type AbstractBoundaryConditions end


"""
    NoBC <: AbstractBoundaryConditions

This type is used for dispatch when no boundary condition needs
to be enforced - for example, if no prognostic variables are included
for a subcomponent, or for lateral faces when a 1D vertical setup
is used.
"""
struct NoBC <: AbstractBoundaryConditions end

"""
    Dirichlet{Fs} <: AbstractBoundaryConditions

A concrete type to hold the state variable 
function, if Dirichlet boundary conditions are desired.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Dirichlet{Fs} <: AbstractBoundaryConditions
    "state boundary condition"
    state_bc::Fs
end

"""
    Neumann{Ff} <: AbstractBoundaryConditions

A concrete type for specifying the magnitude of an inward diffusive flux,
normal to the domain boundary, from which a Neumann boundary condition 
is determined and applied.

For the soil water model, the Darcy flux is defined as F⃗(Darcy) = -K∇h. 
With the `Neumann` boundary
condition type, the user supplies a scalar `f` such that the applied flux
boundary conditon is F⃗(Darcy, boundary) = -f n̂, 
where n̂ is the normal vector pointing out of the domain. For the soil heat model,
the same rule applies, except now we have a conductive heat flux 
F⃗(Conductive, boundary) = -κ∇T = -f n̂.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Neumann{Ff} <: AbstractBoundaryConditions
    "Scalar flux boundary condition"
    scalar_flux_bc::Ff
end

"""
    SurfaceDrivenWaterBoundaryConditions{FT, PD, RD} <: AbstractBoundaryConditions

Boundary condition type to be used when the user wishes to 
apply physical fluxes of water at the top of the domain, normal to the surface.
The flux can be inwards or outwards, depending on the magnitude of infiltration,
evaporation, and precipitation.

Precipitation is assumed to be 
in the ẑ direction, P⃗ = Pẑ, while evaporation is along the normal 
direction, E⃗ = En̂. The applied boundary flux is F⃗ = f n̂, and the magnitude f is
determined internally from P⃗, E⃗, and soil conditions.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SurfaceDrivenWaterBoundaryConditions{FT, PD, RD} <:
       AbstractBoundaryConditions where {FT, PD, RD}
    "Precipitation model"
    precip_model::PD
    "Runoff model"
    runoff_model::RD
end

"""
    SurfaceDrivenWaterBoundaryConditions(
        ::Type{FT};
        precip_model::AbstractPrecipModel{FT} = DrivenConstantPrecip{FT}(),
        runoff_model::AbstractSurfaceRunoffModel{FT} = NoRunoff(),
    ) where {FT}

Constructor for the SurfaceDrivenWaterBoundaryConditions object. The default
is a constant precipitation rate on the subgrid scale, and no runoff.
"""
function SurfaceDrivenWaterBoundaryConditions(
    ::Type{FT};
    precip_model::AbstractPrecipModel{FT} = DrivenConstantPrecip{FT}(
        (t) -> (0.0),
    ),
    runoff_model::AbstractSurfaceRunoffModel = NoRunoff(),
) where {FT}
    args = (precip_model, runoff_model)
    return SurfaceDrivenWaterBoundaryConditions{FT, typeof.(args)...}(args...)
end

"""
    SurfaceDrivenHeatBoundaryConditions{FT, SWD} <: AbstractBoundaryConditions

Boundary condition type to be used when the user wishes to 
apply normal fluxes of heat at the top of the domain 
(according to radiative energy fluxes).

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SurfaceDrivenHeatBoundaryConditions{FT, SWD} <:
       AbstractBoundaryConditions where {FT, SWD}
    "Net short wave flux model"
    nswf_model::SWD
end

"""
    SurfaceDrivenHeatBoundaryConditions(
        ::Type{FT};
        nswf_model::AbstractNetSwFluxModel{FT} = PrescribedNetSwFlux{FT}()
    ) where {FT}

Constructor for the SurfaceDrivenHeatBoundaryConditions object. The default
is no shortwave flux.
"""
function SurfaceDrivenHeatBoundaryConditions(
    ::Type{FT};
    nswf_model::AbstractNetSwFluxModel{FT} = PrescribedNetSwFlux(
        FT;
        nswf = t -> eltype(t)(0),
    ),
) where {FT}
    return SurfaceDrivenHeatBoundaryConditions{FT, typeof(nswf_model)}(
        nswf_model,
    )
end

"""
    LandDomainBC{TBC, BBC, LBC}

A container for the land boundary conditions, with options for surface
boundary conditions, bottom boundary conditions, or lateral face
boundary conditions.

The user should supply an instance of `LandComponentBC` for each piece, as needed.
If none is supplied, the default for is `NoBC` for each subcomponent. At a minimum, both top and bottom boundary conditions should be supplied. Whether or not to include the lateral faces depends on the configuration of the domain.
"""
Base.@kwdef struct LandDomainBC{TBC, BBC, MinXBC, MaxXBC, MinYBC, MaxYBC}
    "surface boundary conditions"
    surface_bc::TBC = LandComponentBC()
    "bottom boundary conditions"
    bottom_bc::BBC = LandComponentBC()
    "lateral boundary conditions"
    minx_bc::MinXBC = LandComponentBC()
    maxx_bc::MaxXBC = LandComponentBC()
    miny_bc::MinYBC = LandComponentBC()
    maxy_bc::MaxYBC = LandComponentBC()
end

"""
    LandComponentBC{SW, SH, SF}

An object that holds the boundary conditions for each of the subcomponents
of the land model. 

The boundary conditions supplied should be of type `AbstractBoundaryConditions`. 
The default is `NoBC` for each component, so that the user only
needs to define the BC for the components they wish to model.
"""
Base.@kwdef struct LandComponentBC{
    SW <: AbstractBoundaryConditions,
    SH <: AbstractBoundaryConditions,
    SF <: AbstractBoundaryConditions,
}
    soil_water::SW = NoBC()
    soil_heat::SH = NoBC()
    surface::SF = NoBC()
end

"""
    function boundary_conditions(land::LandModel)

Unpacks the `boundary_conditions` field of the land model, and
puts into the correct order based on the integers used to identify
faces, as defined in the Driver configuration.
"""
function boundary_conditions(land::LandModel)
    bc = land.boundary_conditions
    mytuple = (
        bc.bottom_bc,
        bc.surface_bc,
        bc.minx_bc,
        bc.maxx_bc,
        bc.miny_bc,
        bc.maxy_bc,
    )
    # faces labeled integer 1,2 are bottom, top, lateral sides are 3, 4, 5, 6
    return mytuple
end

function boundary_state!(
    nf,
    bc::LandComponentBC,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    n,
    state⁻,
    aux⁻,
    t,
    args...,
)
    land_boundary_state!(
        nf,
        bc,
        land,
        state⁺,
        aux⁺,
        n,
        state⁻,
        aux⁻,
        t,
        args...,
    )
end

function land_boundary_state!(nf, bc::LandComponentBC, land, args...)
    soil_boundary_state!(nf, bc.soil_water, land.soil.water, land, args...)
    soil_boundary_state!(nf, bc.soil_heat, land.soil.heat, land, args...)
    surface_boundary_state!(nf, bc.surface, land.surface, land, args...)
end

function boundary_state!(
    nf,
    bc::LandComponentBC,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    hyperdiff⁺::Vars,
    aux⁺::Vars,
    n,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    t,
    args...,
)
    land_boundary_flux!(
        nf,
        bc,
        land,
        state⁺,
        diff⁺,
        aux⁺,
        n,
        state⁻,
        diff⁻,
        aux⁻,
        t,
        args...,
    )
end

function land_boundary_flux!(nf, bc::LandComponentBC, land, args...)
    soil_boundary_flux!(nf, bc.soil_water, land.soil.water, land, args...)
    soil_boundary_flux!(nf, bc.soil_heat, land.soil.heat, land, args...)
    surface_boundary_flux!(nf, bc.surface, land.surface, land, args...)
end
