diff --git a/examples/weather_v2.cpp b/examples/weather_v2.cpp
index b6409b30..d513ab93 100644
--- a/examples/weather_v2.cpp
+++ b/examples/weather_v2.cpp
@@ -44,6 +44,20 @@ auto make_perlin(random::generator & rng, int min_octave, int max_octave, float
 	return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights));
 }
 
+void add_random_scalar_field(random::generator & rng, util::ndarray<float, 2> & result, int min_octave, int max_octave, float scale)
+{
+	auto noise = make_perlin(rng, min_octave, max_octave, 2.f);
+
+	for (int y = 0; y < result.height(); ++y)
+	{
+		for (int x = 0; x < result.width(); ++x)
+		{
+			math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
+			result(x, y) += noise(p) * scale;
+		}
+	}
+}
+
 void add_random_vector_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, int min_octave, int max_octave, float scale)
 {
 	auto noise_1 = make_perlin(rng, min_octave, max_octave, 2.f);
@@ -59,6 +73,21 @@ void add_random_vector_field(random::generator & rng, util::ndarray<math::vector
 	}
 }
 
+template <typename T>
+void mult_random_scalar_field(random::generator & rng, util::ndarray<T, 2> & result, int min_octave, int max_octave, float min, float max)
+{
+	auto noise = make_perlin(rng, min_octave, max_octave, 2.f);
+
+	for (int y = 0; y < result.height(); ++y)
+	{
+		for (int x = 0; x < result.width(); ++x)
+		{
+			math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
+			result(x, y) *= math::lerp(min, max, noise(p));
+		}
+	}
+}
+
 template <typename T>
 void add(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src)
 {
@@ -83,7 +112,11 @@ float do_lerp(float x, float y, float t)
 
 math::vector<float, 2> do_lerp(math::vector<float, 2> const & x, math::vector<float, 2> const & y, float t)
 {
-	return math::lerp(length(x), length(y), t) * math::slerp(math::normalized(x), math::normalized(y), t);
+	auto lx = length(x);
+	auto ly = length(y);
+	if (lx > 0.f && ly > 0.f)
+		return math::lerp(lx, ly, t) * math::slerp(x / lx, y / ly, t);
+	return math::lerp(x, y, t);
 }
 
 template <typename T>
@@ -160,7 +193,7 @@ struct solver
 	climate_snapshot & snapshot;
 
 	int solve_velocity_iterations = 0;
-	int max_solve_velocity_iterations = N / 2;
+	int max_solve_velocity_iterations = 1;
 
 	int solve_temperature_iterations = 0;
 	int max_solve_temperature_iterations = N / 4;
@@ -169,6 +202,7 @@ struct solver
 	int max_solve_humidity_iterations = 3 * N;
 
 	util::ndarray<math::vector<float, 2>, 2> initial_random_wind_velocity;
+	util::ndarray<float, 2> initial_random_wind_velocity_potential;
 	util::ndarray<float, 2> pressure;
 	util::ndarray<float, 2> new_temperature;
 	util::ndarray<float, 2> new_humidity;
@@ -177,14 +211,14 @@ struct solver
 	bool finished = false;
 
 	solver(util::ndarray<float, 2> const & terrain,	random::generator & rng, int season, int day_night, climate_snapshot & snapshot,
-		util::ndarray<math::vector<float, 2>, 2> const & base_wind_velocity)
+		util::ndarray<float, 2> const & base_wind_velocity_potential)
 		: terrain(terrain)
 		, rng(rng)
 		, season(season)
 		, day_night(day_night)
 		, snapshot(snapshot)
 	{
-		initial_random_wind_velocity.resize({N, N});
+		initial_random_wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
 		snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
 		pressure.resize({N, N}, 0.f);
 		snapshot.temperature.resize({N, N}, 0.f);
@@ -192,8 +226,23 @@ struct solver
 		snapshot.humidity.resize({N, N}, 0.f);
 		new_humidity.resize({N, N}, 0.f);
 
-		initial_random_wind_velocity = base_wind_velocity.copy();
-		add_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.002f);
+		initial_random_wind_velocity_potential = base_wind_velocity_potential.copy();
+		add_random_scalar_field(rng, initial_random_wind_velocity_potential, 2, 5, 0.4f);
+		// add_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.002f);
+
+		for (int y = 1; y + 1 < N; ++y)
+			for (int x = 1; x + 1 < N; ++x)
+			{
+				math::vector potential_gradient
+				{
+					(initial_random_wind_velocity_potential(x + 1, y) - initial_random_wind_velocity_potential(x - 1, y)) * 0.5f,
+					(initial_random_wind_velocity_potential(x, y + 1) - initial_random_wind_velocity_potential(x, y - 1)) * 0.5f,
+				};
+
+				initial_random_wind_velocity(x, y) = math::ort(potential_gradient) * 0.5f;
+			}
+
+		// mult_random_scalar_field(rng, initial_random_wind_velocity, 2, 5, 0.01f, 1.f);
 
 		for (int y = 0; y < N; ++y)
 			for (int x = 0; x < N; ++x)
