def run_one_step_catchment(area_m2, time_gap_sec, c_in_rain, c_in_peva,
    c_p_t, c_p_c, c_p_h, c_p_d, c_p_s, c_p_z, c_p_sk, c_p_fk, c_p_gk,
    c_s_v_ove, c_s_v_dra, c_s_v_int, c_s_v_sgw, c_s_v_dgw, c_s_v_ly1,
    c_s_v_ly2, c_s_v_ly3, c_s_v_ly4, c_s_v_ly5, c_s_v_ly6):
    """
    This function was written by Thibault Hallouin but is largely inspired by the work of Eva Mockler, namely for
    the work published in: Mockler, E., O’Loughlin, F., and Bruen, M.: Understanding hydrological flow paths in
    conceptual catchment models using uncertainty and sensitivity analysis, Computers & Geosciences, 90, 66–77,
    doi:10.1016/j.cageo.2015.08.015, 2016.

    Catchment model * c_ *
    _ Hydrology
    ___ Inputs * in_ *
    _____ c_in_rain         precipitation as rain [mm/time step]
    _____ c_in_peva         potential evapotranspiration [mm/time step]
    ___ Parameters * p_ *
    _____ c_p_t             T: rainfall aerial correction coefficient
    _____ c_p_c             C: evaporation decay parameter
    _____ c_p_h             H: quick runoff coefficient
    _____ c_p_d             D: drain flow parameter - fraction of saturation excess diverted to drain flow
    _____ c_p_s             S: soil outflow coefficient
    _____ c_p_z             Z: effective soil depth [mm]
    _____ c_p_sk            SK: surface routing parameter [hours]
    _____ c_p_fk            FK: inter flow routing parameter [hours]
    _____ c_p_gk            GK: groundwater routing parameter [hours]
    ___ States * s_ *
    _____ c_s_v_ove         volume of water in overland store [m3]
    _____ c_s_v_dra         volume of water in drain store [m3]
    _____ c_s_v_int         volume of water in inter store [m3]
    _____ c_s_v_sgw         volume of water in shallow groundwater store [m3]
    _____ c_s_v_dgw         volume of water in deep groundwater store [m3]
    _____ c_s_v_ly1         volume of water in first soil layer store [m3]
    _____ c_s_v_ly2         volume of water in second soil layer store [m3]
    _____ c_s_v_ly3         volume of water in third soil layer store [m3]
    _____ c_s_v_ly4         volume of water in fourth soil layer store [m3]
    _____ c_s_v_ly5         volume of water in fifth soil layer store [m3]
    _____ c_s_v_ly6         volume of water in sixth soil layer store [m3]
    ___ Outputs * out_ *
    _____ c_out_aeva        actual evapotranspiration [m3/s]
    _____ c_out_q_ove       overland flow [m3/s]
    _____ c_out_q_dra       drain flow [m3/s]
    _____ c_out_q_int       inter flow [m3/s]
    _____ c_out_q_sgw       shallow groundwater flow [m3/s]
    _____ c_out_q_dgw       deep groundwater flow [m3/s]
    """
    nb_soil_layers = 6.0
    c_p_sk *= 3600.0
    c_p_fk *= 3600.0
    c_p_gk *= 3600.0
    list_z_lyr = [0.0, c_p_z / nb_soil_layers, c_p_z / nb_soil_layers, 
        c_p_z / nb_soil_layers, c_p_z / nb_soil_layers, c_p_z /
        nb_soil_layers, c_p_z / nb_soil_layers]
    list_lvl_lyr = [0.0, c_s_v_ly1 / area_m2 * 1000.0, c_s_v_ly2 / area_m2 *
        1000.0, c_s_v_ly3 / area_m2 * 1000.0, c_s_v_ly4 / area_m2 * 1000.0,
        c_s_v_ly5 / area_m2 * 1000.0, c_s_v_ly6 / area_m2 * 1000.0]
    lvl_total_start = sum(list_lvl_lyr)
    rain = c_in_rain * c_p_t
    excess_rain = rain - c_in_peva
    aeva = 0.0
    if excess_rain >= 0.0:
        aeva += c_in_peva
        h_prime = c_p_h * (lvl_total_start / c_p_z)
        overland_flow = h_prime * excess_rain
        excess_rain -= overland_flow
        for i in [1, 2, 3, 4, 5, 6]:
            space_in_lyr = list_z_lyr[i] - list_lvl_lyr[i]
            if excess_rain <= space_in_lyr:
                list_lvl_lyr[i] += excess_rain
                excess_rain = 0.0
            else:
                list_lvl_lyr[i] = list_z_lyr[i]
                excess_rain -= space_in_lyr
        drain_flow = c_p_d * excess_rain
        inter_flow = (1.0 - c_p_d) * excess_rain
        s_prime = c_p_s * (lvl_total_start / c_p_z)
        for i in [1, 2, 3, 4, 5, 6]:
            leak_interflow = list_lvl_lyr[i] * s_prime ** i
            if leak_interflow < list_lvl_lyr[i]:
                inter_flow += leak_interflow
                list_lvl_lyr[i] -= leak_interflow
        shallow_flow = 0.0
        for i in [1, 2, 3, 4, 5, 6]:
            leak_shallow_flow = list_lvl_lyr[i] * (s_prime / i)
            if leak_shallow_flow < list_lvl_lyr[i]:
                shallow_flow += leak_shallow_flow
                list_lvl_lyr[i] -= leak_shallow_flow
        deep_flow = 0.0
        for i in [6, 5, 4, 3, 2, 1]:
            leak_deep_flow = list_lvl_lyr[i] * s_prime ** (7 - i)
            if leak_deep_flow < list_lvl_lyr[i]:
                deep_flow += leak_deep_flow
                list_lvl_lyr[i] -= leak_deep_flow
    else:
        overland_flow = 0.0
        drain_flow = 0.0
        inter_flow = 0.0
        shallow_flow = 0.0
        deep_flow = 0.0
        deficit_rain = excess_rain * -1.0
        aeva += rain
        for i in [1, 2, 3, 4, 5, 6]:
            if list_lvl_lyr[i] >= deficit_rain:
                list_lvl_lyr[i] -= deficit_rain
                aeva += deficit_rain
                deficit_rain = 0.0
            else:
                aeva += list_lvl_lyr[i]
                deficit_rain = c_p_c * (deficit_rain - list_lvl_lyr[i])
                list_lvl_lyr[i] = 0.0
    c_out_aeva = aeva / 1000.0 * area_m2 / time_gap_sec
    c_out_q_h2o_ove = c_s_v_ove / c_p_sk
    c_s_v_ove += (overland_flow / 1000.0 * area_m2 - c_out_q_h2o_ove *
        time_gap_sec)
    if c_s_v_ove < 0.0:
        c_s_v_ove = 0.0
    c_out_q_h2o_dra = c_s_v_dra / c_p_sk
    c_s_v_dra += drain_flow / 1000.0 * area_m2 - c_out_q_h2o_dra * time_gap_sec
    if c_s_v_dra < 0.0:
        c_s_v_dra = 0.0
    c_out_q_h2o_int = c_s_v_int / c_p_fk
    c_s_v_int += inter_flow / 1000.0 * area_m2 - c_out_q_h2o_int * time_gap_sec
    if c_s_v_int < 0.0:
        c_s_v_int = 0.0
    c_out_q_h2o_sgw = c_s_v_sgw / c_p_gk
    c_s_v_sgw += (shallow_flow / 1000.0 * area_m2 - c_out_q_h2o_sgw *
        time_gap_sec)
    if c_s_v_sgw < 0.0:
        c_s_v_sgw = 0.0
    c_out_q_h2o_dgw = c_s_v_dgw / c_p_gk
    c_s_v_dgw += deep_flow / 1000.0 * area_m2 - c_out_q_h2o_dgw * time_gap_sec
    if c_s_v_dgw < 0.0:
        c_s_v_dgw = 0.0
    return (c_out_aeva, c_out_q_h2o_ove, c_out_q_h2o_dra, c_out_q_h2o_int,
        c_out_q_h2o_sgw, c_out_q_h2o_dgw, c_s_v_ove, c_s_v_dra, c_s_v_int,
        c_s_v_sgw, c_s_v_dgw, list_lvl_lyr[1] / 1000.0 * area_m2, 
        list_lvl_lyr[2] / 1000.0 * area_m2, list_lvl_lyr[3] / 1000.0 *
        area_m2, list_lvl_lyr[4] / 1000.0 * area_m2, list_lvl_lyr[5] / 
        1000.0 * area_m2, list_lvl_lyr[6] / 1000.0 * area_m2)
