def apply_reflective(U1, U2):
        U1[0]  = U1[1]
        U1[-1] = U1[-2]
        U2[0]  = -U2[1]
        U2[-1] = -U2[-2]
        return U1, U2
    # Storage
    h_store = []
    v_store = []
    t_store = []
    plt.ion()
    pbar = tqdm()

    while t < runningtime: # Time loop
        U1, U2 = apply_reflective(U1, U2) # enforce reflective BC every step 
        h = U1
        v = U2 / U1
        dt = CFL * dx / np.max(np.abs(v) + np.sqrt(g * h))
        t += dt
        c = dt / dx
        pbar.set_description(f"t={t:.2f}")
        pbar.update(1)
        # WELL-BALANCED FLUX
        F1 = np.zeros(N+1)
        F2L = np.zeros(N+1)
        F2R = np.zeros(N+1)
        for i in range(0, N+1):
            # Left state for interface 'i'
            hL, huL, bL = U1[i], U2[i], b[i]
            # Right state for interface 'i'
            hR, huR, bR = U1[i+1], U2[i+1], b[i+1]
            # Free surface elevations
            etaL = hL + bL
            etaR = hR + bR
            # Interface bed elevation for hydrostatic reconstruction
            b_int = max(bL, bR)
            # Reconstructed depths at the interface (h_star)
            hL_star = max(0.0, etaL - b_int)
            hR_star = max(0.0, etaR - b_int)
            # Velocities based on original momentum and depth
            uL = huL / hL if hL > 1e-12 else 0.0
            uR = huR / hR if hR > 1e-12 else 0.0
            # Reconstructed momentum at the interface
            huL_star = hL_star * uL
            huR_star = hR_star * uR
            # Wave speeds
            cL = np.sqrt(g*hL_star)
            cR = np.sqrt(g*hR_star)
            a = max(abs(uL)+cL, abs(uR)+cR)
            # Homogeneous part of the fluxes
            f1L_hom = huL_star
            f2L_hom = huL_star*uL + 0.5*g*hL_star**2
            f1R_hom = huR_star
            f2R_hom = huR_star*uR + 0.5*g*hR_star**2
            f1 = 0.5*(f1L_hom + f1R_hom) - 0.5*a*(hR_star - hL_star)
            f2 = 0.5*(f2L_hom + f2R_hom) - 0.5*a*(huR_star - huL_star)
            # Pressure corrections for well-balanced property
            corrL = 0.5 * g * (hL**2 - hL_star**2) 
            corrR = 0.5 * g * (hR**2 - hR_star**2)
            F1[i] = f1
            F2L[i] = f2 + corrL # Pressure correction
            F2R[i] = f2 + corrR