@@ -220,7 +269,7 @@ struct solver
 	float expected_humidity(int x, int y)
 	{
 		if (terrain(x, y) < 0.f)
-			return 1.f;
+			return 1.5f;
 		else
 			return 0.f;
 	}
@@ -229,37 +278,37 @@ struct solver
 	{
 		if (stage == 0)
 		{
-			for (int y = 1; y + 1 < N; ++y)
-			{
-				for (int x = 1; x + 1 < N; ++x)
-				{
-					float divergence = 0.f
-						+ (initial_random_wind_velocity(x + 1, y)[0] - initial_random_wind_velocity(x - 1, y)[0]) * 0.5f
-						+ (initial_random_wind_velocity(x, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f
-						;
+			// for (int y = 1; y + 1 < N; ++y)
+			// {
+			// 	for (int x = 1; x + 1 < N; ++x)
+			// 	{
+			// 		float divergence = 0.f
+			// 			+ (initial_random_wind_velocity(x + 1, y)[0] - initial_random_wind_velocity(x - 1, y)[0]) * 0.5f
+			// 			+ (initial_random_wind_velocity(x, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f
+			// 			;
 
-					// (sum - 4.f * p(x,y)) = divergence
-					// sum - 4.f * p = divergence
-					// p = (sum - divergence) / 4.f
-					pressure(x, y) = (pressure(x - 1, y) + pressure(x + 1, y) + pressure(x, y - 1) + pressure(x, y + 1) - divergence) / 4.f;
-				}
-			}
+			// 		// (sum - 4.f * p(x,y)) = divergence
+			// 		// sum - 4.f * p = divergence
+			// 		// p = (sum - divergence) / 4.f
+			// 		pressure(x, y) = (pressure(x - 1, y) + pressure(x + 1, y) + pressure(x, y - 1) + pressure(x, y + 1) - divergence) / 4.f;
+			// 	}
+			// }
 
-			for (int y = 1; y + 1 < N; ++y)
-			{
-				for (int x = 1; x + 1 < N; ++x)
-				{
-					math::vector pressure_gradient
-					{
-						(pressure(x + 1, y) - pressure(x - 1, y)) * 0.5f,
-						(pressure(x, y + 1) - pressure(x, y - 1)) * 0.5f,
-					};
-					snapshot.wind_velocity(x, y) = initial_random_wind_velocity(x, y) - pressure_gradient;
+			// for (int y = 1; y + 1 < N; ++y)
+			// {
+			// 	for (int x = 1; x + 1 < N; ++x)
+			// 	{
+			// 		math::vector pressure_gradient
+			// 		{
+			// 			(pressure(x + 1, y) - pressure(x - 1, y)) * 0.5f,
+			// 			(pressure(x, y + 1) - pressure(x, y - 1)) * 0.5f,
+			// 		};
+			// 		snapshot.wind_velocity(x, y) = initial_random_wind_velocity(x, y) - pressure_gradient;
 
-					if (terrain(x, y) > 0.f)
-						snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
-				}
-			}
+			// 		if (terrain(x, y) > 0.f)
+			// 			snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
+			// 	}
+			// }
 
 			++solve_velocity_iterations;
 			if (solve_velocity_iterations == max_solve_velocity_iterations)
@@ -319,10 +368,17 @@ struct solver
 					math::point p{x + 0.5f, y + 0.5f};
 					float advection = sample_bilinear(snapshot.humidity, p - snapshot.wind_velocity(x, y) * step_time);
 
-					float evaporation_speed = 0.0001f;
-					float evaporation_delta = (expected_humidity(x, y) - snapshot.humidity(x, y)) * (- std::expm1(- evaporation_speed * step_time));
+					float income = 0.f;
+					if (terrain(x, y) < 0.f)
+						income = step_time * 0.00012f;
+					else
+						// income = snapshot.humidity(x, y) * std::expm1(- 0.00006f * step_time);
+						income = -std::min(snapshot.humidity(x, y), step_time * 0.00001f);
 
-					new_humidity(x, y) = evaporation_delta + advection;
+					// float evaporation_speed = 0.0001f;
+					// float evaporation_delta = (expected_humidity(x, y) - snapshot.humidity(x, y)) * (- std::expm1(- evaporation_speed * step_time));
+
+					new_humidity(x, y) = income + advection;
 				}
 			}
 			std::swap(snapshot.humidity, new_humidity);
@@ -386,6 +442,7 @@ struct weather_app
 	bool show_humidity = false;
 	bool show_biomes = true;
 	bool show_rivers = true;
+	bool show_clouds = true;
 
 	gfx::pixmap_rgba biomes_map;
 	util::ndarray<float, 2> terrain;
@@ -401,17 +458,27 @@ struct weather_app
 	int day_night = 0;
 
 	util::ndarray<math::vector<float, 2>, 2> base_wind_velocity;
+	util::ndarray<float, 2> base_wind_velocity_potential;
 
 	climate_snapshot snapshots[4][2];
 	climate_snapshot average;
 	climate_snapshot display_snapshot;
 
-	bool need_update_display_snapshot = false;
-
 	std::optional<struct solver> solver;
 
 	river_net rivers;
 
+	util::ndarray<float, 2> cloud_density;
+	util::ndarray<float, 2> new_cloud_density;
+
+	struct cloud_particle
+	{
+		math::point<float, 2> position;
+		int lifetime;
+	};
+
+	std::vector<cloud_particle> cloud_particles;
+
 	float display_season = 0.f;
 
 	util::clock<> frame_clock;
@@ -428,11 +495,16 @@ struct weather_app
 
 		biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});
 
