Couple of notes, I think you forgot to apply timestep when adding rocket exhaust velocity, pretty sure it should be

  u[idx] += backward.x * flame_velocity_amount * falloff * delta
  v[idx] += backward.y * flame_velocity_amount * falloff * delta

I think it would help if you put the advection step in a separate function, given it's identical between the density and velocity advection. Also if you put the grid interpolation in a separate function, you can easily use the midpoint method for advection:

  func advect(i: int, j: int, grid: PackedFloat32Array, grid_prev: PackedFloat32Array, vel_u: PackedFloat32Array, vel_v: PackedFloat32Array, dt: float) -> void:
    # use midpoint method
    var idx := IX(i, j)
    # first take half-step using current velocity vector
    var xh := i - 0.5 * dt * vel_u[idx]
    var yh := j - 0.5 * dt * vel_v[idx]
    #  get interpolated velocity vector from midpoint position
    var uh := sample_grid(xh, yh, vel_u)
    var vh := sample_grid(xh, yh, vel_v)
    # use midpoint velocity vector to do a regular full step
    var x := i - dt * uh
    var y := j - dt * vh
    grid[idx] = sample_grid(x, y, grid_prev)

  for j in range(1, N + 1):
    for i in range(1, N + 1):
      advect(i, j, u, u_prev, u_prev, v_prev, dt0)
      advect(i, j, v, v_prev, u_prev, v_prev, dt0)