@staticmethod
def _run(waterbody, datetime_time_step, logger, time_gap_sec, r_in_temp,
    r_in_c_no3, r_in_c_nh4, r_in_c_dph, r_in_c_pph, r_in_c_sed, r_p_att_no3,
    r_p_att_nh4, r_p_att_dph, r_p_att_pph, r_p_att_sed, r_s_m_no3,
    r_s_m_nh4, r_s_m_dph, r_s_m_pph, r_s_m_sed, r_cst_c_dn, r_cst_c_ni,
    r_cst_flow_tolerance, r_cst_vol_tolerance, r_in_q_h2o, r_s_v_h2o_old,
    r_s_v_h2o, r_out_q_h2o):
    """
        This function was written by Thibault Hallouin but is largely inspired by the work of Eva Mockler and
        Michael Bruen, namely for the work published in: Mockler, E., Bruen, M., Desta, M., Misstear, B., Environmental
        Protection Agency, Ireland, 2014. Pathways Project Final Report Volume 4: Catchment Modelling Tool
        (STRIVE Report). Dublin.

        The INCA model was developed by researchers at the Department of Geography and Environmental Science of the
        University of Reading, United Kingdom: http://www.reading.ac.uk/geographyandenvironmentalscience/research/INCA/

        Wade, A.J., Durand, P., Beaujouan, V., Wessel, W.W., Raat, K.J., Whitehead, P.G., Butterfield, D., Rankinen, K.
        and Lepisto, A. (2002) A nitrogen model for European catchments: INCA, new model structure and equations,
        Hydrology and Earth System Sciences, 6, 559-582.

        Wade, A. J., Whitehead, P. G. and Butterfield, D. (2002) The Integrated Catchments model of Phosphorus dynamics
        (INCA-P), a new approach for multiple source assessment in heterogeneous river systems: model structure and
        equations, Hydrology and Earth Systems Sciences, 6, 583-606.

        River Constants
        _ time_gap_sec        time gap between two simulation time steps [seconds]

        River Model * r_ *
        _ Water Quality
        ___ Inputs * in_ *
        _____ r_in_temp       water temperature [degree celsius]
        _____ r_in_c_no3      nitrate concentration at inlet [kg/m3]
        _____ r_in_c_nh4      ammonia concentration at inlet [kg/m3]
        _____ r_in_c_dph      dissolved phosphorus concentration at inlet [kg/m3]
        _____ r_in_c_pph      particulate phosphorus concentration at inlet [kg/m3]
        _____ r_in_c_sed      sediment concentration at inlet [kg/m3]
        ___ Parameters * p_ *
        _____ r_p_att_no3     daily attenuation factor for nitrate [-]
        _____ r_p_att_nh4     daily attenuation factor for ammonia [-]
        _____ r_p_att_dph     daily attenuation factor for dissolved phosphorus [-]
        _____ r_p_att_pph     daily attenuation factor for particulate phosphorus [-]
        _____ r_p_att_sed     daily attenuation factor for sediment [-]
        ___ States * s_ *
        _____ r_s_m_no3       quantity of nitrate in store [kg]
        _____ r_s_m_nh4       quantity of ammonia in store [kg]
        _____ r_s_m_dph       quantity of dissolved phosphorus in store [kg]
        _____ r_s_m_pph       quantity of particulate phosphorus in store [kg]
        _____ r_s_m_sed       quantity of sediment in store [kg]
        ___ Constants * cst_ *
        _____ r_cst_c_dn      denitrification rate constant [-]
        _____ r_cst_c_ni      nitrification rate constant [-]
        _____ r_cst_flow_tolerance      minimum flow to consider equal to no flow [m3/s]
        _____ r_cst_vol_tolerance       minimum volume to consider equal to empty [m3]
        ___ Outputs * out_ *
        _____ r_out_c_no3     nitrate concentration at outlet [kg/m3]
        _____ r_out_c_nh4     ammonia concentration at outlet [kg/m3]
        _____ r_out_c_dph     dissolved phosphorus concentration at outlet [kg/m3]
        _____ r_out_c_pph     particulate phosphorus concentration at outlet [kg/m3]
        _____ r_out_c_sed     sediment concentration at outlet [kg/m3]
        """
    if r_in_q_h2o < r_cst_flow_tolerance:
        logger.debug(''.join(['INCAS # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Inflow to River Store too low, inflow concentrations set to zero.'
            ]))
        r_in_c_no3 = 0.0
        r_in_c_nh4 = 0.0
        r_in_c_dph = 0.0
        r_in_c_pph = 0.0
        r_in_c_sed = 0.0
    if r_s_v_h2o_old < r_cst_vol_tolerance:
        logger.debug(''.join(['INCAS # ', waterbody, ': ',
            datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
            ' - Volume in River Store too low, in-store contaminant quantities and outflow concentrations set to zero.'
            ]))
        r_s_m_no3 = 0.0
        r_s_m_nh4 = 0.0
        r_s_m_dph = 0.0
        r_s_m_pph = 0.0
        r_s_m_sed = 0.0
        r_out_c_no3 = 0.0
        r_out_c_nh4 = 0.0
        r_out_c_dph = 0.0
        r_out_c_pph = 0.0
        r_out_c_sed = 0.0
    else:
        r_s_m_no3_old = r_s_m_no3
        concentration_no3 = r_s_m_no3_old / r_s_v_h2o_old
        if r_in_temp < 0.0:
            c10 = 0.0
            c11 = 0.0
        else:
            c10 = r_cst_c_ni * 1.047 ** (r_in_temp - 20.0)
            c11 = r_cst_c_dn * 1.047 ** (r_in_temp - 20.0)
            if c10 < 0.0:
                c10 = 0.0
            elif c10 > 1.0:
                c10 = 1.0
            if c11 < 0.0:
                c11 = 0.0
            elif c11 > 1.0:
                c11 = 1.0
        rni = c10 * r_s_m_nh4
        rdn = c11 * r_s_m_no3
        r_s_m_no3 = r_s_m_no3_old + rni - rdn + (r_in_c_no3 * r_in_q_h2o - 
            concentration_no3 * r_out_q_h2o) * time_gap_sec
        if r_s_m_no3 < 0.0:
            logger.debug(''.join(['INCAS # ', waterbody, ': ',
                datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
                ' - NO3 Quantity has gone negative in River Store, quantity reset to zero.'
                ]))
            r_s_m_no3 = 0.0
        if (r_s_v_h2o > r_cst_vol_tolerance and r_out_q_h2o >
            r_cst_flow_tolerance):
            r_out_c_no3 = r_s_m_no3 / r_s_v_h2o
        else:
            logger.debug(''.join(['INCAS # ', waterbody, ': ',
                datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
                ' - Volume/Flow in River Store too low, outflow NO3 concentration set to zero.'
                ]))
            r_out_c_no3 = 0.0
        r_s_m_nh4_old = r_s_m_nh4
        concentration_nh4 = r_s_m_nh4_old / r_s_v_h2o_old
        r_s_m_nh4 = r_s_m_nh4_old - rni + (r_in_c_nh4 * r_in_q_h2o - 
            concentration_nh4 * r_out_q_h2o) * time_gap_sec
        if r_s_m_nh4 < 0.0:
            r_s_m_nh4 = 0.0
        if (r_s_v_h2o > r_cst_vol_tolerance and r_out_q_h2o >
            r_cst_flow_tolerance):
            r_out_c_nh4 = r_s_m_nh4 / r_s_v_h2o
        else:
            logger.debug(''.join(['INCAS # ', waterbody, ': ',
                datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
                ' - Volume/Flow in River Store too low, outflow NH4 concentration set to zero.'
                ]))
            r_out_c_nh4 = 0.0
        r_s_m_dph_old = r_s_m_dph
        concentration_dph = r_s_m_dph_old / r_s_v_h2o_old
        r_s_m_dph = r_s_m_dph_old + (r_in_c_dph * r_in_q_h2o - 
            concentration_dph * r_out_q_h2o) * time_gap_sec
        if r_s_m_dph < 0.0:
            r_s_m_dph = 0.0
        r_s_m_dph = r_s_m_dph * r_p_att_dph
        if (r_s_v_h2o > r_cst_vol_tolerance and r_out_q_h2o >
            r_cst_flow_tolerance):
            r_out_c_dph = r_s_m_dph / r_s_v_h2o
        else:
            logger.debug(''.join(['INCAS # ', waterbody, ': ',
                datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
                ' - Volume/Flow in River Store too low, outflow DPH concentration set to zero.'
                ]))
            r_out_c_dph = 0.0
        r_s_m_pph_old = r_s_m_pph
        concentration_pph = r_s_m_pph_old / r_s_v_h2o_old
        r_s_m_pph = r_s_m_pph_old + (r_in_c_pph * r_in_q_h2o - 
            concentration_pph * r_out_q_h2o) * time_gap_sec
        if r_s_m_pph < 0.0:
            r_s_m_pph = 0.0
        r_s_m_pph = r_s_m_pph * r_p_att_pph
        if (r_s_v_h2o > r_cst_vol_tolerance and r_out_q_h2o >
            r_cst_flow_tolerance):
            r_out_c_pph = r_s_m_pph / r_s_v_h2o
        else:
            logger.debug(''.join(['INCAS # ', waterbody, ': ',
                datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
                ' - Volume/Flow in River Store too low, outflow PPH concentration set to zero.'
                ]))
            r_out_c_pph = 0.0
        r_s_m_sed_old = r_s_m_sed
        concentration_sed = r_s_m_sed_old / r_s_v_h2o_old
        r_s_m_sed = r_s_m_sed_old + (r_in_c_sed * r_in_q_h2o - 
            concentration_sed * r_out_q_h2o) * time_gap_sec
        if r_s_m_sed < 0.0:
            r_s_m_sed = 0.0
        r_s_m_sed = r_s_m_sed * r_p_att_sed
        if (r_s_v_h2o > r_cst_vol_tolerance and r_out_q_h2o >
            r_cst_flow_tolerance):
            r_out_c_sed = r_s_m_sed / r_s_v_h2o
        else:
            logger.debug(''.join(['INCAS # ', waterbody, ': ',
                datetime_time_step.strftime('%d/%m/%Y %H:%M:%S'),
                ' - Volume/Flow in River Store too low, outflow SED concentration set to zero.'
                ]))
            r_out_c_sed = 0.0
    return (r_out_c_no3, r_out_c_nh4, r_out_c_dph, r_out_c_pph, r_out_c_sed,
        r_s_m_no3, r_s_m_nh4, r_s_m_dph, r_s_m_pph, r_s_m_sed)
