def flux(h, hu):
    u = np.where(h > 1e-12, hu/h, 0.0)
    f1 = hu
    f2 = hu*u + 0.5*g*h*h
    return f1, f2
def rusanov_flux(hL, huL, hR, huR):
    uL = huL/hL if hL > 1e-12 else 0.0
    uR = huR/hR if hR > 1e-12 else 0.0
    cL = np.sqrt(g*hL)
    cR = np.sqrt(g*hR)
    smax = max(abs(uL)+cL, abs(uR)+cR)
    f1L, f2L = flux(hL, huL)
    f1R, f2R = flux(hR, huR)
    f1 = 0.5*(f1L + f1R) - 0.5*smax*(hR - hL)
    f2 = 0.5*(f2L + f2R) - 0.5*smax*(huR - huL)
    return f1, f2

def step_fv(h, hu, dt):
    hg = np.zeros(nx+2)
    hug = np.zeros(nx+2)
    hg[1:-1] = h
    hug[1:-1] = hu
    hg, hug = fill_reflective_ghosts(hg, hug)
    Fh = np.zeros(nx+1)
    Fhu = np.zeros(nx+1)
    for i in range(nx+1):
        hL, huL = hg[i], hug[i]
        hR, huR = hg[i+1], hug[i+1]
        f1, f2 = rusanov_flux(hL, huL, hR, huR)
        Fh[i] = f1
        Fhu[i] = f2
    h_new = h.copy()
    hu_new = hu.copy()
    for j in range(nx):
        h_new[j]  = h[j]  - (dt/dx)*(Fh[j+1]  - Fh[j])
        hu_new[j] = hu[j] - (dt/dx)*(Fhu[j+1] - Fhu[j])
    # wind forcing
    if use_wind:
        hu_new += dt * (h_new * a_wind)
    h_new = np.maximum(h_new, 1e-8)
    return h_new, hu_new
# Run simulation
t = 0.0
while t < t_end:
    u_now = np.where(h > 1e-12, hu/h, 0.0)
    c_now = np.sqrt(g*h)
    smax = np.max(np.abs(u_now) + c_now)
    dt = CFL * dx / (smax + 1e-12)
    t += dt
    h, hu = step_fv(h, hu, dt)
    u = np.where(h > 1e-12, hu/h, 0.0)
    eta = h - H0