take absolute values of dz in "stability_constraint_advection"
In my case the dz had negative values resulting in a wrong dt timestep in the "stability_constraint_advection".
current code = # dt of constituents (d) dt_x = R / (abs_v_x / top_bot.dx) dt_y = R / (abs_v_y / np.abs(top_bot.dy)) dt_z = R / (abs_v_z / dz_m)
suggested = dt_x = R / (abs_v_x / top_bot.dx) dt_y = R / (abs_v_y / np.abs(top_bot.dy)) dt_z = R / (abs_v_z / np.abs(dz_m))