def generate_single_return_period(args):
    """
    This function calculates a single return period for a single reach
    """
    qout_file, return_period_file, rivid_index_list, step, num_years, \
        method, mp_lock = args

    skewvals = [-3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2,
                -1.0, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0,
                1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0]
    kfac2 = [0.396, 0.384, 0.368, 0.351, 0.33, 0.307, 0.282, 0.254, 0.225,
             0.195, 0.164, 0.132, 0.099, 0.066, 0.033, 0, -0.033, -0.066,
             -0.099, -0.132, -0.164, -0.195, -0.225, -0.254, -0.282, -0.307,
             -0.33, -0.351, -0.368, -0.384, -0.396]
    kfac10 = [0.66, 0.702, 0.747, 0.795, 0.844, 0.895, 0.945, 0.994, 1.041,
              1.086, 1.128, 1.166, 1.2, 1.231, 1.258, 1.282, 1.301, 1.317,
              1.328, 1.336, 1.34, 1.34, 1.337, 1.329, 1.318, 1.302, 1.284,
              1.262, 1.238, 1.21, 1.18]
    kfac25 = [.666, .712, .764, .823, .888, .959, 1.035, 1.116, 1.198, 1.282,
              1.366, 1.448, 1.528, 1.606, 1.680, 1.751, 1.818, 1.880, 1.939,
              1.993, 2.043, 2.087, 2.128, 2.163, 2.193, 2.219, 2.240, 2.256,
              2.267, 2.275, 2.278]
    kfac50 = [0.666, 0.714, 0.768, 0.83, 0.9, 0.98, 1.069, 1.166, 1.27, 1.379,
              1.492, 1.606, 1.72, 1.834, 1.945, 2.054, 2.159, 2.261, 2.359,
              2.453, 2.542, 2.626, 2.706, 2.78, 2.848, 2.912, 2.97, 3.023,
              3.071, 3.114, 3.152]
    kfac100 = [0.667, 0.714, 0.769, 0.832, 0.905, 0.99, 1.087, 1.197, 1.318,
               1.499, 1.588, 1.733, 1.88, 2.029, 2.178, 2.326, 2.472, 2.615,
               2.755, 2.891, 3.022, 3.149, 3.271, 3.388, 3.499, 3.605, 3.705,
               3.8, 3.889, 3.973, 4.051]

    with RAPIDDataset(qout_file) as qout_nc_file:
        # get index of return period data
        if method == 'weibull':
            rp_index_20 = int((num_years + 1)/20.0)
            rp_index_10 = int((num_years + 1)/10.0)
            rp_index_2 = int((num_years + 1)/2.0)

        if method == 'weibull':
            return_20_array = np.zeros(len(rivid_index_list))
        elif method == 'gumble':
            return_100_array = np.zeros(len(rivid_index_list))
            return_50_array = np.zeros(len(rivid_index_list))
            return_20_array = np.zeros(len(rivid_index_list))
        elif method == 'log_pearson':
            return_100_array = np.zeros(len(rivid_index_list))
            return_50_array = np.zeros(len(rivid_index_list))
            return_25_array = np.zeros(len(rivid_index_list))
        return_10_array = np.zeros(len(rivid_index_list))
        return_2_array = np.zeros(len(rivid_index_list))
        max_flow_array = np.zeros(len(rivid_index_list))

        # iterate through rivids to generate return periods
        for iter_idx, rivid_index in enumerate(rivid_index_list):
            filtered_flow_data = qout_nc_file.get_qout_index(
                rivid_index,
                pd_filter="{0}D".format(step),
                filter_mode="max")
            sorted_flow_data = np.sort(filtered_flow_data)[:num_years:-1]
            max_flow = sorted_flow_data[0]
            if max_flow < 0.01:
                log("Return period data < 0.01 generated for rivid {0}"
                    .format(qout_nc_file.qout_nc.variables[
                        qout_nc_file.river_id_dimension][rivid_index]),
                    "WARNING")
            max_flow_array[iter_idx] = max_flow

            if method == 'weibull':
                return_20_array[iter_idx] = sorted_flow_data[rp_index_20]
                return_10_array[iter_idx] = sorted_flow_data[rp_index_10]
                return_2_array[iter_idx] = sorted_flow_data[rp_index_2]

            elif method == 'gumble':
                mean_flow = np.mean(filtered_flow_data)
                stddev = np.std(filtered_flow_data)
                return_100_array[iter_idx] = mean_flow + 3.14*stddev
                return_50_array[iter_idx] = mean_flow + 2.59*stddev
                return_20_array[iter_idx] = mean_flow + 1.87*stddev
                return_10_array[iter_idx] = mean_flow + 1.3*stddev
                return_2_array[iter_idx] = mean_flow - .164*stddev

            elif method == 'log_pearson':
                log_flow = np.log10(filtered_flow_data[filtered_flow_data > 0])
                if len(log_flow) <= 0:
                    continue
                mean_log_flow = np.mean(log_flow)
                std_log_flow = np.std(log_flow)
                log_flow_array = np.array(log_flow)
                skew = (num_years * (np.sum(
                    np.power((log_flow_array - mean_log_flow), 3)))) / \
                    ((num_years - 1) * (num_years - 2) * std_log_flow ** 3)
                k2 = np.interp(skew, skewvals, kfac2)
                k10 = np.interp(skew, skewvals, kfac10)
                k25 = np.interp(skew, skewvals, kfac25)
                k50 = np.interp(skew, skewvals, kfac50)
                k100 = np.interp(skew, skewvals, kfac100)
                return_100_array[iter_idx] = \
                    np.power(10, (mean_log_flow + k100*std_log_flow))
                return_50_array[iter_idx] = \
                    np.power(10, (mean_log_flow + k50*std_log_flow))
                return_25_array[iter_idx] = \
                    np.power(10, (mean_log_flow + k25*std_log_flow))
                return_10_array[iter_idx] = \
                    np.power(10, (mean_log_flow + k10*std_log_flow))
                return_2_array[iter_idx] = \
                    np.power(10, (mean_log_flow + k2*std_log_flow))

        mp_lock.acquire()
        return_period_nc = Dataset(return_period_file, 'a')
        return_period_nc.variables['max_flow'][rivid_index_list] = \
            max_flow_array
        if method == 'weibull':
            return_period_nc.variables['return_period_20'][
                rivid_index_list] = return_20_array
        elif method in 'gumble':
            return_period_nc.variables['return_period_100'][
                rivid_index_list] = return_100_array
            return_period_nc.variables['return_period_50'][
                rivid_index_list] = return_50_array
            return_period_nc.variables['return_period_20'][
                rivid_index_list] = return_20_array
        elif method == 'log_pearson':
            return_period_nc.variables['return_period_100'][
                rivid_index_list] = return_100_array
            return_period_nc.variables['return_period_50'][
                rivid_index_list] = return_50_array
            return_period_nc.variables['return_period_25'][
                rivid_index_list] = return_25_array
        return_period_nc.variables['return_period_10'][
            rivid_index_list] = return_10_array
        return_period_nc.variables['return_period_2'][
            rivid_index_list] = return_2_array
        return_period_nc.close()
        mp_lock.release()
This function calculates a single return period for a single reach