@staticmethod
def _run(waterbody, datetime_time_step, logger, 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_h2o_ove, c_s_v_h2o_dra, c_s_v_h2o_int,
    c_s_v_h2o_sgw, c_s_v_h2o_dgw, c_s_v_h2o_ly1, c_s_v_h2o_ly2,
    c_s_v_h2o_ly3, c_s_v_h2o_ly4, c_s_v_h2o_ly5, c_s_v_h2o_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 Constants
        _ area_m2                   catchment area [m2]
        _ time_gap_sec              time gap between two simulation time steps [seconds]

        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_h2o_ove         volume of water in overland store [m3]
        _____ c_s_v_h2o_dra         volume of water in drain store [m3]
        _____ c_s_v_h2o_int         volume of water in inter store [m3]
        _____ c_s_v_h2o_sgw         volume of water in shallow groundwater store [m3]
        _____ c_s_v_h2o_dgw         volume of water in deep groundwater store [m3]
        _____ c_s_v_h2o_ly1         volume of water in first soil layer store [m3]
        _____ c_s_v_h2o_ly2         volume of water in second soil layer store [m3]
        _____ c_s_v_h2o_ly3         volume of water in third soil layer store [m3]
        _____ c_s_v_h2o_ly4         volume of water in fourth soil layer store [m3]
        _____ c_s_v_h2o_ly5         volume of water in fifth soil layer store [m3]
        _____ c_s_v_h2o_ly6         volume of water in sixth soil layer store [m3]

        ___ Processes * pr_ *
        _____ c_pr_eff_rain_to_ove  effective rainfall converted into overland flow runoff [mm]
        _____ c_pr_eff_rain_to_dra  effective rainfall converted into to drain flow runoff [mm]
        _____ c_pr_eff_rain_to_int  effective rainfall converted into to interflow runoff [mm]
        _____ c_pr_eff_rain_to_sgw  effective rainfall converted into to shallow groundwater flow runoff [mm]
        _____ c_pr_eff_rain_to_dgw  effective rainfall converted into to deep groundwater flow runoff [mm]
        ___ Outputs * out_ *
        _____ c_out_aeva            actual evapotranspiration [m3/s]
        _____ c_out_q_h2o_ove       overland flow [m3/s]
        _____ c_out_q_h2o_dra       drain flow [m3/s]
        _____ c_out_q_h2o_int       inter flow [m3/s]
        _____ c_out_q_h2o_sgw       shallow groundwater flow [m3/s]
        _____ c_out_q_h2o_dgw       deep groundwater flow [m3/s]
        _____ c_out_q_h2o           total outflow [m3/s]
        """
    nb_soil_layers = 6.0
    c_p_sk *= 3600.0
    c_p_fk *= 3600.0
    c_p_gk *= 3600.0
    dict_z_lyr = dict()
    for i in [1, 2, 3, 4, 5, 6]:
        dict_z_lyr[i] = c_p_z / nb_soil_layers
    dict_lvl_lyr = dict()
    dict_lvl_lyr[1] = c_s_v_h2o_ly1 / area_m2 * 1000.0
    dict_lvl_lyr[2] = c_s_v_h2o_ly2 / area_m2 * 1000.0
    dict_lvl_lyr[3] = c_s_v_h2o_ly3 / area_m2 * 1000.0
    dict_lvl_lyr[4] = c_s_v_h2o_ly4 / area_m2 * 1000.0
    dict_lvl_lyr[5] = c_s_v_h2o_ly5 / area_m2 * 1000.0
    dict_lvl_lyr[6] = c_s_v_h2o_ly6 / area_m2 * 1000.0
    lvl_total_start = 0.0
    for i in [1, 2, 3, 4, 5, 6]:
        lvl_total_start += dict_lvl_lyr[i]
    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)
        c_pr_eff_rain_to_ove = h_prime * excess_rain
        excess_rain -= c_pr_eff_rain_to_ove
        for i in [1, 2, 3, 4, 5, 6]:
            space_in_lyr = dict_z_lyr[i] - dict_lvl_lyr[i]
            if excess_rain <= space_in_lyr:
                dict_lvl_lyr[i] += excess_rain
                excess_rain = 0.0
            else:
                dict_lvl_lyr[i] = dict_z_lyr[i]
                excess_rain -= space_in_lyr
        c_pr_eff_rain_to_dra = c_p_d * excess_rain
        c_pr_eff_rain_to_int = (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 = dict_lvl_lyr[i] * s_prime ** i
            if leak_interflow < dict_lvl_lyr[i]:
                c_pr_eff_rain_to_int += leak_interflow
                dict_lvl_lyr[i] -= leak_interflow
        c_pr_eff_rain_to_sgw = 0.0
        for i in [1, 2, 3, 4, 5, 6]:
            leak_shallow_flow = dict_lvl_lyr[i] * (s_prime / i)
            if leak_shallow_flow < dict_lvl_lyr[i]:
                c_pr_eff_rain_to_sgw += leak_shallow_flow
                dict_lvl_lyr[i] -= leak_shallow_flow
        c_pr_eff_rain_to_dgw = 0.0
        for i in [6, 5, 4, 3, 2, 1]:
            leak_deep_flow = dict_lvl_lyr[i] * s_prime ** (7 - i)
            if leak_deep_flow < dict_lvl_lyr[i]:
                c_pr_eff_rain_to_dgw += leak_deep_flow
                dict_lvl_lyr[i] -= leak_deep_flow
    else:
        c_pr_eff_rain_to_ove = 0.0
        c_pr_eff_rain_to_dra = 0.0
        c_pr_eff_rain_to_int = 0.0
        c_pr_eff_rain_to_sgw = 0.0
        c_pr_eff_rain_to_dgw = 0.0
        deficit_rain = excess_rain * -1.0
        aeva += rain
        for i in [1, 2, 3, 4, 5, 6]:
            if dict_lvl_lyr[i] >= deficit_rain:
                dict_lvl_lyr[i] -= deficit_rain
                aeva += deficit_rain
                deficit_rain = 0.0
            else:
                aeva += dict_lvl_lyr[i]
                deficit_rain = c_p_c * (deficit_rain - dict_lvl_lyr[i])
                dict_lvl_lyr[i] = 0.0
    lvl_total_end = 0.0
    for i in [1, 2, 3, 4, 5, 6]:
        lvl_total_end += dict_lvl_lyr[i]
    c_out_aeva = aeva / 1000.0 * area_m2 / time_gap_sec
    c_out_q_h2o_ove = c_s_v_h2o_ove / c_p_sk
    c_s_v_h2o_ove += (c_pr_eff_rain_to_ove / 1000.0 * area_m2 - 
        c_out_q_h2o_ove * time_gap_sec)
    if c_s_v_h2o_ove < 0.0:
        logger.debug(''.join(['SMART # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Volume in OVE Store has gone negative, volume reset to zero.'])
            )
        c_s_v_h2o_ove = 0.0
    c_out_q_h2o_dra = c_s_v_h2o_dra / c_p_sk
    c_s_v_h2o_dra += (c_pr_eff_rain_to_dra / 1000.0 * area_m2 - 
        c_out_q_h2o_dra * time_gap_sec)
    if c_s_v_h2o_dra < 0.0:
        logger.debug(''.join(['SMART # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Volume in DRA Store has gone negative, volume reset to zero.'])
            )
        c_s_v_h2o_dra = 0.0
    c_out_q_h2o_int = c_s_v_h2o_int / c_p_fk
    c_s_v_h2o_int += (c_pr_eff_rain_to_int / 1000.0 * area_m2 - 
        c_out_q_h2o_int * time_gap_sec)
    if c_s_v_h2o_int < 0.0:
        logger.debug(''.join(['SMART # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Volume in INT Store has gone negative, volume reset to zero.'])
            )
        c_s_v_h2o_int = 0.0
    c_out_q_h2o_sgw = c_s_v_h2o_sgw / c_p_gk
    c_s_v_h2o_sgw += (c_pr_eff_rain_to_sgw / 1000.0 * area_m2 - 
        c_out_q_h2o_sgw * time_gap_sec)
    if c_s_v_h2o_sgw < 0.0:
        logger.debug(''.join(['SMART # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Volume in SGW Store has gone negative, volume reset to zero.'])
            )
        c_s_v_h2o_sgw = 0.0
    c_out_q_h2o_dgw = c_s_v_h2o_dgw / c_p_gk
    c_s_v_h2o_dgw += (c_pr_eff_rain_to_dgw / 1000.0 * area_m2 - 
        c_out_q_h2o_dgw * time_gap_sec)
    if c_s_v_h2o_dgw < 0.0:
        logger.debug(''.join(['SMART # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Volume in DGW Store has gone negative, volume reset to zero.'])
            )
        c_s_v_h2o_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_h2o_ove, c_s_v_h2o_dra,
        c_s_v_h2o_int, c_s_v_h2o_sgw, c_s_v_h2o_dgw, dict_lvl_lyr[1] / 
        1000.0 * area_m2, dict_lvl_lyr[2] / 1000.0 * area_m2, dict_lvl_lyr[
        3] / 1000.0 * area_m2, dict_lvl_lyr[4] / 1000.0 * area_m2, 
        dict_lvl_lyr[5] / 1000.0 * area_m2, dict_lvl_lyr[6] / 1000.0 *
        area_m2, c_pr_eff_rain_to_ove, c_pr_eff_rain_to_dra,
        c_pr_eff_rain_to_int, c_pr_eff_rain_to_sgw, c_pr_eff_rain_to_dgw)
