def GetITout(all_output_times, time_stations_out, tol=1.0e-8, with_count=False):
    """
    Get indices and output times
    ============================
    INPUT:
        all_output_times  : array with all output times. ex: [0,0.1,0.2,0.22,0.3,0.4]
        time_stations_out : time stations for output: ex: [0,0.2,0.4]  # must be sorted ascending
        tol               : tolerance to compare times
    RETURNS:
        iout : indices of times in all_output_times
        tout : times corresponding to iout
    """
    I, T        = [], []                          # indices and times for output
    lower_index = 0                               # lower index in all_output_times
    len_aotimes = len(all_output_times)           # length of a_o_times
    for t in time_stations_out:                   # for all requested output times
        if t < 0:                                 # final time
            k = len(all_output_times)-1           # last output time
            I.append(k)                           # append index to iout
            T.append(all_output_times[k])         # append last output time
            continue                              # skip search
        for k in range(lower_index, len_aotimes): # search within a_o_times
            if abs(t-all_output_times[k]) < tol:  # found near match
                lower_index += 1                  # update index
                I.append(k)                       # add index to iout
                T.append(all_output_times[k])     # add time to tout
                break                             # stop searching for this 't'
            if all_output_times[k] > t:           # failed to search for 't'
                lower_index = k                   # update idx to start from here on
                break                             # skip this 't' and try the next one
    if with_count: return zip(range(len(I)), I, T)
    else:          return zip(I, T)