mutable struct PlanktonInput
    temp::AbstractArray{Float64,4}      # temperature
    PARF::AbstractArray{Float64,3}      # PARF
    vels::NamedTuple                    # velocity fields for nutrients and individuals
    ΔT_vel::Int                         # time step of velocities provided
    ΔT_PAR::Int                         # time step of surface PAR provided
    ΔT_temp::Int                        # time step of temperature provided
end

mutable struct PlanktonSimulation
    model::PlanktonModel                       # Model object
    input::PlanktonInput                       # model input, temp, PAR, and velocities
    diags::Union{PlanktonDiagnostics,Nothing}  # diagnostics
    ΔT::Int64                                  # model time step
    iterations::Int64                          # run the simulation for this number of iterations
    output_writer::Union{PlanktonOutputWriter,Nothing} # Output writer
end

"""
    PlanktonSimulation(model; ΔT, iterations,
                       PARF = default_PARF(model.grid),
                       temp = default_temperature(model.grid),
                       diags = nothing,
                       vels = (;),
                       ΔT_vel = ΔT,
                       ΔT_PAR::Int64 = 3600,
                       ΔT_temp::Int64 = 3600,
                       output_writer = nothing,
                       )
Generate a `PlanktonSimulation` data structure.

Keyword Arguments (Required)
============================
- `ΔT` : time step in second.
- `iterations` : run the simulation for this number of iterations.

Keyword Arguments (Optional)
============================
- `PARF` : External forcings of surface PAR. Hourly PAR of a single day is provided by default.
- `temp` : External forcings of temperature. Hourly data of a single day is provided by default.
- `diags` : Diagnostics of the simulation generated by `PlanktonDiagnostics`.
- `vels` : The velocity fields for nutrient fields and individuals. `nothing` means no velocities
                     will be applied in the simulation. Otherwise, `vels` mush be a `NamedTuple` containing
                     all `u`, `v`, and `w`. Each of `u`, `v`, and `w` must be an 4D-`Array` of 
                     `(Nx, Ny, Nz, nΔT)` elements, excluding halo points. `N+1` is required for bounded direction.
- `ΔT_vel` : time step of velocities provided externally (in seconds).
- `ΔT_PAR` : time step of surface PAR provided externally (in seconds).
- `ΔT_temp` : time step of temperature provided externally (in seconds).
- `output_writer` : Output writer of the simulation generated by `PlanktonOutputWriter`.
"""
function PlanktonSimulation(model::PlanktonModel; ΔT::Int64, iterations::Int64,
                            PARF = default_PARF(model.grid, ΔT, iterations),
                            temp = default_temperature(model.grid, ΔT, iterations),
                            diags = nothing,
                            vels = (;),
                            ΔT_vel::Int64 = ΔT,
                            ΔT_PAR::Int64 = 3600,
                            ΔT_temp::Int64 = 3600,
                            output_writer = nothing
                            )

    if vels ≠ (;)
        if validate_velocity(vels, model.grid)
            vel_copy!(model.timestepper.vel₀, vels.u[:,:,:,1], vels.v[:,:,:,1], vels.w[:,:,:,1], model.grid)
        end

        @assert ΔT_vel % ΔT == 0.0

        if size(vels.u)[4] < floor(Int,iterations*ΔT/ΔT_vel)
            throw(ArgumentError("Velocities provided cannot cover the time period of (iterations+1)*ΔT ($(iterations*ΔT) seconds) with ΔT_vel = $(ΔT_vel)."))
        end
    end

    if size(PARF)[3] < floor(Int,iterations*ΔT/ΔT_PAR)
        throw(ArgumentError("Surface PAR provided cannot cover the time period of (iterations+1)*ΔT ($(iterations*ΔT) seconds) with ΔT_PAR = $(ΔT_PAR)."))
    end

    if size(temp)[4] < floor(Int,iterations*ΔT/ΔT_temp)
        throw(ArgumentError("Temperature provided cannot cover the time period of (iterations+1)*ΔT ($(iterations*ΔT) seconds) with ΔT_temp = $(ΔT_temp)."))
    end

    input = PlanktonInput(temp, PARF, vels, ΔT_vel, ΔT_PAR, ΔT_temp)

    validate_bcs(model.nutrients, model.grid, iterations)

    sim = PlanktonSimulation(model, input, diags, ΔT, iterations, output_writer)

    return sim
end

function vels_string(sim::PlanktonSimulation)
    if sim.input.vels == (;)
        return "no velocity provided"
    else
        return "velocity provided, one time step of velocity covers $(sim.input.ΔT_vel/sim.ΔT) time steps of simulation"
    end
end

function show(io::IO, sim::PlanktonSimulation)
    print(io, "PlanktonSimulation:\n",
              "├── ΔT: $(sim.ΔT)s\n",
              "├── model time: $(sim.model.t)s\n",
              "├── simulation stops after: $(sim.iterations) iterations\n",
              "├── $(vels_string(sim))\n",
              "├── diagnostics: $(short_show(sim.diags))\n",
              "└── output writer: $(short_show(sim.output_writer))"
              )

end

function validate_velocity(vels, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
    u_size = (g.Nx, g.Ny, g.Nz)
    v_size = (g.Nx, g.Ny, g.Nz)
    w_size = (g.Nx, g.Ny, g.Nz)
    validation = true

    if isa(TX(), Bounded)
        u_size = (g.Nx+1, g.Ny, g.Nz)
    end

    if isa(TY(), Bounded)
        v_size = (g.Nx, g.Ny+1, g.Nz)
    end

    if isa(TZ(), Bounded)
        w_size = (g.Nx, g.Ny, g.Nz+1)
    end

    if size(vels.u)[1:3] ≠ u_size
        validation = false
        throw(ArgumentError("Dimension mismatch: the size of u must be $(u_size), for $(TX) topology."))
    end

    if size(vels.v)[1:3] ≠ v_size
        validation = false
        throw(ArgumentError("Dimension mismatch: the size of v must be $(v_size), for $(TY) topology."))
    end

    if size(vels.w)[1:3] ≠ w_size
        validation = false
        throw(ArgumentError("Dimension mismatch: the size of w must be $(w_size), for $(TZ) topology."))
    end

    return validation
end

function short_show(diags::PlanktonDiagnostics)
    return string("\n", "│   ├── diagnostics of tracers: $(keys(diags.tracer))\n", 
                        "│   ├── diagnostics of individuals: $(keys(diags.plankton.sp1))\n",
                        "│   └── save averaged diagnostics every $(diags.time_interval) seconds")
end
function short_show(writer::PlanktonOutputWriter)
    return string("\n", "    ├── files are saved at $(writer.filepath)\n", 
                        "    └── Maximum file size: $(humanize_filesize(writer.max_filesize))")
end
function short_show(nothing::Nothing)
    return string("nothing")
end
