diff --git a/.gitignore b/.gitignore
index 4422efb240..0ddebe7c2a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,6 +33,7 @@
 *.log
 settings.ron
 run.sh
+maps
 screenshots
 todo.txt
 
diff --git a/server/src/lib.rs b/server/src/lib.rs
index 9abc4c90dd..83ba7b8d81 100644
--- a/server/src/lib.rs
+++ b/server/src/lib.rs
@@ -46,7 +46,7 @@ use std::{
 };
 use uvth::{ThreadPool, ThreadPoolBuilder};
 use vek::*;
-use world::{sim::{WORLD_SIZE, WorldOpts}, World};
+use world::{sim::{FileOpts, WORLD_SIZE, WorldOpts}, World};
 const CLIENT_TIMEOUT: f64 = 20.0; // Seconds
 
 pub enum Event {
@@ -106,6 +106,11 @@ impl Server {
 
         let world = World::generate(settings.world_seed, WorldOpts {
             seed_elements: true,
+            world_file: if let Some(ref file) = settings.map_file {
+                FileOpts::Load(file.clone())
+            } else {
+                FileOpts::Generate
+            },
             ..WorldOpts::default()
         });
 
diff --git a/server/src/settings.rs b/server/src/settings.rs
index 444ba9580a..134fc10e50 100644
--- a/server/src/settings.rs
+++ b/server/src/settings.rs
@@ -15,6 +15,7 @@ pub struct ServerSettings {
     //pub login_server: whatever
     pub start_time: f64,
     pub admins: Vec<String>,
+    pub map_file: Option<PathBuf>,
 }
 
 impl Default for ServerSettings {
@@ -27,6 +28,7 @@ impl Default for ServerSettings {
             server_description: "This is the best Veloren server.".to_owned(),
             max_players: 100,
             start_time: 9.0 * 3600.0,
+            map_file: None,
             admins: [
                 "Pfau",
                 "zesterer",
@@ -97,12 +99,13 @@ impl ServerSettings {
                 [127, 0, 0, 1],
                 pick_unused_port().expect("Failed to find unused port!"),
             )),
-            world_seed: 1337,
             server_name: "Singleplayer".to_owned(),
             server_description: "Who needs friends anyway?".to_owned(),
             max_players: 100,
             start_time: 9.0 * 3600.0,
             admins: vec!["singleplayer".to_string()], // TODO: Let the player choose if they want to use admin commands or not
+            ..
+            Self::load() // Fill in remaining fields from settings.ron.
         }
     }
 
diff --git a/world/Cargo.toml b/world/Cargo.toml
index f9e9a70756..0907d0d53d 100644
--- a/world/Cargo.toml
+++ b/world/Cargo.toml
@@ -5,6 +5,7 @@ authors = ["Joshua Barretto <joshua.s.barretto@gmail.com>"]
 edition = "2018"
 
 [dependencies]
+bincode = "1.2.0"
 common = { package = "veloren-common", path = "../common" }
 bitvec = "0.15.2"
 vek = "0.9.9"
@@ -20,6 +21,7 @@ arr_macro = "0.1.2"
 rayon = "1.2.0"
 roots = "0.0.5"
 serde = "1.0.102"
+serde_derive = "1.0.102"
 ron = "0.5.1"
 pretty_env_logger = "0.3.0"
 
diff --git a/world/examples/water.rs b/world/examples/water.rs
index 818262f183..3f960a433e 100644
--- a/world/examples/water.rs
+++ b/world/examples/water.rs
@@ -1,9 +1,9 @@
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
 // use self::Mode::*;
-use std::{f32, f64};
+use std::{f32, f64, path::PathBuf};
 use vek::*;
 use veloren_world::{
-    sim::{RiverKind, WorldOpts, WORLD_SIZE},
+    sim::{self, RiverKind, WorldOpts, WORLD_SIZE},
     util::Sampler,
     World, CONFIG,
 };
@@ -29,8 +29,17 @@ const H: usize = 1024;
 fn main() {
     pretty_env_logger::init();
 
+    let map_file =
+        // "map_1575990726223.bin";
+        // "map_1575987666972.bin";
+        "map_1576046079066.bin";
+    let mut _map_file = PathBuf::from("./maps");
+    _map_file.push(map_file);
+
     let world = World::generate(1337, WorldOpts {
         seed_elements: false,
+        // world_file: sim::FileOpts::Load(_map_file),
+        world_file: sim::FileOpts::Save,
         ..WorldOpts::default()
     });
 
@@ -102,7 +111,7 @@ fn main() {
                 let pos = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32);
                 let downhill_pos = (downhill
                     .map(|downhill_pos| downhill_pos/*.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| e / sz as i32)*/)
-                    .unwrap_or(pos)
+                    .unwrap_or(pos + Vec2::new(1, 0))
                     - pos)/* * scale*/
                     + pos;
                 let downhill_alt = sampler
diff --git a/world/src/block/mod.rs b/world/src/block/mod.rs
index de87f2d28f..7b1cef372d 100644
--- a/world/src/block/mod.rs
+++ b/world/src/block/mod.rs
@@ -61,9 +61,17 @@ impl<'a> BlockGen<'a> {
                 {
                     let cliff_pos3d = Vec3::from(*cliff_pos);
 
+                    // Conservative range of height: [15.70, 49.33]
                     let height = (RandomField::new(seed + 1).get(cliff_pos3d) % 64) as f32
+                        // [0, 63] / (1 + 3 * [0.12, 1.32]) + 3 =
+                        // [0, 63] / (1 + [0.36, 3.96]) + 3 =
+                        // [0, 63] / [1.36, 4.96] + 3 =
+                        // [0, 63] / [1.36, 4.96] + 3 =
+                        // (height min) [0, 0] + 3 = [3, 3]
+                        // (height max) [12.70, 46.33] + 3 = [15.70, 49.33]
                         / (1.0 + 3.0 * cliff_sample.chaos)
                         + 3.0;
+                    // COnservative range of radius: [8, 47]
                     let radius = RandomField::new(seed + 2).get(cliff_pos3d) % 48 + 8;
 
                     max_height.max(
@@ -203,6 +211,7 @@ impl<'a> BlockGen<'a> {
                             world.sim().gen_ctx.fast_turb_x_nz.get(wposf.div(25.0)) as f32,
                             world.sim().gen_ctx.fast_turb_y_nz.get(wposf.div(25.0)) as f32,
                         ) * 8.0;
+                        // let turb = Lerp::lerp(0.0, turb, warp_factor);
 
                         let wpos_turb = Vec2::from(wpos).map(|e: i32| e as f32) + turb;
                         let cliff_height = Self::get_cliff_height(
diff --git a/world/src/column/mod.rs b/world/src/column/mod.rs
index 5d6719f938..a2e061e6db 100644
--- a/world/src/column/mod.rs
+++ b/world/src/column/mod.rs
@@ -583,7 +583,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             .gen_ctx
             .small_nz
             .get((wposf_turb.div(128.0)).into_array()) as f32)
-            .mul(24.0);
+            .mul(4.0);
 
         let alt_for_river = alt
             + if overlap_count == 0.0 {
@@ -806,7 +806,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         } else {
             (false, alt_for_river, downhill_water_alt, 1.0)
         };
-        let warp_factor = 0.0;
+        // let warp_factor = 0.0;
 
         let riverless_alt_delta = (sim
             .gen_ctx
@@ -1117,7 +1117,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 .powf(2.0)
                 .neg()
                 .add(1.0)
-                .mul((1.15 - chaos).min(1.0))
+                .mul((1.32 - chaos).min(1.0))
         };
         let cave_xy = cave_at(wposf);
         let cave_alt = alt - 24.0
diff --git a/world/src/sim/erosion.rs b/world/src/sim/erosion.rs
index 91bddc4906..9ea7e3dc05 100644
--- a/world/src/sim/erosion.rs
+++ b/world/src/sim/erosion.rs
@@ -9,7 +9,9 @@ use rayon::prelude::*;
 use std::{
     cmp::{Ordering, Reverse},
     collections::BinaryHeap,
-    f32, f64, mem, u32,
+    f32, f64, mem,
+    path::PathBuf,
+    u32,
 };
 use vek::*;
 
@@ -495,9 +497,9 @@ pub fn get_rivers(
 ///
 /// TODO: See if allocating in advance is worthwhile.
 fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
-    const MIN_MAX_ANGLE: f64 = 15.0/*6.0*//*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI;
-    const MAX_MAX_ANGLE: f64 = 60.0/*54.0*//*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI;
-    const MAX_ANGLE_RANGE: f64 = MAX_MAX_ANGLE - MIN_MAX_ANGLE;
+    let min_max_angle = (15.0/*6.0*//*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI).tan();
+    let max_max_angle = (45.0/*54.0*//*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI).tan();
+    let max_angle_range = max_max_angle - min_max_angle;
     let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
     h.par_iter()
         .enumerate()
@@ -505,7 +507,7 @@ fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>>
             let wposf = uniform_idx_as_vec2(posi).map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
             let wposz = z as f64 / height_scale;// * CONFIG.mountain_scale as f64;
             // Normalized to be between 6 and and 54 degrees.
-            let div_factor = /*32.0*//*16.0*//*64.0*//*256.0*/8.0/*8.0*//*1.0*//*4.0*/ * /*1.0*/16.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */;
+            let div_factor = /*32.0*//*16.0*//*64.0*//*256.0*/8.0/*8.0*//*1.0*//*4.0*//* * /*1.0*/16.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */*/;
             let rock_strength = rock_strength_nz
                 .get([
                     wposf.x, /* / div_factor*/
@@ -548,7 +550,7 @@ fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>>
                     + 1.0 * log_odds((wposz / CONFIG.mountain_scale as f64).abs().min(dmax).max(dmin)),
             );
             // let rock_strength = 0.5;
-            let max_slope = (rock_strength * MAX_ANGLE_RANGE + MIN_MAX_ANGLE).tan();
+            let max_slope = rock_strength * max_angle_range + min_max_angle;
             // let max_slope = /*30.0.to_radians().tan();*/3.0.sqrt() / 3.0;
             max_slope
         })
@@ -607,12 +609,22 @@ fn get_max_slope(h: &[/*f32*/Alt], rock_strength_nz: &(impl NoiseFn<Point3<f64>>
 ///
 /// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90
 ///
+/// The  approximate equation for soil production function (predictng the rate at which bedrock
+/// turns into soil, increasing the distance between the basement and altitude) is taken from
+/// equation (11) from [2].  This (among numerous other sources) also includes at least one
+/// prediction that hillslope diffusion should be nonlinear, which we sort of attempt to
+/// approximate.
+///
 /// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al..
 ///     Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
 ///     Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175.
 ///     ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩
 ///
-///
+/// [2] William E. Dietrich, Dino G. Bellugi, Leonard S. Sklar, Jonathan D. Stock
+///     Geomorphic Transport Laws for Predicting Landscape Form and Dynamics.
+///     Prediction in Geomorphology, Geophysical Monograph 135.
+///     Copyright 2003 by the American Geophysical Union
+///     10.1029/135GM09
 fn erode(
     h: &mut [Alt],
     b: &mut [Alt],
@@ -626,14 +638,17 @@ fn erode(
     _seed: &RandomField,
     rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
     uplift: impl Fn(usize) -> f32 + Sync,
-    kf: impl Fn(usize) -> f32,
-    kd: impl Fn(usize) -> f32,
+    n_f: impl Fn(usize) -> f32 + Sync,
+    m_f: impl Fn(usize) -> f32 + Sync,
+    kf: impl Fn(usize) -> f64 + Sync,
+    kd: impl Fn(usize) -> f64,
     g: impl Fn(usize) -> f32 + Sync,
     is_ocean: impl Fn(usize) -> bool + Sync,
 ) {
     let compute_stats = true;
     log::debug!("Done draining...");
     let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
+    let min_erosion_height = 0.0;//-f64::INFINITY as Alt;
     let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
     // Minimum sediment thickness before we treat erosion as sediment based.
     let sediment_thickness = 1.0;
@@ -656,11 +671,58 @@ fn erode(
     // 5e-7
     let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / /*5.010e-4*/1e-3;
     println!("dt={:?}", dt);
+    let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+    let chunk_area = neighbor_coef.x * neighbor_coef.y;
+    let min_length = neighbor_coef.reduce_partial_min();
+    let max_stable = /*max_slope * */min_length * min_length / (dt/* / 2.0*/);//1.0/* + /*max_uplift as f64 / dt*/sed / dt*/;
     // Landslide constant: ideally scaled to 10e-2 m / y^-1
-    let l = /*200.0 * max_uplift as f64;*/1.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt;
+    let l = /*200.0 * max_uplift as f64;*/(1.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale);
+    let l_tot = l * dt;
+    // ε₀ = 0.000268 m/y, α = 0.03 (1/m).  This is part of the soil production approximate
+    // equation:
+    //
+    // -∂z_b / ∂t = ε₀ * e^(-αH)
+    //
+    // where
+    //    z_b is the elevation of the soil-bedrock interface (i.e. the basement),
+    //    ε₀ is the production rate of exposed bedrock (H = 0),
+    //    H is the soil thickness normal to the ground surface,
+    //    and α is a parameter (units of 1 / length).
+    //
+    // Note that normal depth at i, for us, will be interpreted as the soil depth vector,
+    //   sed_i = (0, h_i - b_i),
+    // projected onto the bedrock slope vector,
+    //   bedrock_surface_i = (||p_i - p_j||, b_i - b_j),
+    // yielding the soil depth vector
+    //   H_i = sed_i - sed_i ⋅ bedrock_surface_i / (bedrock_surface_i ⋅ bedrock_surface_i) * bedrock_surface_i
+    //
+    //       = (0, h_i - b_i) -
+    //         (0 * ||p_i - p_j|| + (h_i - b_i) * (b_i - b_j)) / (||p_i - p_j||^2 + (b_i - b_j)^2)
+    //         * (||p_i - p_j||, b_i - b_j)
+    //       = (0, h_i - b_i) -
+    //         ((h_i - b_i) * (b_i - b_j)) / (||p_i - p_j||^2 + (b_i - b_j)^2)
+    //         * (||p_i - p_j||, b_i - b_j)
+    //       = (h_i - b_i) *
+    //         ((0, 1) - (b_i - b_j) / (||p_i - p_j||^2 + (b_i - b_j)^2) * (||p_i - p_j||, b_i - b_j))
+    //   H_i_fact = (b_i - b_j) / (||p_i - p_j||^2 + (b_i - b_j)^2)
+    //   H_i = (h_i - b_i) * (((0, 1) - H_i_fact * (||p_i - p_j||, b_i - b_j)))
+    //       = (h_i - b_i) * (-H_i_fact * ||p_i - p_j||, 1 - H_i_fact * (b_i - b_j))
+    //   ||H_i|| = (h_i - b_i) * √((H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact * (b_i - b_j))^2))
+    //
+    // where j is i's receiver and ||p_i - p_j|| is the horizontal displacement between them.  The
+    // idea here is that we first compute the hypotenuse between the horizontal and vertical
+    // displacement of bedrock (getting the horizontal component of the triangle), and then this is
+    // taken as one of the non-hypotenuse sides of the triangle whose other non-hypotenuse side is
+    // the normal height H_i, while their square adds up to the vertical displacement (h_i - b_i).
+    // If h and b have different slopes, this may not work completely correctly, but this is
+    // probably fine as an approximation.
+    let epsilon_0 = 2.68e-4;
+    let alpha = 3e-2;
+    let epsilon_0_tot = epsilon_0 * dt;
     // Net precipitation rate (m / year)
     let p = 1.0 * height_scale;
-    let m = 0.4;
+    /* let n = 2.4;// 1.0;//1.5;//2.4;//1.0;
+    let m = n * 0.5;// n * 0.4;// 0.96;// 0.4;//0.6;//0.96;//0.4; */
     // Stream power erosion constant (bedrock), in m^(1-2m) / year  (times dt).
     let k_fb = // erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
         // 2.244*(5.010e-4)/512*5- (1.097e-5)
@@ -711,9 +773,38 @@ fn erode(
     );
 
     assert!(h.len() == dh.len() && dh.len() == area.len());
+
+    // Precompute factors for Stream Power Law.
+    let k_fact = dh.par_iter().enumerate()
+        .map(|(posi, &posj)| {
+            if posj < 0 {
+                // Egress with no outgoing flows, no stream power.
+                0.0
+            } else {
+                let posj = posj as usize;
+                let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
+                let neighbor_distance = (neighbor_coef * dxy).magnitude();
+                let old_b_i = b[posi];
+                let sed = (ht[posi] - old_b_i) as f64;
+                let k = if sed > sediment_thickness {
+                    // Sediment
+                    // k_fs
+                    k_fs_mult * kf(posi)
+                } else {
+                    // Bedrock
+                    // k_fb
+                    kf(posi)
+                } * dt;
+                let n = n_f(posi) as f64;
+                let m = m_f(posi) as f64;
+
+                k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance.powf(n)
+            }
+        })
+        .collect::<Vec<f64>>();
+    log::info!("Computed stream power factors...");
+
     // max angle of slope depends on rock strength, which is computed from noise function.
-    let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
-    let chunk_area = neighbor_coef.x * neighbor_coef.y;
     // TODO: Make more principled.
     let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan();//1.0;
 
@@ -727,6 +818,7 @@ fn erode(
     let mut err = 2.0 * tol;
 
     // Some variables for tracking statistics, currently only for debugging purposes.
+    let mut minh = f64::INFINITY as Alt;
     let mut maxh = 0.0;
     let mut nland = 0usize;
     let mut sums = 0.0;
@@ -747,6 +839,7 @@ fn erode(
 
         // Reset statistics in each loop.
         maxh = 0.0;
+        minh = f64::INFINITY as Alt;
         nland = 0usize;
         sums = 0.0;
         sumh = 0.0;
@@ -774,6 +867,9 @@ fn erode(
         );
 
         // sum the erosion in stack order
+        //
+        // After:
+        // deltah_i = Σ{j ∈ {i} ∪ upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_j(t + Δt, k))
         for &posi in newh.iter().rev() {
             let posi = posi as usize;
             let posj = dh[posi];
@@ -820,12 +916,19 @@ fn erode(
             } else {
                 let uplift_i = uplift(posi) as Alt;
                 assert!(uplift_i.is_normal() && uplift_i > 0.0 || uplift_i == 0.0);
+                // One side of nonlinear equation (23):
+                //
+                // h_i(t) + U_i * Δt + G / (p̃ * Ã_i) * Σ{j ∈ upstream_i(t)}(h_j(t, FINAL) + U_j * Δt - h_j(t + Δt, k))
+                //
+                // where
+                //
+                // Ã_i = A_i / (∆x∆y) = N_i, number of cells upstream of cell i.
                 *elev = (ht[posi] + uplift_i) as Compute + (deltah[posi] - ((ht[posi] + uplift_i) as Compute - hp[posi])) * g(posi) as Compute / area[posi] as Compute;
             }
         });
 
         // Iterate in ascending height order.
-        let mut sum_err = 0.0f64;
+        let mut sum_err = 0.0 as Compute;
         for &posi in &*newh {
             let posi = posi as usize;
             let old_h_i = /*h*/elev[posi] as f64;
@@ -837,8 +940,11 @@ fn erode(
                 if posj == -1 {
                     panic!("Disconnected lake!");
                 }
-                // wh for oceans is always 0.
-                wh[posi] = 0.0;
+                if ht[posi] > 0.0 {
+                    log::warn!("Ocean above zero?");
+                }
+                // wh for oceans is always at least min_erosion_height.
+                wh[posi] = min_erosion_height.max(ht[posi]);
                 lake_sill[posi] = posi as isize;
                 lake_water_volume[posi] = 0.0;
                 // max_slopes[posi] = kd(posi);
@@ -846,11 +952,11 @@ fn erode(
             } else {
                 // *is_done.at(posi) = done_val;
                 let posj = posj as usize;
-                let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
+                // let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
 
                 // Has an outgoing flow edge (posi, posj).
                 // flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1
-                let neighbor_distance = (neighbor_coef * dxy).magnitude();
+                // let neighbor_distance = (neighbor_coef * dxy).magnitude();
                 // Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5,
                 // we have meters^(1) / meters^(1), so they should cancel out.  Otherwise, we would
                 // want to multiply neighbor_distance by height_scale and area[posi] by
@@ -863,34 +969,82 @@ fn erode(
                 // Therefore, we can rely on wh being set to the water height for that node.
                 let h_j = h[posj] as f64;
                 let wh_j = wh[posj] as f64;
-                let mut new_h_i = old_h_i/* + uplift_i*/;
+                let mut new_h_i = /*old_h_i*/h[posi] as f64/* + uplift_i*/;
                 // Only perform erosion if we are above the water level of the previous node.
-                if new_h_i > wh_j {
+                if old_h_i > wh_j {
                     // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
-                    let k = if sed > sediment_thickness {
+                    /* let k = if sed > sediment_thickness {
                         // Sediment
                         // k_fs
-                        k_fs_mult * kf(posi) as f64
+                        k_fs_mult * kf(posi)
                     } else {
                         // Bedrock
                         // k_fb
-                        kf(posi) as f64
+                        kf(posi)
                     } * dt;
                     // let k = k * uplift_i / max_uplift as f64;
-                    let flux = k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance;
-                    assert!(flux.is_normal() && flux.is_positive());
-                    new_h_i = (new_h_i + flux * h_j) / (1.0 + flux);
+                    let n = n_f(posi) as f64;
+                    let m = m_f(posi) as f64; */
+                    let n = n_f(posi) as f64;
+
+                    if /*n == 1.0*/(n - 1.0).abs() <= 1.0e-3/*f64::EPSILON*/ {
+                        let flux = /*k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance;*/k_fact[posi];
+                        assert!(flux.is_normal() && flux.is_positive() || flux == 0.0);
+                        new_h_i = (/*new_h_i*/old_h_i + flux * h_j) / (1.0 + flux);
+                    } else {
+                        // Local Newton-Raphson
+                        let omega = 0.875f64 / n;
+                        let tolp = 1.0e-3;
+                        let mut errp = 2.0 * tolp;
+                        let h0 = old_h_i;
+                        let fact = k_fact[posi];// k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance.powf(n);
+                        while errp > tolp {
+                            let mut f = new_h_i - h0;
+                            let mut df = 1.0;
+                            // Want: h_i(t+Δt) = h0 - fact * (h_i(t+Δt) - h_j(t+Δt))^n
+                            // Goal: h_i(t+Δt) - h0 + fact * (h_i(t+Δt) - h_j(t+Δt))^n = 0
+                            f += fact * 0.0.max(new_h_i - h_j).powf(n);
+                            // ∂h_i(t+Δt)/∂n = 1 + fact * n * (h_i(t+Δt) - h_j(t+Δt))^(n - 1)
+                            df += fact * n * 0.0.max(new_h_i - h_j).powf(n - 1.0);
+                            // hn = h_i(t+Δt, k) - (h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) - h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)
+                            let hn = new_h_i - f / df;
+                            // errp = |(h_i(t+Δt, k) - (h0 - fact * (h_i(t+Δt, k) - h_j(t+Δt))^n)) / ∂h_i/∂n(t+Δt, k)|
+                            errp = (hn - new_h_i).abs();
+                            // h_i(t+∆t, k+1) = ...
+                            new_h_i = new_h_i * (1.0 - omega) + hn * omega;
+                        }
+                        /* omega=0.875d0/n
+                        tolp=1.d-3
+                        errp=2.d0*tolp
+                        h0=elev(ijk)
+                        do while (errp.gt.tolp)
+                          f=h(ijk)-h0
+                          df=1.d0
+                          if (ht(ijk).gt.ht(ijr)) then
+                            fact = kfint(ijk)*dt*a(ijk)**m/length(ijk)**n
+                            f=f+fact*max(0.d0,h(ijk)-h(ijr))**n
+                            df=df+fact*n*max(0.d0,h(ijk)-h(ijr))**(n-1.d0)
+                          endif
+                          hn=h(ijk)-f/df
+                          errp=abs(hn-h(ijk))
+                          h(ijk)=h(ijk)*(1.d0-omega)+hn*omega
+                        enddo */
+                    }
                     lake_sill[posi] = posi as isize;
                     lake_water_volume[posi] = 0.0;
 
                     /* // Thermal erosion (landslide)
-                    let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+                    let dz = (new_h_i - /*h_j*//*h_k*//*wh_j*/h_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                     let mag_slope = dz/*.abs()*/ / neighbor_distance;
                     let max_slope = max_slopes[posi] as f64;
                     if mag_slope > max_slope {
                         let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
-                        new_h_i = (ht[posi] as f64 + l * (mag_slope - max_slope));
-                        if compute_stats && new_h_i > wh_j {
+                        // new_h_i = (ht[posi] as f64 + l_tot * (mag_slope - max_slope));
+                        // new_h_i = new_h_i - l_tot * (mag_slope - max_slope);
+                        // new_h_i = new_h_i - l_tot * (mag_slope - max_slope);
+                        // new_h_i = hp[posj] + dh;
+                        new_h_i = /*old_h_i.max*/(/*wh_j*//*ht[posi] as Compute*//*h_j*/hp[posj]/*ht[posj] as Compute*/ + dh).max(new_h_i - l_tot * (mag_slope - max_slope));
+                        if compute_stats/* && new_h_i > wh_j*/ {
                             ntherm += 1;
                         }
                     } */
@@ -900,7 +1054,9 @@ fn erode(
                     if new_h_i <= wh_j {
                         new_h_i = wh_j;
                     } else {
-                        if compute_stats {
+                        if compute_stats && new_h_i > 0.0 {
+                            let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
+                            let neighbor_distance = (neighbor_coef * dxy).magnitude();
                             let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                             let mag_slope = dz/*.abs()*/ / neighbor_distance;
 
@@ -911,6 +1067,7 @@ fn erode(
                         }
                     }
                 } else {
+                    new_h_i = old_h_i;
                     let lposj = lake_sill[posj];
                     lake_sill[posi] = lposj;
                     if lposj >= 0 {
@@ -963,7 +1120,7 @@ fn erode(
                         // max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
                         // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
                         let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
-                        new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
+                        new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l_tot * (mag_slope - max_slope));
                         let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                         let slope = dz/*.abs()*/ / neighbor_distance;
                         sums += slope;
@@ -984,8 +1141,12 @@ fn erode(
             }
             // *is_done.at(posi) = done_val;
             if compute_stats {
-                maxh = h[posi].max(maxh);
                 sumsed += sed;
+                let h_i = h[posi];
+                if h_i > 0.0 {
+                    minh = h_i.min(minh);
+                }
+                maxh = h_i.max(maxh);
             }
 
             // Add sum of squares of errors.
@@ -1004,12 +1165,58 @@ fn erode(
     //b=min(h,b)
 
     // update the basement
-    b.par_iter_mut().zip(h.par_iter()).enumerate().for_each(|(posi, (mut b, h))| {
+    //
+    // NOTE: Despite this not quite applying since basement order and height order differ, we once
+    // again borrow the implicit FastScape stack order.  If this becomes a problem we can easily
+    // compute a separate stack order just for basement.
+    for &posi in &*newh {
+        let posi = posi as usize;
+        let old_b_i = b[posi];
+        let h_i = h[posi];
+        let uplift_i = uplift(posi) as Alt;
+
+        // First, add uplift...
+        let mut new_b_i = (old_b_i + uplift_i).min(h_i);
+
+        let posj = dh[posi];
+        // Sediment height normal to bedrock.  NOTE: Currently we can actually have sedment and
+        // bedrock slope at different heights, meaning there's no uniform slope.  There are
+        // probably more correct ways to account for this, such as averaging, integrating, or doing
+        // things by mass / volume instead of height, but for now we use the time-honored
+        // technique of ignoring the problem.
+        let h_normal = if posj < 0 {
+            // Egress with no outgoing flows; for now, we assume this means normal and vertical
+            // coincide.
+            (h_i - new_b_i) as f64
+        } else {
+            let posj = posj as usize;
+            let b_j = b[posj];
+            let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
+            let neighbor_distance_squared = (neighbor_coef * dxy).magnitude_squared();
+            let vertical_sed = (h_i - new_b_i) as f64;
+            let db = (new_b_i - b_j) as f64;
+            // H_i_fact = (b_i - b_j) / (||p_i - p_j||^2 + (b_i - b_j)^2)
+            let h_i_fact = db / (neighbor_distance_squared + db * db);
+            let h_i_vertical = 1.0 - h_i_fact * db;
+            // ||H_i|| = (h_i - b_i) * √((H_i_fact^2 * ||p_i - p_j||^2 + (1 - H_i_fact * (b_i - b_j))^2))
+            vertical_sed * (h_i_fact * h_i_fact * neighbor_distance_squared + h_i_vertical * h_i_vertical).sqrt()
+        };
+        // Rate of sediment production: -∂z_b / ∂t = ε₀ * e^(-αH)
+        let p_i = epsilon_0_tot * f64::exp(-alpha * h_normal);
+        // println!("h_normal = {:?}, p_i = {:?}", h_normal, p_i);
+
+        new_b_i -= p_i as Alt;
+
+        b[posi] = new_b_i;
+    }
+    log::info!("Done updating basement and applying soil production...");
+
+    /* b.par_iter_mut().zip(h.par_iter()).enumerate().for_each(|(posi, (mut b, h))| {
         let old_b_i = *b;
         let uplift_i = uplift(posi) as Alt;
 
         *b = (old_b_i + uplift_i).min(*h);
-    });
+    }); */
 
     // update the height to reflect sediment flux.
     h.par_iter_mut().enumerate().for_each(|(posi, mut h)| {
@@ -1035,11 +1242,12 @@ fn erode(
     // enddo
 
     log::info!(
-        "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})\n        \
+        "Done applying stream power (max height: {:?}) (avg height: {:?}) ((min height: {:?}) avg slope: {:?})\n        \
         (old avg sediment thickness [all/land]: {:?} / {:?})\n        \
         (num land: {:?}) (num thermal: {:?})",
         maxh,
         avgz(sumh, nland),
+        minh,
         avgz(sums, nland),
         avgz(sumsed, newh.len()),
         avgz(sumsed_land, nland),
@@ -1049,6 +1257,7 @@ fn erode(
 
     // Apply thermal erosion.
     maxh = 0.0;
+    minh = f64::INFINITY as Alt;
     sumh = 0.0;
     sums = 0.0;
     sumsed = 0.0;
@@ -1063,6 +1272,7 @@ fn erode(
 
         let max_slope = max_slopes[posi];
         // Remember k_d for this chunk in max_slopes.
+        // higher max_slope => much lower kd_factor.
         let kd_factor =
             // 1.0;
             (1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/2.0))/*.min(kdsed)*/;
@@ -1071,7 +1281,7 @@ fn erode(
             kdsed/* * kd_factor*/
         } else {
             // Bedrock
-            kd(posi) as f64 / kd_factor
+            kd(posi) / kd_factor
         };
 
         let posj = dh[posi];
@@ -1079,13 +1289,14 @@ fn erode(
             if posj == -1 {
                 panic!("Disconnected lake!");
             }
-            wh[posi] = 0.0;
+            // wh for oceans is always at least min_erosion_height.
+            wh[posi] = min_erosion_height.max(ht[posi]);
             // Egress with no outgoing flows.
         } else {
             let posj = posj as usize;
             // Find the water height for this chunk's receiver; we only apply thermal erosion
             // for chunks above water.
-            let wh_j = wh[posj] as f64;
+            let mut wh_j = wh[posj] as f64;
             // If you're on the lake bottom and not right next to your neighbor, don't compute a
             // slope.
             let mut new_h_i = old_h_i;/*old_b_i;*/
@@ -1096,7 +1307,8 @@ fn erode(
                 // totally submerged in water, and that talus angle if part of the substance is
                 // in water is 0 (or the same as the dry part, if this is set to wh_j), but
                 // actually that's probably not true.
-                let h_j = /*h[posj] as f64*/wh_j;
+                let old_h_j = h[posj] as f64;
+                let h_j = /*h[posj] as f64*//*wh_j*/old_h_j;
                 // let h_j = b[posj] as f64;
                 /* let indirection_idx = indirection[posi];
                 let is_lake_bottom = indirection_idx < 0;
@@ -1131,12 +1343,22 @@ fn erode(
                     // max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
                     // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
                     let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
-                    new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*/h_k + dh + l * (mag_slope - max_slope));
-                    // new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
+                    // new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*/ht[posi] as f64 + l_tot * (mag_slope - max_slope));
+                    // new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*/h_k + dh + l_tot * (mag_slope - max_slope));
+                    // new_h_i = /*h_j.max*//*(h_k + dh).max*/(new_h_i - l_tot * (mag_slope - max_slope));
+                    let dtherm = (l_tot * (mag_slope - max_slope)).min((dz - dh)/* / 2.0*/);
+                    new_h_i = /*h_j.max*//*(h_k + dh).max*/(/*new_h_i*//*h_k + dh*/new_h_i - dtherm);
+                    /* let new_h_j = (old_h_j + dtherm).min(old_h_j.max(new_h_i));
+                    h[posj] = new_h_j as Alt;
+                    wh_j = wh_j.max(new_h_j);
+                    wh[posj] = wh_j as Alt; */
+                    // No more hillslope processes on newly exposed bedrock.
+                    max_slopes[posi] = 0.0;
+                    // max_slopes[posi] = l;
                     if new_h_i <= wh_j {
                         new_h_i = wh_j;
                     } else {
-                        if compute_stats {
+                        if compute_stats && new_h_i > 0.0 {
                             let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                             let slope = dz/*.abs()*/ / neighbor_distance;
                             sums += slope;
@@ -1161,6 +1383,10 @@ fn erode(
                         ntherm += 1;
                     }
                 } else {
+                    // Poorly emulating nonlinear hillslope transport as described by
+                    // http://eps.berkeley.edu/~bill/papers/112.pdf.
+                    // sqrt(3)/3*32*32/(128000/2)
+                    max_slopes[posi] = (max_slopes[posi] * 1.0 / (1.0 - (mag_slope / max_slope).powi(2)));
                     /*if kd_factor < 1.0 {
                         max_slopes[posi] *= kd_factor;
                     }*/
@@ -1168,7 +1394,7 @@ fn erode(
                         max_slopes[posi] *= kd_factor;
                     } */
                     // max_slopes[posi] *= kd_factor;
-                    if compute_stats {
+                    if compute_stats && new_h_i > 0.0 {
                         sums += mag_slope;
                         // Just use the computed rate.
                         nland += 1;
@@ -1176,28 +1402,35 @@ fn erode(
                         sumsed_land += sed;
                     }
                 }
-                h/*b*/[posi] = old_h_i.min(new_h_i) as Alt;
+                h/*b*/[posi] = /*old_h_i.min(new_h_i)*/new_h_i as Alt;
                 // Make sure to update the basement as well!
                 // b[posi] = old_b_i.min(new_h_i) as f32;
-                b[posi] = old_b_i.min(old_b_i + (new_h_i - old_h_i)) as Alt;
+                b[posi] = old_b_i.min(old_b_i + (/*old_h_i.min(*/new_h_i/*)*/ - old_h_i)) as Alt;
                 // sumh += new_h_i;
             }
             // Set wh to this node's water height (max of receiver's water height and
             // this node's height).
             wh[posi] = wh_j.max(new_h_i) as Alt;
         }
+        max_slopes[posi] = max_slopes[posi].min(max_stable);
+
         // *is_done.at(posi) = done_val;
         if compute_stats {
             sumsed += sed;
-            maxh = h[posi].max(maxh);
+            let h_i = h[posi];
+            if h_i > 0.0 {
+                minh = h_i.min(minh);
+            }
+            maxh = h_i.max(maxh);
         }
     }
     log::debug!(
-        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?})\n        \
+        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (min height: {:?}) (avg slope: {:?})\n        \
         (avg sediment thickness [all/land]: {:?} / {:?})\n        \
         (num land: {:?}) (num thermal: {:?})",
         maxh,
         avgz(sumh, nland),
+        minh,
         avgz(sums, nland),
         avgz(sumsed, newh.len()),
         avgz(sumsed_land, nland),
@@ -1694,17 +1927,20 @@ pub fn get_lakes<F: Float>(h: impl Fn(usize) -> F, downhill: &mut [isize]) -> (u
 pub fn do_erosion(
     erosion_base: f32,
     _max_uplift: f32,
-    n: usize,
+    n_steps: usize,
     seed: &RandomField,
     rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
     oldh: impl Fn(usize) -> f32 + Sync,
     oldb: impl Fn(usize) -> f32 + Sync,
     is_ocean: impl Fn(usize) -> bool + Sync,
     uplift: impl Fn(usize) -> f32 + Sync,
-    kf: impl Fn(usize) -> f32 + Sync,
-    kd: impl Fn(usize) -> f32 + Sync,
+    n: impl Fn(usize) -> f32 + Sync,
+    theta: impl Fn(usize) -> f32 + Sync,
+    kf: impl Fn(usize) -> f64 + Sync,
+    kd: impl Fn(usize) -> f64 + Sync,
     g: impl Fn(usize) -> f32 + Sync,
 ) -> (Box<[Alt]>, Box<[Alt]>) {
+    log::info!("Initializing erosion arrays...");
     let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y)
         .into_par_iter()
         .map(|posi| oldh(posi) as Alt)
@@ -1716,6 +1952,19 @@ pub fn do_erosion(
         .map(|posi| oldb(posi) as Alt)
         .collect::<Vec<_>>()
         .into_boxed_slice();
+    // Stream power law slope exponent--link between channel slope and erosion rate.
+    let n = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
+        .map(|posi| n(posi))
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
+    // Stream power law concavity index (θ = m/n), turned into an exponent on drainage
+    // (which is a proxy for discharge according to Hack's Law).
+    let m = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
+        .map(|posi| theta(posi) * n[posi])
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
     // Stream power law erodability constant for fluvial erosion (bedrock)
     let kf = (0..WORLD_SIZE.x * WORLD_SIZE.y)
         .into_par_iter()
@@ -1760,6 +2009,7 @@ pub fn do_erosion(
         .max_by(|a, b| a.partial_cmp(&b).unwrap())
         .unwrap();
     log::debug!("Max uplift: {:?}", max_uplift);
+    log::debug!("Max g: {:?}", max_g);
     // Height of terrain, including sediment.
     let mut h = oldh_;
     // 0.01 / 2e-5 = 500
@@ -1778,12 +2028,14 @@ pub fn do_erosion(
         /*1.5e-2*//*1e-4*//*1.25e-2*/1.5e-2 / 1.0 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
         /* * k_fb / 2e-5 */;
     // let kd = |posi: usize| kd_bedrock; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
+    let n = |posi: usize| n[posi];
+    let m = |posi: usize| m[posi];
     let kd = |posi: usize| kd[posi]; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
     let kf = |posi: usize| kf[posi];
     let g = |posi: usize| g[posi];
     // Hillslope diffusion coefficient for sediment.
     let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
-    for i in 0..n {
+    for i in 0..n_steps {
         log::debug!("Erosion iteration #{:?}", i);
         erode(
             &mut h,
@@ -1802,6 +2054,8 @@ pub fn do_erosion(
             seed,
             rock_strength_nz,
             |posi| uplift[posi],
+            n,
+            m,
             kf,
             kd,
             g,
diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs
index 5c2a670315..57a55ea995 100644
--- a/world/src/sim/mod.rs
+++ b/world/src/sim/mod.rs
@@ -36,10 +36,14 @@ use num::{Float, Signed};
 use rand::{Rng, SeedableRng};
 use rand_chacha::ChaChaRng;
 use rayon::prelude::*;
+use serde_derive::{Deserialize, Serialize};
 use std::{
     collections::HashMap,
+    io::{BufReader, BufWriter},
     f32, f64,
+    fs::File,
     ops::{Add, Div, Mul, Neg, Sub},
+    path::PathBuf,
     sync::Arc,
 };
 use vek::*;
@@ -107,19 +111,49 @@ pub(crate) struct GenCtx {
     pub uplift_nz: Worley,
 }
 
+pub enum FileOpts {
+    /// If set, generate the world map and do not try to save to or load from file
+    /// (default).
+    Generate,
+    /// If set, generate the world map and save the world file (path is created
+    /// the same way screenshot paths are).
+    Save,
+    /// If set, load the world file from this path (errors if path not found).
+    Load(PathBuf),
+}
+
+impl Default for FileOpts {
+    fn default() -> Self {
+        Self::Generate
+    }
+}
+
 pub struct WorldOpts {
     /// Set to false to disable seeding elements during worldgen.
     pub seed_elements: bool,
+    pub world_file: FileOpts,
 }
 
 impl Default for WorldOpts {
     fn default() -> Self {
         Self {
             seed_elements: true,
+            world_file: Default::default(),
         }
     }
 }
 
+/// A way to store certain components between runs of map generation.  Only intended for
+/// development purposes--no attempt is made to detect map invalidation or make sure that the map
+/// is synchronized with updates to noise-rs, changes to other parameters, etc.
+#[derive(Serialize,Deserialize)]
+pub struct WorldFile {
+    /// Saved altitude height map.
+    pub alt: Box<[Alt]>,
+    /// Saved basement height map.
+    pub basement: Box<[Alt]>,
+}
+
 pub struct WorldSim {
     pub seed: u32,
     pub(crate) chunks: Vec<SimChunk>,
@@ -226,7 +260,7 @@ impl WorldSim {
             .set_seed(rng.gen()),
             uplift_nz: Worley::new()
                 .set_seed(rng.gen())
-                .set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * 256.0))
+                .set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * /*256.0*/128.0))
                 // .set_displacement(/*0.5*/0.0)
                 .set_displacement(/*0.5*/1.0)
                 .set_range_function(RangeFunction::Euclidean)
@@ -248,13 +282,13 @@ impl WorldSim {
             .set_scale(1.0 / rock_strength_scale); */
 
         let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
-        let max_erosion_per_delta_t = /*8.0*/64.0/*128.0*//*32.0*/ * height_scale;
+        let max_erosion_per_delta_t = /*8.0*/32.0/*32.0*//*128.0*//*16.0*//*128.0*//*32.0*/ * height_scale;
         let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002;
         let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002;
         let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
-        let n_steps = 50; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
-        let n_small_steps = 25;//50;//8;//8;//8;//8;//8; // 8
-
+        let n_steps = 100;//100; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
+        let n_small_steps = 25;//25;//8;//50;//50;//8;//8;//8;//8;//8; // 8
+        let n_post_load_steps = 0;//8
 
         // Logistic regression.  Make sure x ∈ (0, 1).
         let logit = |x: f64| x.ln() - (-x).ln_1p();
@@ -359,7 +393,7 @@ impl WorldSim {
                     .add(0.3)
                     .max(0.0);
 
-                    // chaos produces a value in [0.12, 1.24].  It is a meta-level factor intended to
+                    // chaos produces a value in [0.12, 1.32].  It is a meta-level factor intended to
                     // reflect how "chaotic" the region is--how much weird stuff is going on on this
                     // terrain.
                     Some(
@@ -384,7 +418,7 @@ impl WorldSim {
                         )
                         // Chaos is always increased by a little when we're on a hill (but remember
                         // that hill is 0.3 or less about 50% of the time).
-                        // [0, 1] + 0.15 * [0, 1.6] = [0, 1.24]
+                        // [0, 1] + 0.2 * [0, 1.6] = [0, 1.32]
                         .add(0.2 * hill)
                         // We can't have *no* chaos!
                         .max(0.12)) as f32,
@@ -447,13 +481,13 @@ impl WorldSim {
             };
 
             // Now we can compute the final altitude using chaos.
-            // We multiply by chaos clamped to [0.1, 1.24] to get a value between [0.03, 2.232]
+            // We multiply by chaos clamped to [0.1, 1.32] to get a value between [0.03, 2.232]
             // for alt_pre, then multiply by CONFIG.mountain_scale and add to the base and sea
             // level to get an adjusted value, then multiply the whole thing by map_edge_factor
             // (TODO: compute final bounds).
             //
-            // [-.3675, .3325] + [-0.445, 0.565] * [0.12, 1.24]^1.2
-            // ~ [-.3675, .3325] + [-0.445, 0.565] * [_, 1.30]
+            // [-.3675, .3325] + [-0.445, 0.565] * [0.12, 1.32]^1.2
+            // ~ [-.3675, .3325] + [-0.445, 0.565] * [0.07, 1.40]
             // = [-.3675, .3325] + ([-0.5785, 0.7345])
             // = [-0.946, 1.067]
             Some(
@@ -499,67 +533,70 @@ impl WorldSim {
             .clone()
             .enable_range(true);
         // Recalculate altitudes without oceans.
-        // NaNs in these uniform vectors wherever pure_water() returns true.
-        let ((alt_old_no_ocean, alt_old_inverse), (uplift_uniform, _)) = rayon::join(
-            || {
-                uniform_noise(|posi, _| {
-                    if is_ocean_fn(posi) {
-                        None
-                    } else {
-                        Some(old_height(posi) /*.abs()*/)
-                    }
-                })
-            },
-            || {
-                uniform_noise(|posi, wposf| {
-                    if is_ocean_fn(posi) {
-                        None
-                    } else {
-                        let turb_wposf =
-                            wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
-                        let turb = Vec2::new(
-                            gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
-                            gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                        ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
-                        // let turb = Vec2::zero();
-                        let turb_wposf = wposf + turb;
-                        let turb_wposi = turb_wposf
-                            .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
-                            .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
-                        let turb_posi = vec2_as_uniform_idx(turb_wposi);
-                        let udist = uplift_nz_dist.get(turb_wposf.into_array())
-                            .min(1.0)
-                            .max(-1.0)
-                            .mul(0.5)
-                            .add(0.5);
-                        let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
-                            /* .min(0.5)
-                            .max(-0.5)*/
-                            .min(1.0)
-                            .max(-1.0)
-                            .mul(0.5)
-                            .add(0.5);
-                        let oheight = alt_old[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5;
-                        assert!(udist >= 0.0);
-                        assert!(udist <= 1.0);
-                        let uheight_1 = uheight;//.powf(2.0);
-                        let udist_1 = (0.5 - udist).mul(2.0).max(0.0);
-                        let udist_2 = udist.mul(2.0).min(1.0);
-                        let height =
-                            // uheight_1;
-                            // uheight_1 * (/*udist_2*/udist.powf(2.0) * (f64::consts::PI * 2.0 * (1.0 / (1.0 - udist).max(f64::EPSILON)).min(2.5)/*udist * 5.0*/ * 2.0).cos().mul(0.5).add(0.5));
-                            // uheight * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin()
-                            /* (uheight /* * 0.8*//* * udist_1.powf(2.0)*/ +
-                            /*exp_inverse_cdf*/(oheight./*max(0.0).min(max_epsilon).abs()*/)/*.powf(2.0)*/ * 0.2/* * udist_2.powf(2.0)*/; */
-                            (uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
-                            // * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
-                        // let height = uheight * (0.5 - udist) * 0.8 + (oheight.signum() * oheight.max(0.0).abs().powf(2.0)) * 0.2;// * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
+        // NaNs in these uniform vectors wherever is_ocean_fn returns true.
+        let (alt_old_no_ocean, alt_old_inverse) =
+            uniform_noise(|posi, _| {
+                if is_ocean_fn(posi) {
+                    None
+                } else {
+                    Some(old_height(posi) /*.abs()*/)
+                }
+            });
+        let (uplift_uniform, _) =
+            uniform_noise(|posi, wposf| {
+                if is_ocean_fn(posi) {
+                    None
+                } else {
+                    let turb_wposf =
+                        wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+                    let turb = Vec2::new(
+                        gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                        gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                    ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                    // let turb = Vec2::zero();
+                    let turb_wposf = wposf + turb;
+                    let turb_wposi = turb_wposf
+                        .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
+                        .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
+                    let turb_posi = vec2_as_uniform_idx(turb_wposi);
+                    let udist = uplift_nz_dist.get(turb_wposf.into_array())
+                        .min(1.0)
+                        .max(-1.0)
+                        .mul(0.5)
+                        .add(0.5);
+                    let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                        /* .min(0.5)
+                        .max(-0.5)*/
+                        .min(1.0)
+                        .max(-1.0)
+                        .mul(0.5)
+                        .add(0.5);
+                    let oheight = /*alt_old*//*alt_base*/alt_old_no_ocean[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5;
+                    assert!(udist >= 0.0);
+                    assert!(udist <= 1.0);
+                    let uheight_1 = uheight;//.powf(2.0);
+                    let udist_1 = (0.5 - udist).mul(2.0).max(0.0);
+                    let udist_2 = udist.mul(2.0).min(1.0);
+                    let height =
+                        // 1.0;
+                        // uheight_1;
+                        // uheight_1 * (/*udist_2*/udist.powf(2.0) * (f64::consts::PI * 2.0 * (1.0 / (1.0 - udist).max(f64::EPSILON)).min(2.5)/*udist * 5.0*/ * 2.0).cos().mul(0.5).add(0.5));
+                        // uheight * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin()
+                        // uheight;
+                        // (0.8 * uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
+                        // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * udist_2.powf(2.0)).max(0.0).min(1.0);
+                        // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * 0.2).max(0.0).min(1.0);
+                        // (0.8 * uheight * udist_1 + 0.8 * udist_2 + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
+                        /* uheight * 0.8 * udist_1.powf(2.0) +
+                        /*exp_inverse_cdf*/(oheight/*.max(0.0).min(max_epsilon).abs()*/).powf(2.0) * 0.2 * udist_2.powf(2.0); */
+                        // (uheight + oheight.powf(2.0) * 0.05).max(0.0).min(1.0);
+                        (uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0);
+                        // * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
+                    // let height = uheight * (0.5 - udist) * 0.8 + (oheight.signum() * oheight.max(0.0).abs().powf(2.0)) * 0.2;// * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
 
-                        Some(height)
-                    }
-                })
-            },
-        );
+                    Some(height)
+                }
+            });
 
         let old_height_uniform = |posi: usize| alt_old_no_ocean[posi].0;
         let alt_old_min_uniform = 0.0;
@@ -618,10 +655,53 @@ impl WorldSim {
         ).powf(erosion_pow).mul(0.5))*/; */
         let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(/*exp_inverse_cdf*//*logit*/inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
 ).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*//*.powf(2.0)*/;
+        let n_func = |posi| {
+            if is_ocean_fn(posi) {
+                return 1.0;
+            }
+            let wposf = (uniform_idx_as_vec2(posi)
+                * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+            .map(|e| e as f64);
+            let turb_wposf =
+                wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+            let turb = Vec2::new(
+                gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+            ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+            // let turb = Vec2::zero();
+            let turb_wposf = wposf + turb;
+            let turb_wposi = turb_wposf
+                .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
+                .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
+            let turb_posi = vec2_as_uniform_idx(turb_wposi);
+            let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                /* .min(0.5)
+                .max(-0.5)*/
+                .min(1.0)
+                .max(-1.0)
+                .mul(0.5)
+                .add(0.5);
+            /* if uheight > 0.8 {
+                1.5
+            } else {
+                1.0
+            } */
+            // ((1.5 - 0.6) * uheight + 0.6) as f32
+            // ((1.5 - 1.0) * uheight + 1.0) as f32
+            // ((3.5 - 1.5) * (1.0 - uheight) + 1.5) as f32
+            1.0
+        };
+        let theta_func = |posi| {
+            0.4
+        };
         let kf_func = {
             |posi| {
                 if is_ocean_fn(posi) {
                     return 1.0e-4;
+                    // return 2.0e-5;
+                    // return 2.0e-6;
+                    // return 2.0e-10;
+                    // return 0.0;
                 }
                 let wposf = (uniform_idx_as_vec2(posi)
                     * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
@@ -650,7 +730,23 @@ impl WorldSim {
                 // kf = 2e-5: normal (dike [unexposed])
                 // kf = 1e-6: normal-low (dike [exposed])
                 // kf = 2e-6: low (mountain)
-                ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32
+                // --
+                // kf = 2.5e-7 to 8e-7: very low (Cordonnier papers on plate tectonics)
+                // ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32
+                //
+                // ACTUAL recorded values worldwide: much lower...
+                //
+                // Or maybe not?  Getting something like 2e-3...
+                //
+                // ...or 8.345e5.
+                // ((1.0 - uheight) * (5e-5 - 9.88e-15) + 9.88e-15)
+                // ((1.0 - uheight) * (1.5e-4 - 9.88e-15) + 9.88e-15)
+                // ((1.0 - uheight) * (8.345e5 - 2.0e-6) + 2.0e-6) as f32
+                ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6)
+                // 2e-5
+                // 2.9e-10
+                // ((1.0 - uheight) * (5e-5 - 2.9e-10) + 2.9e-10)
+                // ((1.0 - uheight) * (5e-5 - 2.9e-14) + 2.9e-14)
             }
         };
         let kd_func = {
@@ -684,26 +780,38 @@ impl WorldSim {
                 // kd = 1.5e-2: normal-high (plateau [fan sediment])
                 // kd = 1e-2: normal (plateau)
                 1.0e-2
-                // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2) as f32
+                // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2)
             }
         };
         let g_func =
             |posi| {
                 if /*is_ocean_fn(posi)*/map_edge_factor(posi) == 0.0 {
                     return 0.0;
+                    // return 5.0;
                 }
                 let wposf = (uniform_idx_as_vec2(posi)
                     * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
                 .map(|e| e as f64);
-                /* let turb_wposf =
+                let turb_wposf =
                     wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
                 let turb = Vec2::new(
                     gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
                     gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
                 ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
                 // let turb = Vec2::zero();
-                let turb_wposf = wposf + turb; */
-                let turb_wposf = wposf;
+                let turb_wposf = wposf + turb;
+                let turb_wposi = turb_wposf
+                    .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
+                    .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
+                let turb_posi = vec2_as_uniform_idx(turb_wposi);
+                let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                    /* .min(0.5)
+                    .max(-0.5)*/
+                    .min(1.0)
+                    .max(-1.0)
+                    .mul(0.5)
+                    .add(0.5);
+
                 let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
                     .min(1.0)
                     .max(-1.0)
@@ -711,7 +819,7 @@ impl WorldSim {
                     .add(0.5); */
                     chaos[posi].1;
 
-                assert!(uchaos <= 1.24);
+                assert!(uchaos <= 1.32);
 
                 // G = d* v_s / p_0, where
                 //  v_s is the settling velocity of sediment grains
@@ -722,15 +830,21 @@ impl WorldSim {
                 //  and washed loads.
                 //
                 // G is typically on the order of 1 or greater.  However, we are only guaranteed to
-                // converge for G ≤ 1, so we keep it in the chaos range of [0.12, 1.24].
-                ((1.24 - uchaos) / 1.24).powf(0.75) * 1.24
+                // converge for G ≤ 1, so we keep it in the chaos range of [0.12, 1.32].
+                ((1.32 - uchaos) / 1.32).powf(0.75) * 1.32
+                // ((1.32 - 0.12) * (1.0 - uheight) + 0.12) as f32
+                // 1.1 * (1.0 - uheight) as f32
+                // 1.0 * (1.0 - uheight) as f32
                 // 1.0
+                // 0.0
                 // 1.0
                 // 1.5
             };
         let uplift_fn =
             |posi| {
                 if is_ocean_fn(posi) {
+                    /* return 1e-2
+                        .mul(max_erosion_per_delta_t) as f32; */
                     return 0.0;
                 }
                 let wposf = (uniform_idx_as_vec2(posi)
@@ -838,8 +952,8 @@ impl WorldSim {
                     .mul(0.5)
                     .add(0.5);
                 // u = 1e-3: normal-high (dike, mountain)
-                // u = 0.5: normal (mid example in Yuan, average mountain uplift)
-                // u = 0.2: low (low example in Yuan; known that lagoons etc. may have u ~ 0.05).
+                // u = 5e-4: normal (mid example in Yuan, average mountain uplift)
+                // u = 2e-4: low (low example in Yuan; known that lagoons etc. may have u ~ 0.05).
                 // u = 0: low (plateau [fan, altitude = 0.0])
                 // let height = uheight;
                 // let height = 1.0f64;
@@ -855,6 +969,8 @@ impl WorldSim {
                     .add(1.0 / 16.0) */
                     /* .mul(5.0 / 8.0)
                     .add(3.0 / 8.0) */
+                    /* .mul(6.0 / 8.0)
+                    .add(2.0 / 8.0) */
                     /* .mul(7.0 / 8.0)
                     .add(1.0 / 8.0) */
                     .mul(max_erosion_per_delta_t)
@@ -907,8 +1023,10 @@ impl WorldSim {
                             .mul(0.045)*/)
                 };
 
-                /* 0.0 */
-                (old_height_uniform(posi)/*.powf(2.0)*/ - 0.5)/* * CONFIG.mountain_scale as f32*/
+                (old_height_uniform(posi) - 0.5)/* * max_erosion_per_delta_t as f32*/
+                // uplift_fn(posi) * (CONFIG.mountain_scale / max_erosion_per_delta_t as f32)
+                // 0.0
+                // /*-CONFIG.mountain_scale * 0.5 + *//*-CONFIG.mountain_scale/* * 0.75*/ + */(old_height_uniform(posi)/*.powf(2.0)*/ - 0.5)/* * CONFIG.mountain_scale as f32*/
                 // uplift_fn(posi) * (CONFIG.mountain_scale / max_erosion_per_delta_t as f32)
                 // 0.0
                 /* // 0.0
@@ -923,35 +1041,145 @@ impl WorldSim {
                 // 5.0 / CONFIG.mountain_scale */
             }
         };
-        let (alt, basement) = do_erosion(
-            0.0,
-            max_erosion_per_delta_t as f32,
-            n_steps,
-            &river_seed,
-            &rock_strength_nz,
-            |posi| alt_func(posi),// + if is_ocean_fn(posi) { 0.0 } else { 128.0 },
-            |posi| alt_func(posi) - if is_ocean_fn(posi) { 0.0 } else { 1400.0 },// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
-            is_ocean_fn,
-            uplift_fn,
-            |posi| kf_func(posi),
-            |posi| kd_func(posi),
-            |posi| g_func(posi),
-        );
-        // Quick "small scale" erosion cycle in order to lower extreme angles.
-        let (alt, basement) = do_erosion(
-            0.0,
-            (1.0 * height_scale) as f32,
-            n_small_steps,
-            &river_seed,
-            &rock_strength_nz,
-            |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi] as f32,
-            |posi| basement[posi] as f32,
-            is_ocean_fn,
-            |posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
-            |posi| kf_func(posi),
-            |posi| kd_func(posi),
-            |posi| g_func(posi),
-        );
+
+        /* // FIXME: Remove.
+        let is_ocean = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+            .into_par_iter()
+            .map(|i| map_edge_factor(i) == 0.0)
+            .collect::<Vec<_>>();
+        let is_ocean_fn = |posi: usize| is_ocean[posi]; */
+
+        // Load map, if necessary.
+        let parsed_world_file = (|| {
+            if let FileOpts::Load(ref path) = opts.world_file {
+                let file = match File::open(path) {
+                    Ok(file) => file,
+                    Err(err) => {
+                        log::warn!("Couldn't read path for maps: {:?}", err);
+                        return None;
+                    }
+                };
+
+                let reader = BufReader::new(file);
+                let map : WorldFile = match bincode::deserialize_from(reader) {
+                    Ok(map) => map,
+                    Err(err) => {
+                        log::warn!("Couldn't parse map: {:?})", err);
+                        return None;
+                    }
+                };
+
+                if map.alt.len() != map.basement.len() || map.alt.len() != WORLD_SIZE.x as usize * WORLD_SIZE.y as usize {
+                    log::warn!("World size of map is invalid.");
+                    return None;
+                }
+
+                Some(map)
+            } else {
+                None
+            }
+        })();
+
+        let (alt, basement) = if let Some(map) = parsed_world_file {
+            (map.alt, map.basement)
+        } else {
+            let (alt, basement) = do_erosion(
+                0.0,
+                max_erosion_per_delta_t as f32,
+                n_steps,
+                &river_seed,
+                &rock_strength_nz,
+                |posi| alt_func(posi),// + if is_ocean_fn(posi) { 0.0 } else { 128.0 },
+                |posi| alt_func(posi) - if is_ocean_fn(posi) { 0.0 } else { /*1400.0*//*CONFIG.mountain_scale * 0.75*/0.0 },// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
+                is_ocean_fn,
+                uplift_fn,
+                |posi| n_func(posi),
+                |posi| theta_func(posi),
+                |posi| kf_func(posi),
+                |posi| kd_func(posi),
+                |posi| g_func(posi),
+            );
+
+            // Quick "small scale" erosion cycle in order to lower extreme angles.
+            do_erosion(
+                0.0,
+                (1.0 * height_scale) as f32,
+                n_small_steps,
+                &river_seed,
+                &rock_strength_nz,
+                |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi] as f32,
+                |posi| basement[posi] as f32,
+                is_ocean_fn,
+                |posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
+                |posi| n_func(posi),
+                |posi| theta_func(posi),
+                |posi| kf_func(posi),
+                |posi| kd_func(posi),
+                |posi| g_func(posi),
+            )
+        };
+
+        // Save map, if necessary.
+        let map = WorldFile {
+            alt,
+            basement,
+        };
+        (|| {
+            if let FileOpts::Save = opts.world_file {
+                use std::{time::SystemTime};
+                // Check if folder exists and create it if it does not
+                let mut path = PathBuf::from("./maps");
+                if !path.exists() {
+                    if let Err(err) = std::fs::create_dir(&path) {
+                        log::warn!("Couldn't create folder for map: {:?}", err);
+                        return;
+                    }
+                }
+                path.push(format!(
+                    // TODO: Work out a nice bincode file extension.
+                    "map_{}.bin",
+                    SystemTime::now()
+                        .duration_since(SystemTime::UNIX_EPOCH)
+                        .map(|d| d.as_millis())
+                        .unwrap_or(0)
+                ));
+                let file = match File::create(path) {
+                    Ok(file) => file,
+                    Err(err) => {
+                        log::warn!("Couldn't create file for maps: {:?}", err);
+                        return;
+                    }
+                };
+
+                let writer = BufWriter::new(file);
+                if let Err(err) = bincode::serialize_into(writer, &map) {
+                    log::warn!("Couldn't write map: {:?}", err);
+                }
+            }
+        })();
+        let (alt, basement) = (map.alt, map.basement);
+
+        // Additional small-scale eroson after map load, only used during testing.
+        let (alt, basement) = if n_post_load_steps == 0 {
+            (alt, basement)
+        } else {
+            do_erosion(
+                0.0,
+                (1.0 * height_scale) as f32,
+                n_post_load_steps,
+                &river_seed,
+                &rock_strength_nz,
+                |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi] as f32,
+                |posi| basement[posi] as f32,
+                is_ocean_fn,
+                |posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
+                |posi| n_func(posi),
+                |posi| theta_func(posi),
+                |posi| kf_func(posi),
+                |posi| kd_func(posi),
+                |posi| g_func(posi),
+            )
+        };
 
         let is_ocean = get_oceans(|posi| alt[posi] as f32);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
diff --git a/world/src/sim/util.rs b/world/src/sim/util.rs
index ced1fc93fa..698787c8df 100644
--- a/world/src/sim/util.rs
+++ b/world/src/sim/util.rs
@@ -287,7 +287,7 @@ pub fn get_oceans(oldh: impl Fn(usize) -> f32 + Sync) -> BitBox {
         if oldh(posi) <= 0.0 {
             stack.push(posi);
         } else {
-            panic!("Hopefully no border tiles are above sea level.");
+            // panic!("Hopefully no border tiles are above sea level.");
         }
     };
     for x in 0..WORLD_SIZE.x as i32 {