#-------------------------------------------------------------------------------------
# Create a function to run SNEASY-Hector and return model output used in calibration.
#-------------------------------------------------------------------------------------

function construct_run_sneasy_hectorch4(calibration_end_year::Int)

    # Load an instance of SNEASY-Hector.
    m = create_sneasy_hectorch4(start_year=1765, end_year=calibration_end_year, etminan_ch4_forcing=true)

    # F2x_CO₂ has a default value in the component definition and therefore no
    # parameter in the model is created. The following line creates a parameter
    # in the model that can then be updated with update_param! in the inner loop.
    set_param!(m, :doeclim, :F2x_CO₂, 3.7)

    # Get indices needed to normalize temperature anomalies relative to 1861-1880 mean.
    index_1861, index_1880 = findall((in)([1861, 1880]), collect(1765:calibration_end_year))

    mi = Mimi.build(m)

    # Given user settings, create a function to run SNEASY-Hector and return model output used for calibration.
    function run_sneasy_hectorch4!(
        params::Array{Float64,1},
        modeled_CO₂::Vector{Float64},
        modeled_CH₄::Vector{Float64},
        modeled_oceanCO₂_flux::Vector{Float64},
        modeled_temperature::Vector{Float64},
        modeled_ocean_heat::Vector{Float64})

        # Assign names to uncertain model and initial condition parameters for convenience.
        # Note: This assumes "params" is the full vector of uncertain parameters with the same ordering as in "create_log_posterior.jl".
        temperature_0     = params[9]
        ocean_heat_0      = params[10]
        CO₂_0             = params[11]
        CH₄_0             = params[12]
        N₂O_0             = params[13]
        ECS               = params[14]
        heat_diffusivity  = params[15]
        rf_scale_aerosol  = params[16]
        rf_scale_CH₄      = params[17]
        F2x_CO₂           = params[18]
        Q10               = params[19]
        CO₂_fertilization = params[20]
        CO₂_diffusivity   = params[21]
        τ_troposphere     = params[22]
        CH₄_natural       = params[23]
        τ_soil            = params[24]
        τ_stratosphere    = params[25]

        #----------------------------------------------------------
        # Set SNEASY-Hector with sampled parameter values.
        #----------------------------------------------------------

        # ---- Diffusion Ocean Energy balance CLIMate model (DOECLIM) ---- #
        update_param!(mi, :t2co, ECS)
        update_param!(mi, :kappa, heat_diffusivity)
        update_param!(mi, :F2x_CO₂, F2x_CO₂)

        # ---- Carbon Cycle ---- #
        update_param!(mi, :Q10, Q10)
        update_param!(mi, :Beta, CO₂_fertilization)
        update_param!(mi, :Eta, CO₂_diffusivity)
        update_param!(mi, :atmco20, CO₂_0)

        # ---- Tropospheric Sink (OH) Lifetime ---- #
        update_param!(mi, :TOH0, τ_troposphere)

        # ---- Methane Cycle ---- #
        update_param!(mi, :CH4N, CH₄_natural)
        update_param!(mi, :Tsoil, τ_soil)
        update_param!(mi, :Tstrat, τ_stratosphere)

        # ---- Methane Radiative Forcing ---- #
        update_param!(mi, :CH₄_0, CH₄_0)
        update_param!(mi, :scale_CH₄, rf_scale_CH₄)

        # ---- Carbon Dioxide Radiative Forcing ---- #
        update_param!(mi, :CO₂_0, CO₂_0)
        update_param!(mi, :rf_scale_CO₂, co2_rf_scale(F2x_CO₂, CO₂_0, N₂O_0))

        # ---- Total Radiative Forcing ---- #
        update_param!(mi, :α, rf_scale_aerosol)

        # ---- Shared parameters ---- #
        # TODO Ideally we wouldn't have the `M0` parameter here, but
        # instead use the already existing `CH₄_0` for this as well.
        update_param!(mi, :M0, CH₄_0)
        update_param!(mi, :N₂O_0, N₂O_0)

        # Run model.
        run(mi)

        #----------------------------------------------------------
        # Calculate model output being compared to observations.
        #----------------------------------------------------------

        # Atmospheric concentration of CO₂.
        modeled_CO₂[:] = mi[:ccm, :atmco2]

        # Ocean carbon flux (Note: timesteps cause last `atm_oc_flux` value to equal `missing`, so exclude it here).
        modeled_oceanCO₂_flux[1:end-1] = mi[:ccm, :atm_oc_flux][1:end-1]

        # Global surface temperature anomaly (normalized to 1861-1880 mean with initial condition offset).
        modeled_temperature[:] = mi[:doeclim, :temp] .- mean(mi[:doeclim, :temp][index_1861:index_1880]) .+ temperature_0

        # Ocean heat content (with initial condition offset).
        modeled_ocean_heat[:] = mi[:doeclim, :heat_mixed] .+ mi[:doeclim, :heat_interior] .+ ocean_heat_0

        # Atmospheric concentration of CH₄.
        modeled_CH₄[:] = mi[:ch4_cycle, :CH4]

        return
    end

    # Return function.
    return run_sneasy_hectorch4!
end