-		base_wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
+		// base_wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
+		// add_random_vector_field(rng, base_wind_velocity, 3, 6, 0.003f);
 
-		add_random_vector_field(rng, base_wind_velocity, 3, 6, 0.003f);
+		base_wind_velocity_potential.resize({N, N}, 0.f);
+		add_random_scalar_field(rng, base_wind_velocity_potential, 1, 3, 0.6f);
 
-		context.vsync(false);
+		cloud_density.resize({N, N}, 0.f);
+		new_cloud_density.resize({N, N}, 0.f);
+
+		context.vsync(true);
 	}
 
 	void on_event(app::key_event const & event) override
@@ -456,12 +528,17 @@ struct weather_app
 
 		if (event.down && event.key == app::keycode::R)
 			show_rivers ^= true;
+
+		if (event.down && event.key == app::keycode::C)
+			show_clouds ^= true;
 	}
 
 	void update() override
 	{
 		float const dt = frame_clock.restart().count();
 
+		bool need_update_display_snapshot = false;
+
 		if (solver && solver->finished)
 		{
 			solver.reset();
@@ -480,34 +557,6 @@ struct weather_app
 			}
 		}
 
-		if (!solver && season < 4)
-			solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night], base_wind_velocity);
-
-		if (!paused)
-		{
-			if (solver)
-			{
-				for (int i = 0; i < 64; ++i)
-				{
-					solver->step();
-					if (solver->finished)
-						break;
-				}
-			}
-			else if (season == 4 && !rivers.queue.empty())
-			{
-				for (int i = 0; i < 64; ++i)
-				{
-					propagate_rivers();
-					if (rivers.queue.empty())
-					{
-						compute_river_flow();
-						break;
-					}
-				}
-			}
-		}
-
 		if (state().key_down.contains(app::keycode::LEFT))
 		{
 			display_season -= 0.5f * dt;
@@ -524,6 +573,144 @@ struct weather_app
 			display_season += 1.f;
 		while (display_season >= 1.f)
 			display_season -= 1.f;
+
+		if (season == 4 && need_update_display_snapshot)
+		{
+			update_display_snapshot();
+			need_update_display_snapshot = false;
+		}
+
+		if (!solver && season < 4)
+			solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night], base_wind_velocity_potential);
+
+		if (!paused)
+		{
+			if (solver)
+			{
+				for (int i = 0; i < 64; ++i)
+				// for (int i = 0; i < 1; ++i)
+				{
+					solver->step();
+					if (solver->finished)
+						break;
+				}
+			}
+			else if (season == 4 && !rivers.queue.empty())
+			{
+				for (int i = 0; i < 64; ++i)
+				{
+					propagate_rivers();
+					if (rivers.queue.empty())
+					{
+						compute_river_flow();
+						paused = true;
+						break;
+					}
+				}
+			}
+			else
+			{
+				float const dt = 10.f;
+
+				// if (false)
+				// for (int y = 0; y < N; ++y)
+				// {
+				// 	for (int x = 0; x < N; ++x)
+				// 	{
+				// 		new_cloud_density(x, y) = sample_bilinear(cloud_density, math::point{x + 0.5f, y + 0.5f} - display_snapshot.wind_velocity(x, y) * dt);
+				// 	}
+				// }
+
+				// for (int y = 0; y < N; ++y)
+				// {
+				// 	for (int x = 0; x < N; ++x)
+				// 	{
+				// 		new_cloud_density(x, y) += (terrain(x, y) < 0.f ? 0.0001f : 0.f) * dt;
+				// 	}
+				// }
+
+				// // if (false)
+				// for (int y = 1; y + 1 < N; ++y)
+				// {
+				// 	for (int x = 1; x + 1 < N; ++x)
+				// 	{
+				// 		auto gradient = math::vector {
+				// 			(new_cloud_density(x + 1, y) - new_cloud_density(x - 1, y)) / 2.f,
+				// 			(new_cloud_density(x, y + 1) - new_cloud_density(x, y - 1)) / 2.f,
+				// 		};
+
+				// 		if (math::length(gradient) > 0.f)
+				// 			cloud_density(x, y) = sample_bilinear(new_cloud_density, math::point{x + 0.5f, y + 0.5f} - math::normalized(gradient) * 1.f);
+				// 		else
+				// 			cloud_density(x, y) = new_cloud_density(x, y);
+				// 	}
+				// }
+
+				// std::swap(new_cloud_density, cloud_density);
+
+				int total_cloud_particles = N * N / 8;
+
+				auto new_cloud_particle = [&]
+				{
+					while (true)
+					{
+						cloud_particle p{{random::uniform<float>(rng, 0.f, N), random::uniform<float>(rng, 0.f, N)}, 0};
+						if (sample_bilinear(terrain, p.position) < 0.f)
+							return p;
+					}
+				};
+
+				if (cloud_particles.size() < total_cloud_particles)
+					for (int i = 0; i < 16; ++i)
+						cloud_particles.push_back(new_cloud_particle());
+
+				util::ndarray<std::vector<cloud_particle *>, 2> cloud_particle_index({N, N});
+
+				for (auto & p : cloud_particles)
+				{
+					p.position += sample_bilinear(display_snapshot.wind_velocity, p.position) * dt;
+
+					p.lifetime += 1;
+
+					auto r = math::distance_sqr(p.position, math::point{N/2.f, N/2.f}) / math::sqr(N/2.f);
+					auto t = std::min(1.f, 2.f * std::exp(- 1.f * r));
+
+					if (p.lifetime >= 1024 || p.position[0] < 0.f || p.position[0] > N || p.position[1] < 0.f || p.position[1] > N || (p.lifetime > 16 && random::uniform<float>(rng) > t))
+					{
+						p = new_cloud_particle();
+					}
+
+					cloud_particle_index(std::min<int>(N-1, std::floor(p.position[0])), std::min<int>(N-1, std::floor(p.position[1]))).push_back(&p);
+				}
+
+				for (auto & p : cloud_particles)
+				{
+					int x = std::min<int>(N-1, std::floor(p.position[0]));
+					int y = std::min<int>(N-1, std::floor(p.position[1]));
+
+					for (int tx = -1; tx <= 1; ++tx)
+					{
+						for (int ty = -1; ty <= 1; ++ty)
+						{
+							int xx = x + tx;
+							int yy = y + ty;
+							if (xx < 0 || xx + 1 >= N || yy < 0 || yy + 1 >= N) continue;
+
+							for (auto q : cloud_particle_index(xx, yy))
+							{
+								if (q == &p) continue;
+								auto r = p.position - q->position;
+								if (auto l = math::length(r); l < 1.f)
+								{
+									auto n = r / l;
+									p.position += n * dt * 0.02f;
+								}
+							}
+						}
+					}
+				}
+			}
+		}
 	}
 
 	void compute_average()
@@ -679,12 +866,6 @@ struct weather_app
 
 	void present() override
 	{
-		if (need_update_display_snapshot)
-		{
-			update_display_snapshot();
-			need_update_display_snapshot = false;
-		}
-
 		gl::ClearColor(0.f, 0.f, 0.f, 0.f);
 		gl::Clear(gl::COLOR_BUFFER_BIT);
 
@@ -780,6 +961,11 @@ struct weather_app
 					color = gfx::blend(color, map_humidity(snapshot.humidity(x, y)));
 				}
 
+				// if (show_clouds)
+				// {
+				// 	color = gfx::blend(color, map_color_positive(cloud_density(x, y), {1.f, 1.f, 1.f, 0.f}, {1.f, 1.f, 1.f, 1.f}));
+				// }
+
 				painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color));
 			}
 		}
@@ -839,6 +1025,12 @@ struct weather_app
 				}
 		}
 
+		if (show_clouds)
+		{
+			for (auto const & p : cloud_particles)
+				painter.circle(p.position, 2.f, {255, 255, 255, std::min(255, p.lifetime)}, {255, 255, 255, 0}, 12);
+		}
+
 		int row = 0;
 		auto push_text = [&](std::string const & text) mutable
 		{