def refold_vdop(vdop_art, v_nyq_vel, rth_position):
    """
    Refold Doppler velocity from PHIDP folding position.

    Parameters:
    ===========
        vdop_art: array
            Doppler velocity
        v_nyq_vel: float
            Nyquist velocity.
        rth_position: list
            Folding position of PHIDP along axis 2 (length of rth_position is
            length of axis 1).

    Returns:
    ========
        new_vdop: array
            Properly folded doppler velocity.
    """
    new_vdop = vdop_art
    for j in range(len(rth_position)):
        i = rth_position[j]
        if i == 0:
            continue
        else:
            new_vdop[j, i:] += v_nyq_vel

    pos = (vdop_art > v_nyq_vel)
    new_vdop[pos] = new_vdop[pos] - 2*v_nyq_vel

    return new_vdop