diff --git a/client/src/lib.rs b/client/src/lib.rs
index 4264dceb80..7a0b633dc9 100644
--- a/client/src/lib.rs
+++ b/client/src/lib.rs
@@ -113,11 +113,8 @@ impl Client {
                 log::debug!("Preparing image...");
                 let world_map = Arc::new(image::DynamicImage::ImageRgba8({
                     // Should not fail if the dimensions are correct.
-                    let world_map = image::ImageBuffer::from_raw(
-                        map_size.x,
-                        map_size.y,
-                        world_map_raw,
-                    );
+                    let world_map =
+                        image::ImageBuffer::from_raw(map_size.x, map_size.y, world_map_raw);
                     world_map.ok_or(Error::Other("Server sent a bad world map image".into()))?
                 }));
                 log::debug!("Done preparing image...");
diff --git a/common/src/msg/server.rs b/common/src/msg/server.rs
index 8636ffa25e..6e3f9ad661 100644
--- a/common/src/msg/server.rs
+++ b/common/src/msg/server.rs
@@ -37,7 +37,7 @@ pub enum ServerMsg {
         entity_package: sync::EntityPackage<EcsCompPacket>,
         server_info: ServerInfo,
         time_of_day: state::TimeOfDay,
-        world_map: (Vec2<u32>, Vec<u32>)
+        world_map: (Vec2<u32>, Vec<u32>),
     },
     PlayerListUpdate(PlayerListUpdate),
     StateAnswer(Result<ClientState, (RequestStateError, ClientState)>),
diff --git a/server/src/lib.rs b/server/src/lib.rs
index 2c42102011..2b04f88b62 100644
--- a/server/src/lib.rs
+++ b/server/src/lib.rs
@@ -46,7 +46,10 @@ use std::{
 };
 use uvth::{ThreadPool, ThreadPoolBuilder};
 use vek::*;
-use world::{sim::{FileOpts, WORLD_SIZE, WorldOpts}, World};
+use world::{
+    sim::{FileOpts, WorldOpts, WORLD_SIZE},
+    World,
+};
 const CLIENT_TIMEOUT: f64 = 20.0; // Seconds
 
 pub enum Event {
@@ -104,15 +107,18 @@ impl Server {
         state.ecs_mut().register::<RegionSubscription>();
         state.ecs_mut().register::<Client>();
 
-        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
+        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()
             },
-            ..WorldOpts::default()
-        });
+        );
 
         let spawn_point = {
             // NOTE: all of these `.map(|e| e as [type])` calls should compile into no-ops,
diff --git a/server/src/settings.rs b/server/src/settings.rs
index 134fc10e50..92d6f7409b 100644
--- a/server/src/settings.rs
+++ b/server/src/settings.rs
@@ -104,8 +104,7 @@ impl ServerSettings {
             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.
+            ..Self::load()                            // Fill in remaining fields from settings.ron.
         }
     }
 
diff --git a/voxygen/src/hud/map.rs b/voxygen/src/hud/map.rs
index fdc11fdf42..e1f4c613e7 100644
--- a/voxygen/src/hud/map.rs
+++ b/voxygen/src/hud/map.rs
@@ -162,7 +162,7 @@ impl<'a> Widget for Map<'a> {
         let (world_map, worldsize) = self.world_map;
         let worldsize = worldsize.map2(TerrainChunkSize::RECT_SIZE, |e, f| e as f64 * f as f64);
 
-        Image::new(world_map/*self.imgs.map_placeholder*/)
+        Image::new(world_map /*self.imgs.map_placeholder*/)
             .middle_of(state.ids.map_bg)
             .color(Some(Color::Rgba(1.0, 1.0, 1.0, fade - 0.1)))
             .w_h(700.0, 700.0)
diff --git a/voxygen/src/hud/minimap.rs b/voxygen/src/hud/minimap.rs
index 2bbf2702b2..a5597661a9 100644
--- a/voxygen/src/hud/minimap.rs
+++ b/voxygen/src/hud/minimap.rs
@@ -137,7 +137,7 @@ impl<'a> Widget for MiniMap<'a> {
             // Map Image
             let (world_map, worldsize) = self.world_map;
             let worldsize = worldsize.map2(TerrainChunkSize::RECT_SIZE, |e, f| e as f64 * f as f64);
-            Image::new(world_map/*self.imgs.map_placeholder*/)
+            Image::new(world_map /*self.imgs.map_placeholder*/)
                 .middle_of(state.ids.mmap_frame_bg)
                 .w_h(92.0 * 4.0 * zoom, 82.0 * 4.0 * zoom)
                 .parent(state.ids.mmap_frame_bg)
diff --git a/voxygen/src/hud/mod.rs b/voxygen/src/hud/mod.rs
index c6b04c2bb0..2bf5e07b25 100644
--- a/voxygen/src/hud/mod.rs
+++ b/voxygen/src/hud/mod.rs
@@ -454,7 +454,10 @@ impl Hud {
         // Generate ids.
         let ids = Ids::new(ui.id_generator());
         // Load world map
-        let world_map = (ui.add_graphic(Graphic::Image(client.world_map.0.clone())), client.world_map.1);
+        let world_map = (
+            ui.add_graphic(Graphic::Image(client.world_map.0.clone())),
+            client.world_map.1,
+        );
         // Load images.
         let imgs = Imgs::load(&mut ui).expect("Failed to load images!");
         // Load rotation images.
diff --git a/voxygen/src/scene/mod.rs b/voxygen/src/scene/mod.rs
index f0691d18af..4faa0f7bf7 100644
--- a/voxygen/src/scene/mod.rs
+++ b/voxygen/src/scene/mod.rs
@@ -20,7 +20,7 @@ use client::Client;
 use common::{
     comp,
     terrain::{BlockKind, TerrainChunk, TerrainChunkSize},
-    vol::{RectVolSize, ReadVol},
+    vol::{ReadVol, RectVolSize},
 };
 use specs::{Join, WorldExt};
 use vek::*;
diff --git a/voxygen/src/scene/terrain.rs b/voxygen/src/scene/terrain.rs
index 5726f2bc28..536714f650 100644
--- a/voxygen/src/scene/terrain.rs
+++ b/voxygen/src/scene/terrain.rs
@@ -175,7 +175,9 @@ fn mesh_worker<V: BaseVol<Vox = Block> + RectRasterableVol + ReadVol + Debug>(
                         let kind = volume.get(wpos).unwrap_or(&Block::empty()).kind();
 
                         if let Some(cfg) = sprite_config_for(kind) {
-                            let seed = wpos.x as u64 * 3 + wpos.y as u64 * 7 + wpos.x as u64 * wpos.y as u64; // Awful PRNG
+                            let seed = wpos.x as u64 * 3
+                                + wpos.y as u64 * 7
+                                + wpos.x as u64 * wpos.y as u64; // Awful PRNG
 
                             let instance = SpriteInstance::new(
                                 Mat4::identity()
diff --git a/world/examples/view.rs b/world/examples/view.rs
index eab8ef1398..6ff94ee124 100644
--- a/world/examples/view.rs
+++ b/world/examples/view.rs
@@ -6,10 +6,13 @@ const W: usize = 640;
 const H: usize = 480;
 
 fn main() {
-    let world = World::generate(0, WorldOpts {
-        seed_elements: true,
-        ..WorldOpts::default()
-    });
+    let world = World::generate(
+        0,
+        WorldOpts {
+            seed_elements: true,
+            ..WorldOpts::default()
+        },
+    );
 
     let sampler = world.sample_columns();
 
diff --git a/world/examples/water.rs b/world/examples/water.rs
index 7b3c977ebe..68a7d48aa3 100644
--- a/world/examples/water.rs
+++ b/world/examples/water.rs
@@ -36,12 +36,15 @@ fn main() {
     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()
-    });
+    let world = World::generate(
+        1337,
+        WorldOpts {
+            seed_elements: false,
+            // world_file: sim::FileOpts::Load(_map_file),
+            world_file: sim::FileOpts::Save,
+            ..WorldOpts::default()
+        },
+    );
 
     let sampler = world.sim();
 
@@ -87,25 +90,35 @@ fn main() {
 
         for i in 0..W {
             for j in 0..H {
-                let pos = (focus_rect + Vec2::new(i as f64, j as f64) * scale).map(|e: f64| e as i32);
+                let pos =
+                    (focus_rect + Vec2::new(i as f64, j as f64) * scale).map(|e: f64| e as i32);
                 /* let top_left = pos;
                 let top_right = focus + Vec2::new(i as i32 + light_res, j as i32) * scale;
                 let bottom_left = focus + Vec2::new(i as i32, j as i32 + light_res) * scale; */
 
-                let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) = sampler
-                    .get(pos)
-                    .map(|sample| {
-                        (
-                            sample.alt,
-                            sample.basement,
-                            sample.water_alt,
-                            sample.humidity,
-                            sample.temp,
-                            sample.downhill,
-                            sample.river.river_kind,
-                        )
-                    })
-                    .unwrap_or((CONFIG.sea_level, CONFIG.sea_level, CONFIG.sea_level, 0.0, 0.0, None, None));
+                let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) =
+                    sampler
+                        .get(pos)
+                        .map(|sample| {
+                            (
+                                sample.alt,
+                                sample.basement,
+                                sample.water_alt,
+                                sample.humidity,
+                                sample.temp,
+                                sample.downhill,
+                                sample.river.river_kind,
+                            )
+                        })
+                        .unwrap_or((
+                            CONFIG.sea_level,
+                            CONFIG.sea_level,
+                            CONFIG.sea_level,
+                            0.0,
+                            0.0,
+                            None,
+                            None,
+                        ));
                 let humidity = humidity.min(1.0).max(0.0);
                 let temperature = temperature.min(1.0).max(-1.0) * 0.5 + 0.5;
                 let pos = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32);
@@ -177,15 +190,9 @@ fn main() {
 
                 let true_water_alt = (alt.max(water_alt) as f64 - focus.z) / gain as f64;
                 let true_alt = (alt as f64 - focus.z) / gain as f64;
-                let water_depth = (true_water_alt - true_alt)
-                    .min(1.0)
-                    .max(0.0);
-                let water_alt = true_water_alt
-                    .min(1.0)
-                    .max(0.0);
-                let alt = true_alt
-                    .min(1.0)
-                    .max(0.0);
+                let water_depth = (true_water_alt - true_alt).min(1.0).max(0.0);
+                let water_alt = true_water_alt.min(1.0).max(0.0);
+                let alt = true_alt.min(1.0).max(0.0);
                 let quad =
                     |x: f32| ((x as f64 * QUADRANTS as f64).floor() as usize).min(QUADRANTS - 1);
                 if river_kind.is_none() || humidity != 0.0 {
@@ -205,17 +212,29 @@ fn main() {
                 }
 
                 buf[j * W + i] = match (river_kind, (is_water, true_alt >= true_sea_level)) {
-                    (_, (false, _)) | ( None, (_, true)) => {
+                    (_, (false, _)) | (None, (_, true)) => {
                         let (r, g, b) = (
-                            (if is_shaded { alt } else { alt } * if is_temperature { temperature as f64 } else if is_shaded { alt } else { 0.0 }).sqrt(),
+                            (if is_shaded { alt } else { alt }
+                                * if is_temperature {
+                                    temperature as f64
+                                } else if is_shaded {
+                                    alt
+                                } else {
+                                    0.0
+                                })
+                            .sqrt(),
                             if is_shaded { 0.2 + (alt * 0.8) } else { alt },
-                            (if is_shaded { alt } else { alt } * if is_humidity { humidity as f64 } else if is_shaded { alt } else { 0.0 }).sqrt(),
+                            (if is_shaded { alt } else { alt }
+                                * if is_humidity {
+                                    humidity as f64
+                                } else if is_shaded {
+                                    alt
+                                } else {
+                                    0.0
+                                })
+                            .sqrt(),
                         );
-                        let light = if is_shaded {
-                            light
-                        } else {
-                            1.0
-                        };
+                        let light = if is_shaded { light } else { 1.0 };
                         u32::from_le_bytes([
                             (b * light * 255.0) as u8,
                             (g * light * 255.0) as u8,
@@ -228,7 +247,7 @@ fn main() {
                             (/*alt*//*alt * *//*(1.0 - humidity)*/(alt * temperature).sqrt() * 255.0) as u8,
                             255,
                         ]) */
-                    },
+                    }
                     (Some(RiverKind::Ocean), _) => u32::from_le_bytes([
                         ((64.0 - water_depth * 64.0) * 1.0) as u8,
                         ((32.0 - water_depth * 32.0) * 1.0) as u8,
@@ -242,8 +261,8 @@ fn main() {
                         255,
                     ]),
                     (None, _) | (Some(RiverKind::Lake { .. }), _) => u32::from_le_bytes([
-                        (((64.0 + water_alt * 191.0) + (- water_depth * 64.0)) * 1.0) as u8,
-                        (((32.0 + water_alt * 95.0) + (- water_depth * 32.0)) * 1.0) as u8,
+                        (((64.0 + water_alt * 191.0) + (-water_depth * 64.0)) * 1.0) as u8,
+                        (((32.0 + water_alt * 95.0) + (-water_depth * 32.0)) * 1.0) as u8,
                         0,
                         255,
                     ]),
@@ -280,7 +299,8 @@ fn main() {
         }
         if win.get_mouse_down(minifb::MouseButton::Left) {
             if let Some((mx, my)) = win.get_mouse_pos(minifb::MouseMode::Clamp) {
-                let pos = (focus_rect + (Vec2::new(mx as f64, my as f64) * scale)).map(|e| e as i32);
+                let pos =
+                    (focus_rect + (Vec2::new(mx as f64, my as f64) * scale)).map(|e| e as i32);
                 println!(
                     "Chunk position: {:?}",
                     pos.map2(TerrainChunkSize::RECT_SIZE, |e, f| e * f as i32)
diff --git a/world/src/block/mod.rs b/world/src/block/mod.rs
index 0aefcd7b72..83d6553ff7 100644
--- a/world/src/block/mod.rs
+++ b/world/src/block/mod.rs
@@ -265,7 +265,8 @@ impl<'a> BlockGen<'a> {
                     sub_surface_color,
                     stone_col.map(|e| e as f32 / 255.0),
                     (height - grass_depth - wposf.z as f32) * 0.15,
-                ).map(|e| (e * 255.0) as u8);
+                )
+                .map(|e| (e * 255.0) as u8);
 
                 // Underground
                 if (wposf.z as f32) > alt - 32.0 * chaos {
diff --git a/world/src/column/mod.rs b/world/src/column/mod.rs
index bb006f148e..0a2f859572 100644
--- a/world/src/column/mod.rs
+++ b/world/src/column/mod.rs
@@ -437,13 +437,14 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
          downhill_alt.sub(CONFIG.sea_level) >= CONFIG.mountain_scale * 0.25)*/
         is_rocky {
             sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
-            // sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
-            // sim.get_interpolated(wpos, |chunk| chunk.alt)?
+        // sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
+        // sim.get_interpolated(wpos, |chunk| chunk.alt)?
         } else {
             sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
             // sim.get_interpolated(wpos, |chunk| chunk.alt)?
         };
-        let basement = alt + sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?;
+        let basement = alt
+            + sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?;
 
         // Find the average distance to each neighboring body of water.
         let mut river_count = 0.0f64;
@@ -830,7 +831,6 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 .mul(1.0 - humidity) */
                 /* .mul(32.0) */;
 
-
         let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor);
         let alt = alt_ + riverless_alt_delta;
         let basement = basement.min(alt);
@@ -895,7 +895,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         );
         let tundra = Lerp::lerp(snow, Rgb::new(0.01, 0.3, 0.0), 0.4 + marble * 0.6);
         let dead_tundra = Lerp::lerp(warm_stone, Rgb::new(0.3, 0.12, 0.2), marble);
-        let cliff = Rgb::lerp(cold_stone, /*warm_stone*/hot_stone, marble);
+        let cliff = Rgb::lerp(cold_stone, /*warm_stone*/ hot_stone, marble);
 
         let grass = Rgb::lerp(
             cold_grass,
@@ -953,11 +953,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 .mul(1.0),
         );
 
-        let sub_surface_color = Lerp::lerp(
-            cliff,
-            ground,
-            alt.sub(basement).mul(0.25)
-        );
+        let sub_surface_color = Lerp::lerp(cliff, ground, alt.sub(basement).mul(0.25));
 
         /*  let ground = Rgb::lerp(
             dead_tundra,
@@ -1087,15 +1083,18 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         ); */
 
         // Snow covering
-        let snow_cover =
-            temp.sub(CONFIG.snow_temp)
-                .max(-humidity.sub(CONFIG.desert_hum))
-                .mul(16.0)
-                .add((marble_small - 0.5) * 0.5);
+        let snow_cover = temp
+            .sub(CONFIG.snow_temp)
+            .max(-humidity.sub(CONFIG.desert_hum))
+            .mul(16.0)
+            .add((marble_small - 0.5) * 0.5);
         let (alt, ground, sub_surface_color) = if snow_cover /*< 0.1*/<= 0.5 && alt > water_level {
             // Allow snow cover.
-            (alt + 1.0 - snow_cover.max(0.0), Rgb::lerp(snow, ground, snow_cover),
-                Lerp::lerp(sub_surface_color, ground, alt.sub(basement).mul(0.15)))
+            (
+                alt + 1.0 - snow_cover.max(0.0),
+                Rgb::lerp(snow, ground, snow_cover),
+                Lerp::lerp(sub_surface_color, ground, alt.sub(basement).mul(0.15)),
+            )
         } else {
             (alt, ground, sub_surface_color)
         };
@@ -1165,7 +1164,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 // Beach
                 ((ocean_level - 1.0) / 2.0).max(0.0),
             ),
-            sub_surface_color,// /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)),
+            sub_surface_color, // /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)),
             // No growing directly on bedrock.
             tree_density: Lerp::lerp(0.0, tree_density, alt.sub(2.0).sub(basement).mul(0.5)),
             forest_kind: sim_chunk.forest_kind,
@@ -1183,7 +1182,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             humidity,
             spawn_rate,
             location: sim_chunk.location.as_ref(),
-	    stone_col,
+            stone_col,
 
             chunk: sim_chunk,
             spawn_rules: sim_chunk
diff --git a/world/src/sim/diffusion.rs b/world/src/sim/diffusion.rs
index eebeb5e229..5e7cc82575 100644
--- a/world/src/sim/diffusion.rs
+++ b/world/src/sim/diffusion.rs
@@ -1,5 +1,5 @@
-use rayon::prelude::*;
 use super::Alt;
+use rayon::prelude::*;
 
 /// From https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
 ///
@@ -36,68 +36,76 @@ use super::Alt;
 
   implicit none
 */
-pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
-                 h: &mut [Alt], b: &mut [Alt],
-                 kd: impl Fn(usize) -> f64, kdsed: f64,
-                 ) {
+pub fn diffusion(
+    nx: usize,
+    ny: usize,
+    xl: f64,
+    yl: f64,
+    dt: f64,
+    _ibc: (),
+    h: &mut [Alt],
+    b: &mut [Alt],
+    kd: impl Fn(usize) -> f64,
+    kdsed: f64,
+) {
     let aij = |i: usize, j: usize| j * nx + i;
-/*
-  double precision, dimension(:), allocatable :: f,diag,sup,inf,res
-  double precision, dimension(:,:), allocatable :: zint,kdint,zintp
-  integer i,j,ij
-  double precision factxp,factxm,factyp,factym,dx,dy
-*/
-    let mut f : Vec<f64>;
-    let mut diag : Vec<f64>;
-    let mut sup : Vec<f64>;
-    let mut inf : Vec<f64>;
-    let mut res : Vec<f64>;
-    let mut zint : Vec<f64>;
-    let mut kdint : Vec<f64>;
-    let mut zintp : Vec<f64>;
-    let mut i : usize;
-    let mut j : usize;
-    let mut ij : usize;
-    let mut factxp : f64;
-    let mut factxm : f64;
-    let mut factyp : f64;
-    let mut factym : f64;
-    let mut dx : f64;
-    let mut dy : f64;
-/*
-  character cbc*4
+    /*
+      double precision, dimension(:), allocatable :: f,diag,sup,inf,res
+      double precision, dimension(:,:), allocatable :: zint,kdint,zintp
+      integer i,j,ij
+      double precision factxp,factxm,factyp,factym,dx,dy
+    */
+    let mut f: Vec<f64>;
+    let mut diag: Vec<f64>;
+    let mut sup: Vec<f64>;
+    let mut inf: Vec<f64>;
+    let mut res: Vec<f64>;
+    let mut zint: Vec<f64>;
+    let mut kdint: Vec<f64>;
+    let mut zintp: Vec<f64>;
+    let mut i: usize;
+    let mut j: usize;
+    let mut ij: usize;
+    let mut factxp: f64;
+    let mut factxm: f64;
+    let mut factyp: f64;
+    let mut factym: f64;
+    let mut dx: f64;
+    let mut dy: f64;
+    /*
+      character cbc*4
 
-  !print*,'Diffusion'
+      !print*,'Diffusion'
 
-  write (cbc,'(i4)') ibc
+      write (cbc,'(i4)') ibc
 
-  dx=xl/(nx-1)
-  dy=yl/(ny-1)
-*/
+      dx=xl/(nx-1)
+      dy=yl/(ny-1)
+    */
     // 2048*32/2048/2048
     // 1 / 64 m
     dx = xl / /*(nx - 1)*/nx as f64;
     dy = yl / /*(ny - 1)*/ny as f64;
-/*
-  ! creates 2D internal arrays to store topo and kd
+    /*
+      ! creates 2D internal arrays to store topo and kd
 
-  allocate (zint(nx,ny),kdint(nx,ny),zintp(nx,ny))
-*/
+      allocate (zint(nx,ny),kdint(nx,ny),zintp(nx,ny))
+    */
     zint = vec![Default::default(); nx * ny];
     kdint = vec![Default::default(); nx * ny];
     zintp = vec![Default::default(); nx * ny];
-/*
-  do j=1,ny
-    do i=1,nx
-      ij=(j-1)*nx+i
-      zint(i,j)=h(ij)
-      kdint(i,j)=kd(ij)
-      if (kdsed.gt.0.d0 .and. (h(ij)-b(ij)).gt.1.d-6) kdint(i,j)=kdsed
-    enddo
-  enddo
+    /*
+      do j=1,ny
+        do i=1,nx
+          ij=(j-1)*nx+i
+          zint(i,j)=h(ij)
+          kdint(i,j)=kd(ij)
+          if (kdsed.gt.0.d0 .and. (h(ij)-b(ij)).gt.1.d-6) kdint(i,j)=kdsed
+        enddo
+      enddo
 
-  zintp = zint
-*/
+      zintp = zint
+    */
     for j in 0..ny {
         for i in 0..nx {
             // ij = vec2_as_uniform_idx(i, j);
@@ -111,63 +119,62 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
     }
 
     zintp = zint.clone();
-/*
-  ! first pass along the x-axis
+    /*
+      ! first pass along the x-axis
 
-  allocate (f(nx),diag(nx),sup(nx),inf(nx),res(nx))
-  f=0.d0
-  diag=0.d0
-  sup=0.d0
-  inf=0.d0
-  res=0.d0
-  do j=2,ny-1
-*/
+      allocate (f(nx),diag(nx),sup(nx),inf(nx),res(nx))
+      f=0.d0
+      diag=0.d0
+      sup=0.d0
+      inf=0.d0
+      res=0.d0
+      do j=2,ny-1
+    */
     f = vec![0.0; nx];
     diag = vec![0.0; nx];
     sup = vec![0.0; nx];
     inf = vec![0.0; nx];
     res = vec![0.0; nx];
-    for j in 1..ny-1 {
-/*
-    do i=2,nx-1
-    factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
-    factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
-      factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
-    factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
-    diag(i)=1.d0+factxp+factxm
-    sup(i)=-factxp
-    inf(i)=-factxm
-    f(i)=zintp(i,j)+factyp*zintp(i,j+1)-(factyp+factym)*zintp(i,j)+factym*zintp(i,j-1)
-    enddo
-*/
-        for i in 1..nx-1 {
-            factxp = (kdint[aij(i+1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxm = (kdint[aij(i-1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factyp = (kdint[aij(i, j+1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
-            factym = (kdint[aij(i, j-1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+    for j in 1..ny - 1 {
+        /*
+            do i=2,nx-1
+            factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+            factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+              factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+            factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+            diag(i)=1.d0+factxp+factxm
+            sup(i)=-factxp
+            inf(i)=-factxm
+            f(i)=zintp(i,j)+factyp*zintp(i,j+1)-(factyp+factym)*zintp(i,j)+factym*zintp(i,j-1)
+            enddo
+        */
+        for i in 1..nx - 1 {
+            factxp = (kdint[aij(i + 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i - 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, j + 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(i, j - 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[i] = 1.0 + factxp + factxm;
             sup[i] = -factxp;
             inf[i] = -factxm;
-            f[i] =
-                zintp[aij(i, j)] + factyp * zintp[aij(i, j+1)] -
-                (factyp + factym) *
-                zintp[aij(i, j)] + factym * zintp[aij(i, j-1)];
+            f[i] = zintp[aij(i, j)] + factyp * zintp[aij(i, j + 1)]
+                - (factyp + factym) * zintp[aij(i, j)]
+                + factym * zintp[aij(i, j - 1)];
         }
-/*
-! left bc
-    if (cbc(4:4).eq.'1') then
-    diag(1)=1.
-    sup(1)=0.
-    f(1)=zintp(1,j)
-    else
-    factxp=(kdint(2,j)+kdint(1,j))/2.d0*(dt/2.)/dx**2
-    factyp=(kdint(1,j+1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
-    factym=(kdint(1,j-1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
-    diag(1)=1.d0+factxp
-    sup(1)=-factxp
-    f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1)
-    endif
-*/
+        /*
+        ! left bc
+            if (cbc(4:4).eq.'1') then
+            diag(1)=1.
+            sup(1)=0.
+            f(1)=zintp(1,j)
+            else
+            factxp=(kdint(2,j)+kdint(1,j))/2.d0*(dt/2.)/dx**2
+            factyp=(kdint(1,j+1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
+            factym=(kdint(1,j-1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
+            diag(1)=1.d0+factxp
+            sup(1)=-factxp
+            f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1)
+            endif
+        */
         if true {
             diag[0] = 1.0;
             sup[0] = 0.0;
@@ -175,119 +182,118 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
         } else {
             // reflective boundary
             factxp = (kdint[aij(1, j)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factyp = (kdint[aij(0, j+1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
-            factym = (kdint[aij(0, j-1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factyp = (kdint[aij(0, j + 1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(0, j - 1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[0] = 1.0 + factxp;
             sup[0] = -factxp;
-            f[0] =
-                zintp[aij(0, j)] + factyp * zintp[aij(0, j+1)] -
-                (factyp + factym) *
-                zintp[aij(0, j)] + factym * zintp[aij(0, j-1)];
+            f[0] = zintp[aij(0, j)] + factyp * zintp[aij(0, j + 1)]
+                - (factyp + factym) * zintp[aij(0, j)]
+                + factym * zintp[aij(0, j - 1)];
         }
-/*
-! right bc
-    if (cbc(2:2).eq.'1') then
-    diag(nx)=1.
-    inf(nx)=0.
-    f(nx)=zintp(nx,j)
-    else
-    factxm=(kdint(nx-1,j)+kdint(nx,j))/2.d0*(dt/2.)/dx**2
-    factyp=(kdint(nx,j+1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
-    factym=(kdint(nx,j-1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
-    diag(nx)=1.d0+factxm
-    inf(nx)=-factxm
-    f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1)
-    endif
-*/
+        /*
+        ! right bc
+            if (cbc(2:2).eq.'1') then
+            diag(nx)=1.
+            inf(nx)=0.
+            f(nx)=zintp(nx,j)
+            else
+            factxm=(kdint(nx-1,j)+kdint(nx,j))/2.d0*(dt/2.)/dx**2
+            factyp=(kdint(nx,j+1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
+            factym=(kdint(nx,j-1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
+            diag(nx)=1.d0+factxm
+            inf(nx)=-factxm
+            f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1)
+            endif
+        */
         if true {
             diag[nx - 1] = 1.0;
             inf[nx - 1] = 0.0;
             f[nx - 1] = zintp[aij(nx - 1, j)];
         } else {
             // reflective boundary
-            factxm = (kdint[aij(nx-2, j)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factyp = (kdint[aij(nx-1, j+1)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
-            factym = (kdint[aij(nx-1, j-1)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factxm = (kdint[aij(nx - 2, j)] + kdint[aij(nx - 1, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp =
+                (kdint[aij(nx - 1, j + 1)] + kdint[aij(nx - 1, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym =
+                (kdint[aij(nx - 1, j - 1)] + kdint[aij(nx - 1, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[nx - 1] = 1.0 + factxm;
             inf[nx - 1] = -factxm;
-            f[nx - 1] =
-                zintp[aij(nx-1, j)] + factyp * zintp[aij(nx-1, j+1)] -
-                (factyp + factym) *
-                zintp[aij(nx-1, j)] + factym * zintp[aij(nx-1, j-1)];
+            f[nx - 1] = zintp[aij(nx - 1, j)] + factyp * zintp[aij(nx - 1, j + 1)]
+                - (factyp + factym) * zintp[aij(nx - 1, j)]
+                + factym * zintp[aij(nx - 1, j - 1)];
         }
-/*
-  call tridag (inf,diag,sup,f,res,nx)
-    do i=1,nx
-    zint(i,j)=res(i)
-    enddo
-*/
+        /*
+          call tridag (inf,diag,sup,f,res,nx)
+            do i=1,nx
+            zint(i,j)=res(i)
+            enddo
+        */
         tridag(&inf, &diag, &sup, &f, &mut res, nx);
         for i in 0..nx {
             zint[aij(i, j)] = res[i];
         }
-/*
-  enddo
-  deallocate (f,diag,sup,inf,res)
-*/
+        /*
+          enddo
+          deallocate (f,diag,sup,inf,res)
+        */
     }
 
-/*
-  ! second pass along y-axis
+    /*
+      ! second pass along y-axis
 
-  allocate (f(ny),diag(ny),sup(ny),inf(ny),res(ny))
-  f=0.d0
-  diag=0.d0
-  sup=0.d0
-  inf=0.d0
-  res=0.d0
-  do i=2,nx-1
-*/
+      allocate (f(ny),diag(ny),sup(ny),inf(ny),res(ny))
+      f=0.d0
+      diag=0.d0
+      sup=0.d0
+      inf=0.d0
+      res=0.d0
+      do i=2,nx-1
+    */
     f = vec![0.0; ny];
     diag = vec![0.0; ny];
     sup = vec![0.0; ny];
     inf = vec![0.0; ny];
     res = vec![0.0; ny];
-    for i in 1..nx-1 {
-/*
-    do j=2,ny-1
-    factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
-    factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
-    factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
-    factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
-    diag(j)=1.d0+factyp+factym
-    sup(j)=-factyp
-    inf(j)=-factym
-    f(j)=zint(i,j)+factxp*zint(i+1,j)-(factxp+factxm)*zint(i,j)+factxm*zint(i-1,j)
-    enddo
-*/
-        for j in 1..ny-1 {
-            factxp = (kdint[aij(i+1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxm = (kdint[aij(i-1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factyp = (kdint[aij(i, j+1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
-            factym = (kdint[aij(i, j-1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+    for i in 1..nx - 1 {
+        /*
+            do j=2,ny-1
+            factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+            factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+            factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+            factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+            diag(j)=1.d0+factyp+factym
+            sup(j)=-factyp
+            inf(j)=-factym
+            f(j)=zint(i,j)+factxp*zint(i+1,j)-(factxp+factxm)*zint(i,j)+factxm*zint(i-1,j)
+            enddo
+        */
+        for j in 1..ny - 1 {
+            factxp = (kdint[aij(i + 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i - 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, j + 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(i, j - 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[j] = 1.0 + factyp + factym;
             sup[j] = -factyp;
             inf[j] = -factym;
-            f[j] =
-                zint[aij(i, j)] + factxp * zint[aij(i+1, j)] -
-                (factxp + factxm) *
-                zint[aij(i, j)] + factxm * zint[aij(i-1, j)];
+            f[j] = zint[aij(i, j)] + factxp * zint[aij(i + 1, j)]
+                - (factxp + factxm) * zint[aij(i, j)]
+                + factxm * zint[aij(i - 1, j)];
         }
-/*
-! bottom bc
-    if (cbc(1:1).eq.'1') then
-    diag(1)=1.
-    sup(1)=0.
-    f(1)=zint(i,1)
-    else
-    factxp=(kdint(i+1,1)+kdint(i,j))/2.d0*(dt/2.)/dx**2
-    factxm=(kdint(i-1,1)+kdint(i,1))/2.d0*(dt/2.)/dx**2
-    factyp=(kdint(i,2)+kdint(i,1))/2.d0*(dt/2.)/dy**2
-    diag(1)=1.d0+factyp
-    sup(1)=-factyp
-    f(1)=zint(i,1)+factxp*zint(i+1,1)-(factxp+factxm)*zint(i,1)+factxm*zint(i-1,1)
-    endif
-*/
+        /*
+        ! bottom bc
+            if (cbc(1:1).eq.'1') then
+            diag(1)=1.
+            sup(1)=0.
+            f(1)=zint(i,1)
+            else
+            factxp=(kdint(i+1,1)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+            factxm=(kdint(i-1,1)+kdint(i,1))/2.d0*(dt/2.)/dx**2
+            factyp=(kdint(i,2)+kdint(i,1))/2.d0*(dt/2.)/dy**2
+            diag(1)=1.d0+factyp
+            sup(1)=-factyp
+            f(1)=zint(i,1)+factxp*zint(i+1,1)-(factxp+factxm)*zint(i,1)+factxm*zint(i-1,1)
+            endif
+        */
         if true {
             diag[0] = 1.0;
             sup[0] = 0.0;
@@ -297,76 +303,76 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
             // TODO: Check whether this j was actually supposed to be a 0 in the original
             // (probably).
             // factxp = (kdint[aij(i+1, 0)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxp = (kdint[aij(i+1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxm = (kdint[aij(i-1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxp = (kdint[aij(i + 1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i - 1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx);
             factyp = (kdint[aij(i, 1)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[0] = 1.0 + factyp;
             sup[0] = -factyp;
-            f[0] =
-                zint[aij(i, 0)] + factxp * zint[aij(i+1, 0)] -
-                (factxp + factxm) *
-                zint[aij(i, 0)] + factxm * zint[aij(i-1, 0)];
+            f[0] = zint[aij(i, 0)] + factxp * zint[aij(i + 1, 0)]
+                - (factxp + factxm) * zint[aij(i, 0)]
+                + factxm * zint[aij(i - 1, 0)];
         }
-/*
-! top bc
-    if (cbc(3:3).eq.'1') then
-    diag(ny)=1.
-    inf(ny)=0.
-    f(ny)=zint(i,ny)
-    else
-    factxp=(kdint(i+1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
-    factxm=(kdint(i-1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
-    factym=(kdint(i,ny-1)+kdint(i,ny))/2.d0*(dt/2.)/dy**2
-    diag(ny)=1.d0+factym
-    inf(ny)=-factym
-    f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny)
-    endif
-*/
+        /*
+        ! top bc
+            if (cbc(3:3).eq.'1') then
+            diag(ny)=1.
+            inf(ny)=0.
+            f(ny)=zint(i,ny)
+            else
+            factxp=(kdint(i+1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
+            factxm=(kdint(i-1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
+            factym=(kdint(i,ny-1)+kdint(i,ny))/2.d0*(dt/2.)/dy**2
+            diag(ny)=1.d0+factym
+            inf(ny)=-factym
+            f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny)
+            endif
+        */
         if true {
             diag[ny - 1] = 1.0;
             inf[ny - 1] = 0.0;
             f[ny - 1] = zint[aij(i, ny - 1)];
         } else {
             // reflective boundary
-            factxp = (kdint[aij(i+1, ny-1)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxm = (kdint[aij(i-1, ny-1)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factym = (kdint[aij(i, ny-2)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factxp =
+                (kdint[aij(i + 1, ny - 1)] + kdint[aij(i, ny - 1)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm =
+                (kdint[aij(i - 1, ny - 1)] + kdint[aij(i, ny - 1)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factym = (kdint[aij(i, ny - 2)] + kdint[aij(i, ny - 1)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[ny - 1] = 1.0 + factym;
             inf[ny - 1] = -factym;
-            f[ny - 1] =
-                zint[aij(i, ny-1)] + factxp * zint[aij(i+1, ny-1)] -
-                (factxp + factxm) *
-                zint[aij(i, ny-1)] + factxm * zint[aij(i-1, ny-1)];
+            f[ny - 1] = zint[aij(i, ny - 1)] + factxp * zint[aij(i + 1, ny - 1)]
+                - (factxp + factxm) * zint[aij(i, ny - 1)]
+                + factxm * zint[aij(i - 1, ny - 1)];
         }
-/*
-  call tridag (inf,diag,sup,f,res,ny)
-    do j=1,ny
-    zintp(i,j)=res(j)
-    enddo
-*/
+        /*
+          call tridag (inf,diag,sup,f,res,ny)
+            do j=1,ny
+            zintp(i,j)=res(j)
+            enddo
+        */
         tridag(&inf, &diag, &sup, &f, &mut res, ny);
         for j in 0..ny {
             zintp[aij(i, j)] = res[j];
         }
-/*
-  enddo
-  deallocate (f,diag,sup,inf,res)
-*/
+        /*
+          enddo
+          deallocate (f,diag,sup,inf,res)
+        */
     }
-/*
-  ! stores result in 1D array
+    /*
+      ! stores result in 1D array
 
-  do j=1,ny
-    do i=1,nx
-      ij=(j-1)*nx+i
-      etot(ij)=etot(ij)+h(ij)-zintp(i,j)
-      erate(ij)=erate(ij)+(h(ij)-zintp(i,j))/dt
-      h(ij)=zintp(i,j)
-    enddo
-  enddo
+      do j=1,ny
+        do i=1,nx
+          ij=(j-1)*nx+i
+          etot(ij)=etot(ij)+h(ij)-zintp(i,j)
+          erate(ij)=erate(ij)+(h(ij)-zintp(i,j))/dt
+          h(ij)=zintp(i,j)
+        enddo
+      enddo
 
-  b=min(h,b)
-*/
+      b=min(h,b)
+    */
     for j in 0..ny {
         for i in 0..nx {
             ij = aij(i, j);
@@ -378,13 +384,13 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
     b.par_iter_mut().zip(h).for_each(|(mut b, h)| {
         *b = h.min(*b);
     });
-/*
-  deallocate (zint,kdint,zintp)
+    /*
+      deallocate (zint,kdint,zintp)
 
-  return
+      return
 
-end subroutine Diffusion
-*/
+    end subroutine Diffusion
+    */
 }
 
 /*
@@ -400,39 +406,39 @@ end subroutine Diffusion
       double precision a(n),b(n),c(n),r(n),u(n)
 */
 pub fn tridag(a: &[f64], b: &[f64], c: &[f64], r: &[f64], u: &mut [f64], n: usize) {
-/*
-      INTEGER j
-      double precision bet
-      double precision,dimension(:),allocatable::gam
+    /*
+          INTEGER j
+          double precision bet
+          double precision,dimension(:),allocatable::gam
 
-      allocate (gam(n))
+          allocate (gam(n))
 
-      if(b(1).eq.0.d0) stop 'in tridag'
-*/
-    let mut j : usize;
-    let mut bet : f64;
-    let mut precision : f64;
+          if(b(1).eq.0.d0) stop 'in tridag'
+    */
+    let mut j: usize;
+    let mut bet: f64;
+    let mut precision: f64;
     let mut gam: Vec<f64>;
 
     gam = vec![Default::default(); n];
 
     assert!(b[0] != 0.0);
-/*
+    /*
 
-! first pass
+    ! first pass
 
-bet=b(1)
-u(1)=r(1)/bet
-do 11 j=2,n
-  gam(j)=c(j-1)/bet
-  bet=b(j)-a(j)*gam(j)
-  if(bet.eq.0.) then
-    print*,'tridag failed'
-    stop
-  endif
-  u(j)=(r(j)-a(j)*u(j-1))/bet
-  11    continue
-*/
+    bet=b(1)
+    u(1)=r(1)/bet
+    do 11 j=2,n
+      gam(j)=c(j-1)/bet
+      bet=b(j)-a(j)*gam(j)
+      if(bet.eq.0.) then
+        print*,'tridag failed'
+        stop
+      endif
+      u(j)=(r(j)-a(j)*u(j-1))/bet
+      11    continue
+    */
     bet = b[0];
     u[0] = r[0] / bet;
     for j in 1..n {
@@ -448,21 +454,21 @@ do 11 j=2,n
         //               = r'[j] / b'[j]
         u[j] = (r[j] - a[j] * u[j - 1]) / bet;
     }
-/*
-  ! second pass
+    /*
+      ! second pass
 
-  do 12 j=n-1,1,-1
-    u(j)=u(j)-gam(j+1)*u(j+1)
-    12    continue
-*/
+      do 12 j=n-1,1,-1
+        u(j)=u(j)-gam(j+1)*u(j+1)
+        12    continue
+    */
     for j in (0..n - 1).rev() {
         u[j] = u[j] - gam[j + 1] * u[j + 1];
     }
-/*
-    deallocate (gam)
+    /*
+        deallocate (gam)
 
-    return
+        return
 
-    END
-*/
+        END
+    */
 }
diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs
index 9fddd9bd25..b8fb96acff 100644
--- a/world/src/sim/mod.rs
+++ b/world/src/sim/mod.rs
@@ -12,9 +12,9 @@ pub use self::erosion::{
 pub use self::location::Location;
 pub use self::settlement::Settlement;
 pub use self::util::{
-    cdf_irwin_hall, downhill, get_oceans, HybridMulti as HybridMulti_, local_cells, map_edge_factor, neighbors,
-    NEIGHBOR_DELTA, ScaleBias,
-    uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, InverseCdf,
+    cdf_irwin_hall, downhill, get_oceans, local_cells, map_edge_factor, neighbors,
+    uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, HybridMulti as HybridMulti_,
+    InverseCdf, ScaleBias, NEIGHBOR_DELTA,
 };
 
 use crate::{
@@ -29,8 +29,8 @@ use common::{
     vol::RectVolSize,
 };
 use noise::{
-    BasicMulti, Billow, Fbm, HybridMulti, MultiFractal, NoiseFn, RangeFunction,
-    RidgedMulti, Seedable, SuperSimplex, Worley,
+    BasicMulti, Billow, Fbm, HybridMulti, MultiFractal, NoiseFn, RangeFunction, RidgedMulti,
+    Seedable, SuperSimplex, Worley,
 };
 use num::{Float, Signed};
 use rand::{Rng, SeedableRng};
@@ -39,9 +39,9 @@ use rayon::prelude::*;
 use serde_derive::{Deserialize, Serialize};
 use std::{
     collections::HashMap,
-    io::{BufReader, BufWriter},
     f32, f64,
     fs::File,
+    io::{BufReader, BufWriter},
     ops::{Add, Div, Mul, Neg, Sub},
     path::PathBuf,
     sync::Arc,
@@ -54,10 +54,7 @@ use vek::*;
 // cleanly representable in f32 (that stops around 1024 * 4 * 1024 * 4, for signed floats anyway)
 // but I think that is probably less important since I don't think we actually cast a chunk id to
 // float, just coordinates... could be wrong though!
-pub const WORLD_SIZE: Vec2<usize> = Vec2 {
-    x: 1024,
-    y: 1024,
-};
+pub const WORLD_SIZE: Vec2<usize> = Vec2 { x: 1024, y: 1024 };
 
 /// A structure that holds cached noise values and cumulative distribution functions for the input
 /// that led to those values.  See the definition of InverseCdf for a description of how to
@@ -107,7 +104,7 @@ pub(crate) struct GenCtx {
     pub town_gen: StructureGen2d,
 
     pub river_seed: RandomField,
-    pub rock_strength_nz: /*HybridMulti_*/Fbm,
+    pub rock_strength_nz: Fbm,
     pub uplift_nz: Worley,
 }
 
@@ -146,7 +143,7 @@ impl Default for WorldOpts {
 /// 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)]
+#[derive(Serialize, Deserialize)]
 pub struct WorldFile {
     /// Saved altitude height map.
     pub alt: Box<[Alt]>,
@@ -166,7 +163,9 @@ pub struct WorldSim {
 impl WorldSim {
     pub fn generate(seed: u32, opts: WorldOpts) -> Self {
         let mut rng = ChaChaRng::from_seed(seed_expan::rng_state(seed));
-        let continent_scale = 5_000.0f64/*32768.0*/.div(32.0).mul(TerrainChunkSize::RECT_SIZE.x as f64);
+        let continent_scale = 5_000.0f64 /*32768.0*/
+            .div(32.0)
+            .mul(TerrainChunkSize::RECT_SIZE.x as f64);
         let rock_lacunarity = /*0.5*/2.0/*HybridMulti::DEFAULT_LACUNARITY*/;
         let uplift_scale = /*512.0*//*256.0*/128.0;
         let uplift_turb_scale = uplift_scale / 4.0/*32.0*//*64.0*/;
@@ -290,9 +289,9 @@ impl WorldSim {
         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 = /*200*//*10_000*//*1000*//*50*//*100*/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 = 0;//25;//8;//50;//50;//8;//8;//8;//8;//8; // 8
-        let n_post_load_steps = 0;//25;//8
+        let n_steps = /*200*//*10_000*//*1000*//*50*//*100*/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 = 0; //25;//8;//50;//50;//8;//8;//8;//8;//8; // 8
+        let n_post_load_steps = 0; //25;//8
 
         // Logistic regression.  Make sure x ∈ (0, 1).
         let logit = |x: f64| x.ln() - (-x).ln_1p();
@@ -304,7 +303,8 @@ impl WorldSim {
 
         let exp_inverse_cdf = |x: f64/*, pow: f64*/| -(-x).ln_1p()/* / ln(pow)*/;
         // 2 / pi * ln(tan(pi/2 * p))
-        let hypsec_inverse_cdf = |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln());
+        let hypsec_inverse_cdf =
+            |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln());
 
         let min_epsilon =
             1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64).max(f64::EPSILON as f64 * 0.5);
@@ -364,10 +364,9 @@ impl WorldSim {
                             .alt_nz
                             .get((wposf.div(10_000.0)).into_array())
                             .min(1.0)
-                            .max(-1.0)
-                            /* .mul(0.25)
-                            .add(0.125) */)
-                            // .add(0.5)
+                            .max(-1.0)/* .mul(0.25)
+                        .add(0.125) */)
+                        // .add(0.5)
                         .sub(0.05)
                         // .add(0.05)
                         // .add(0.075)
@@ -525,102 +524,110 @@ impl WorldSim {
         });
 
         // Calculate oceans.
-        let old_height = |posi: usize| alt_old[posi].1 * CONFIG.mountain_scale * height_scale as f32;
+        let old_height =
+            |posi: usize| alt_old[posi].1 * CONFIG.mountain_scale * height_scale as f32;
         /* let is_ocean = (0..WORLD_SIZE.x * WORLD_SIZE.y)
-            .into_par_iter()
-            .map(|i| map_edge_factor(i) == 0.0)
-            .collect::<Vec<_>>(); */
+        .into_par_iter()
+        .map(|i| map_edge_factor(i) == 0.0)
+        .collect::<Vec<_>>(); */
         let is_ocean = get_oceans(old_height);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
         let turb_wposf_div = 8.0/*64.0*/;
 
-        let uplift_nz_dist = gen_ctx.uplift_nz
-            .clone()
-            .enable_range(true);
+        let uplift_nz_dist = gen_ctx.uplift_nz.clone().enable_range(true);
         // Recalculate altitudes without oceans.
         // 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(turb_wposf_div);
-                    let turb = Vec2::new(
-                        gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
-                        gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                    ) * uplift_turb_scale * 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 uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
+        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(turb_wposf_div);
+                let turb = Vec2::new(
+                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                ) * uplift_turb_scale
+                    * 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 uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
                         .min(1.0)
                         .max(-1.0)
                         .mul(0.5)
                         .add(0.5); */
                         chaos[posi].1;
 
-                    let uchaos_1 = (uchaos as f64) / 1.32;
+                let uchaos_1 = (uchaos as f64) / 1.32;
 
-                    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 udist_3 = (1.0 - udist).max(0.0);
-                    let udist_4 = udist.min(1.0);
-                    let variation = 1.0.min(64.0 * 64.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64 * (TerrainChunkSize::RECT_SIZE.x as f64 * TerrainChunkSize::RECT_SIZE.y as f64 / 128.0 / 128.0)));
-                    let variation_1 = (uheight * /*udist_2*/udist_4).min(variation);
-                    let height =
-                        (oheight + 0.5).powf(2.0);
-                        // 1.0 - variation + variation * uchaos_1;
-                        // uheight * /*udist_2*/udist_4 - variation_1 + variation_1 * uchaos_1;
-                        // uheight * (0.5 + 0.5 * ((uchaos as f64) / 1.32)) - 0.125;
-                        // 0.2;
-                        // 1.0;
-                        // uheight_1;
-                        // uheight_1 * (0.8 + 0.2 * oheight.signum() * oheight.abs().powf(0.25));
-                        // 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)
-                }
-            });
+                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 udist_3 = (1.0 - udist).max(0.0);
+                let udist_4 = udist.min(1.0);
+                let variation = 1.0.min(
+                    64.0 * 64.0
+                        / (WORLD_SIZE.x as f64
+                            * WORLD_SIZE.y as f64
+                            * (TerrainChunkSize::RECT_SIZE.x as f64
+                                * TerrainChunkSize::RECT_SIZE.y as f64
+                                / 128.0
+                                / 128.0)),
+                );
+                let variation_1 = (uheight * /*udist_2*/udist_4).min(variation);
+                let height = (oheight + 0.5).powf(2.0);
+                // 1.0 - variation + variation * uchaos_1;
+                // uheight * /*udist_2*/udist_4 - variation_1 + variation_1 * uchaos_1;
+                // uheight * (0.5 + 0.5 * ((uchaos as f64) / 1.32)) - 0.125;
+                // 0.2;
+                // 1.0;
+                // uheight_1;
+                // uheight_1 * (0.8 + 0.2 * oheight.signum() * oheight.abs().powf(0.25));
+                // 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)
+            }
+        });
 
         let old_height_uniform = |posi: usize| alt_old_no_ocean[posi].0;
         let alt_old_min_uniform = 0.0;
@@ -684,22 +691,25 @@ impl WorldSim {
             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(turb_wposf_div);
+            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(turb_wposf_div);
             let turb = Vec2::new(
                 gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
                 gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-            ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+            ) * uplift_turb_scale
+                * 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())
+            let uheight = gen_ctx
+                .uplift_nz
+                .get(turb_wposf.into_array())
                 /* .min(0.5)
                 .max(-0.5)*/
                 .min(1.0)
@@ -716,9 +726,7 @@ impl WorldSim {
             // ((3.5 - 1.5) * (1.0 - uheight) + 1.5) as f32
             1.0
         };
-        let theta_func = |posi| {
-            0.4
-        };
+        let theta_func = |posi| 0.4;
         let kf_func = {
             |posi| {
                 if is_ocean_fn(posi) {
@@ -731,19 +739,23 @@ impl WorldSim {
                 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(turb_wposf_div);
+                let turb_wposf = wposf
+                    .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
+                    .div(turb_wposf_div);
                 let turb = Vec2::new(
                     gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
                     gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                ) * uplift_turb_scale
+                    * 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())
+                let uheight = gen_ctx
+                    .uplift_nz
+                    .get(turb_wposf.into_array())
                     /* .min(0.5)
                     .max(-0.5)*/
                     .min(1.0)
@@ -784,7 +796,7 @@ impl WorldSim {
                 // ((1.0 - uheight) * (0.5 - 0.5 * /*((1.32 - uchaos as f64) / 1.32)*/oheight_2) * (1.5e-4 - 2.0e-6) + 2.0e-6)
                 // ((1.0 - uheight) * (0.5 - 0.5 * /*((1.32 - uchaos as f64) / 1.32)*/oheight) * (1.5e-4 - 2.0e-6) + 2.0e-6)
                 // 2e-5
-                2.5e-6 * 16.0.powf(0.4)/* / 4.0 * 0.25 *//* * 4.0*/
+                2.5e-6 * 16.0.powf(0.4) /* / 4.0 * 0.25 *//* * 4.0*/
                 // 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)
@@ -798,19 +810,23 @@ impl WorldSim {
                 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(turb_wposf_div);
+                let turb_wposf = wposf
+                    .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
+                    .div(turb_wposf_div);
                 let turb = Vec2::new(
                     gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
                     gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                ) * uplift_turb_scale
+                    * 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())
+                let uheight = gen_ctx
+                    .uplift_nz
+                    .get(turb_wposf.into_array())
                     /* .min(0.5)
                     .max(-0.5)*/
                     .min(1.0)
@@ -821,209 +837,215 @@ impl WorldSim {
                 // kd = 1.5e-2: normal-high (plateau [fan sediment])
                 // kd = 1e-2: normal (plateau)
                 1.0e-2 * (1.0 / 16.0) // m_old^2 / y * (1 m_new / 4 m_old)^2
-                // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2)
+                                      // (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))
+        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 =
-                    wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div);
-                let turb = Vec2::new(
-                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
-                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                ) * uplift_turb_scale * 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);
+            let turb_wposf = wposf
+                .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
+                .div(turb_wposf_div);
+            let turb = Vec2::new(
+                gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+            ) * uplift_turb_scale
+                * 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);
 
-                let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
+            let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array())
                     .min(1.0)
                     .max(-1.0)
                     .mul(0.5)
                     .add(0.5); */
                     chaos[posi].1;
 
-                assert!(uchaos <= 1.32);
+            assert!(uchaos <= 1.32);
 
-                // G = d* v_s / p_0, where
-                //  v_s is the settling velocity of sediment grains
-                //  p_0 is the mean precipitation rate
-                //  d* is the sediment concentration ratio (between concentration near riverbed
-                //  interface, and average concentration over the water column).
-                //  d* varies with Rouse number which defines relative contribution of bed, suspended,
-                //  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.32].
-                // (((1.32 - uchaos) / 1.32).powf(0.75) * 1.32).min(/*1.1*/1.0)
-                // ((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
-                // 5.0
-                // 10.0
-                // 2.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)
-                    * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+            // G = d* v_s / p_0, where
+            //  v_s is the settling velocity of sediment grains
+            //  p_0 is the mean precipitation rate
+            //  d* is the sediment concentration ratio (between concentration near riverbed
+            //  interface, and average concentration over the water column).
+            //  d* varies with Rouse number which defines relative contribution of bed, suspended,
+            //  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.32].
+            // (((1.32 - uchaos) / 1.32).powf(0.75) * 1.32).min(/*1.1*/1.0)
+            // ((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
+            // 5.0
+            // 10.0
+            // 2.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) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
                 .map(|e| e as f64);
-                let alt_main = {
-                    // Extension upwards from the base.  A positive number from 0 to 1 curved to be
-                    // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
-                    let alt_main = (gen_ctx
-                        .alt_nz
-                        .get((wposf.div(2_000.0)).into_array())
+            let alt_main = {
+                // Extension upwards from the base.  A positive number from 0 to 1 curved to be
+                // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
+                let alt_main = (gen_ctx
+                    .alt_nz
+                    .get((wposf.div(2_000.0)).into_array())
+                    .min(1.0)
+                    .max(-1.0))
+                .abs()
+                .powf(1.35);
+
+                fn spring(x: f64, pow: f64) -> f64 {
+                    x.abs().powf(pow) * x.signum()
+                }
+
+                (0.0 + alt_main
+                    + (gen_ctx
+                        .small_nz
+                        .get((wposf.div(300.0)).into_array())
                         .min(1.0)
                         .max(-1.0))
-                    .abs()
-                    .powf(1.35);
-
-                    fn spring(x: f64, pow: f64) -> f64 {
-                        x.abs().powf(pow) * x.signum()
-                    }
-
-                    (0.0 + alt_main
-                        + (gen_ctx
-                            .small_nz
-                            .get((wposf.div(300.0)).into_array())
-                            .min(1.0)
-                            .max(-1.0))
-                        .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
-                        .mul(0.3)
-                        .add(1.0)
-                        .mul(0.4)
-                        /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
-                            .mul(0.045)*/)
-                };
-                let height =
+                    .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
+                    .mul(0.3)
+                    .add(1.0)
+                    .mul(0.4)/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
+                .mul(0.045)*/)
+            };
+            let height =
                     ((/*old_height_uniform*/uplift_uniform[posi]./*0*/1 - alt_old_min_uniform) as f64
                     / (alt_old_max_uniform - alt_old_min_uniform) as f64)
                     /*((old_height(posi) - alt_old_min) as f64
                     / (alt_old_max - alt_old_min) as f64)*/
                 ;
 
-                let height = height.mul(max_epsilon - min_epsilon).add(min_epsilon);
-                /*.max(1e-7 / CONFIG.mountain_scale as f64)
-                .min(1.0f64 - 1e-7);*/
-                /* let alt_main = {
-                    // Extension upwards from the base.  A positive number from 0 to 1 curved to be
-                    // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
-                    let alt_main = (gen_ctx
-                        .alt_nz
-                        .get((wposf.div(2_000.0)).into_array())
+            let height = height.mul(max_epsilon - min_epsilon).add(min_epsilon);
+            /*.max(1e-7 / CONFIG.mountain_scale as f64)
+            .min(1.0f64 - 1e-7);*/
+            /* let alt_main = {
+                // Extension upwards from the base.  A positive number from 0 to 1 curved to be
+                // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
+                let alt_main = (gen_ctx
+                    .alt_nz
+                    .get((wposf.div(2_000.0)).into_array())
+                    .min(1.0)
+                    .max(-1.0))
+                .abs()
+                .powf(1.35);
+
+                fn spring(x: f64, pow: f64) -> f64 {
+                    x.abs().powf(pow) * x.signum()
+                }
+
+                (0.0 + alt_main
+                    + (gen_ctx
+                        .small_nz
+                        .get((wposf.div(300.0)).into_array())
                         .min(1.0)
                         .max(-1.0))
-                    .abs()
-                    .powf(1.35);
-
-                    fn spring(x: f64, pow: f64) -> f64 {
-                        x.abs().powf(pow) * x.signum()
-                    }
-
-                    (0.0 + alt_main
-                        + (gen_ctx
-                            .small_nz
-                            .get((wposf.div(300.0)).into_array())
-                            .min(1.0)
-                            .max(-1.0))
-                        .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
-                        .mul(0.3)
-                        .add(1.0)
-                        .mul(0.4)
-                        + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0).mul(0.045))
-                }; */
-                // let height = height + (alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64 * ((1.0 / CONFIG.mountain_scale as f64).powf(1.0 / erosion_pow_low));
-                let height = erosion_factor(height);
-                assert!(height >= 0.0);
-                assert!(height <= 1.0);
-                // assert!(alt_main >= 0.0);
-                let (bump_factor, bump_max) = if
-                /*height < f32::EPSILON as f64 * 0.5*//*false*/
-                /*true*/false {
-                    (
-                        /*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/
-                        (alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64,
-                        (f32::EPSILON * 0.5) as f64,
-                    )
-                } else {
-                    (0.0, 0.0)
-                };
-                // tan(6/360*2*pi)*32 ~ 3.4
-                // 3.4/32*512 ~ 54
-                // 18/32*512 ~ 288
-                // tan(pi/6)*32 ~ 18
-                // tan(54/360*2*pi)*32
-                // let height = 1.0f64;
-                let turb_wposf =
-                    wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div);
-                let turb = Vec2::new(
-                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
-                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
-                ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
-                let turb_wposf = wposf + turb;
-                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);
-                // u = 1e-3: normal-high (dike, mountain)
-                // 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;
-
-                // let height = 1.0 / 7.0f64;
-                // let height = 0.0 / 31.0f64;
-                let bfrac = /*erosion_factor(0.5);*/0.0;
-                let height = (height - bfrac).abs().div(1.0 - bfrac);
-                let height = height
-                    /* .mul(31.0 / 32.0)
-                    .add(1.0 / 32.0) */
-                    /* .mul(15.0 / 16.0)
-                    .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)
-                    .sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max)
-                    .add(bump_factor);
-                /* .sub(/*1.0 / CONFIG.mountain_scale as f64*/(f32::EPSILON * 0.5) as f64)
-                .add(bump_factor); */
-                height as f32
+                    .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
+                    .mul(0.3)
+                    .add(1.0)
+                    .mul(0.4)
+                    + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0).mul(0.045))
+            }; */
+            // let height = height + (alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64 * ((1.0 / CONFIG.mountain_scale as f64).powf(1.0 / erosion_pow_low));
+            let height = erosion_factor(height);
+            assert!(height >= 0.0);
+            assert!(height <= 1.0);
+            // assert!(alt_main >= 0.0);
+            let (bump_factor, bump_max) = if
+            /*height < f32::EPSILON as f64 * 0.5*//*false*/
+            /*true*/
+            false {
+                (
+                    /*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/
+                    (alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64,
+                    (f32::EPSILON * 0.5) as f64,
+                )
+            } else {
+                (0.0, 0.0)
             };
+            // tan(6/360*2*pi)*32 ~ 3.4
+            // 3.4/32*512 ~ 54
+            // 18/32*512 ~ 288
+            // tan(pi/6)*32 ~ 18
+            // tan(54/360*2*pi)*32
+            // let height = 1.0f64;
+            let turb_wposf = wposf
+                .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
+                .div(turb_wposf_div);
+            let turb = Vec2::new(
+                gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+            ) * uplift_turb_scale
+                * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+            let turb_wposf = wposf + turb;
+            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);
+            // u = 1e-3: normal-high (dike, mountain)
+            // 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;
+
+            // let height = 1.0 / 7.0f64;
+            // let height = 0.0 / 31.0f64;
+            let bfrac = /*erosion_factor(0.5);*/0.0;
+            let height = (height - bfrac).abs().div(1.0 - bfrac);
+            let height = height
+                /* .mul(31.0 / 32.0)
+                .add(1.0 / 32.0) */
+                /* .mul(15.0 / 16.0)
+                .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)
+                .sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max)
+                .add(bump_factor);
+            /* .sub(/*1.0 / CONFIG.mountain_scale as f64*/(f32::EPSILON * 0.5) as f64)
+            .add(bump_factor); */
+            height as f32
+        };
         let alt_func = |posi| {
             if is_ocean_fn(posi) {
                 // -max_erosion_per_delta_t as f32
@@ -1062,9 +1084,8 @@ impl WorldSim {
                         .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
                         .mul(0.3)
                         .add(1.0)
-                        .mul(0.4)
-                        /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
-                            .mul(0.045)*/)
+                        .mul(0.4)/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
+                    .mul(0.045)*/)
                 };
 
                 // (kf_func(posi) / 1.5e-4 * CONFIG.mountain_scale as f64) as f32
@@ -1112,7 +1133,7 @@ impl WorldSim {
                 };
 
                 let reader = BufReader::new(file);
-                let map : WorldFile = match bincode::deserialize_from(reader) {
+                let map: WorldFile = match bincode::deserialize_from(reader) {
                     Ok(map) => map,
                     Err(err) => {
                         log::warn!("Couldn't parse map: {:?})", err);
@@ -1120,7 +1141,9 @@ impl WorldSim {
                     }
                 };
 
-                if map.alt.len() != map.basement.len() || map.alt.len() != WORLD_SIZE.x as usize * WORLD_SIZE.y as usize {
+                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;
                 }
@@ -1149,8 +1172,16 @@ impl WorldSim {
                 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 },
+                |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),
@@ -1180,13 +1211,10 @@ impl WorldSim {
         };
 
         // Save map, if necessary.
-        let map = WorldFile {
-            alt,
-            basement,
-        };
+        let map = WorldFile { alt, basement };
         (|| {
             if let FileOpts::Save = opts.world_file {
-                use std::{time::SystemTime};
+                use std::time::SystemTime;
                 // Check if folder exists and create it if it does not
                 let mut path = PathBuf::from("./maps");
                 if !path.exists() {
@@ -1243,7 +1271,10 @@ impl WorldSim {
 
         let is_ocean = get_oceans(|posi| alt[posi]);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
-        let mut dh = downhill(|posi| alt[posi] as f32/*&alt*/, /*old_height*/ is_ocean_fn);
+        let mut dh = downhill(
+            |posi| alt[posi] as f32, /*&alt*/
+            /*old_height*/ is_ocean_fn,
+        );
         let (boundary_len, indirection, water_alt_pos, _) = get_lakes(/*&/*water_alt*/alt*/|posi| alt[posi] as f32, &mut dh);
         let flux_old = get_drainage(&water_alt_pos, &dh, boundary_len);
 
@@ -1258,7 +1289,9 @@ impl WorldSim {
             /* // Find the pass this lake is flowing into (i.e. water at the lake bottom gets
             // pushed towards the point identified by pass_idx).
             let neighbor_pass_idx = dh[lake_idx]; */
-            let chunk_water_alt = if /*neighbor_pass_idx*/dh[lake_idx] < 0 {
+            let chunk_water_alt = if
+            /*neighbor_pass_idx*/
+            dh[lake_idx] < 0 {
                 // This is either a boundary node (dh[chunk_idx] == -2, i.e. water is at sea level)
                 // or part of a lake that flows directly into the ocean.  In the former case, water
                 // is at sea level so we just return 0.0.  In the latter case, the lake bottom must
@@ -1293,9 +1326,9 @@ impl WorldSim {
 
         let water_alt = fill_sinks(water_height_initial, is_ocean_fn);
         /* let water_alt = (0..WORLD_SIZE.x * WORLD_SIZE.y)
-            .into_par_iter()
-            .map(|posi| water_height_initial(posi))
-            .collect::<Vec<_>>(); */
+        .into_par_iter()
+        .map(|posi| water_height_initial(posi))
+        .collect::<Vec<_>>(); */
 
         let rivers = get_rivers(&water_alt_pos, &water_alt, &dh, &indirection, &flux_old);
 
@@ -1312,7 +1345,9 @@ impl WorldSim {
                 /* // Find the pass this lake is flowing into (i.e. water at the lake bottom gets
                 // pushed towards the point identified by pass_idx).
                 let neighbor_pass_idx = dh[lake_idx]; */
-                if /*neighbor_pass_idx*/dh[lake_idx] < 0 {
+                if
+                /*neighbor_pass_idx*/
+                dh[lake_idx] < 0 {
                     // This is either a boundary node (dh[chunk_idx] == -2, i.e. water is at sea level)
                     // or part of a lake that flows directly into the ocean.  In the former case, water
                     // is at sea level so we just return 0.0.  In the latter case, the lake bottom must
@@ -1970,9 +2005,7 @@ impl SimChunk {
         let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale;
         let mut alt = CONFIG.sea_level.add(alt_pre.div(height_scale));
         let mut basement = CONFIG.sea_level.add(basement_pre.div(height_scale));
-        let water_alt = CONFIG
-            .sea_level
-            .add(water_alt_pre.div(height_scale));
+        let water_alt = CONFIG.sea_level.add(water_alt_pre.div(height_scale));
         let downhill = if downhill_pre == -2 {
             None
         } else if downhill_pre < 0 {
@@ -2030,7 +2063,9 @@ impl SimChunk {
         }
 
         // No trees in the ocean, with zero humidity (currently), or directly on bedrock.
-        let tree_density = if is_underwater/* || alt - basement.min(alt) < 2.0 */ {
+        let tree_density = if is_underwater
+        /* || alt - basement.min(alt) < 2.0 */
+        {
             0.0
         } else {
             let tree_density = (gen_ctx.tree_nz.get((wposf.div(1024.0)).into_array()))
diff --git a/world/src/sim/util.rs b/world/src/sim/util.rs
index fd587de445..65bf6886de 100644
--- a/world/src/sim/util.rs
+++ b/world/src/sim/util.rs
@@ -218,7 +218,7 @@ pub fn local_cells(posi: usize) -> impl Clone + Iterator<Item = usize> {
 }
 
 // NOTE: want to keep this such that the chunk index is in ascending order!
-pub const NEIGHBOR_DELTA : [(i32, i32); 8] = [
+pub const NEIGHBOR_DELTA: [(i32, i32); 8] = [
     (-1, -1),
     (0, -1),
     (1, -1),
@@ -233,12 +233,12 @@ pub const NEIGHBOR_DELTA : [(i32, i32); 8] = [
 pub fn neighbors(posi: usize) -> impl Clone + Iterator<Item = usize> {
     let pos = uniform_idx_as_vec2(posi);
     NEIGHBOR_DELTA
-    .iter()
-    .map(move |&(x, y)| Vec2::new(pos.x + x, pos.y + y))
-    .filter(|pos| {
-        pos.x >= 0 && pos.y >= 0 && pos.x < WORLD_SIZE.x as i32 && pos.y < WORLD_SIZE.y as i32
-    })
-    .map(vec2_as_uniform_idx)
+        .iter()
+        .map(move |&(x, y)| Vec2::new(pos.x + x, pos.y + y))
+        .filter(|pos| {
+            pos.x >= 0 && pos.y >= 0 && pos.x < WORLD_SIZE.x as i32 && pos.y < WORLD_SIZE.y as i32
+        })
+        .map(vec2_as_uniform_idx)
 }
 
 // Note that we should already have okay cache locality since we have a grid.
@@ -249,10 +249,14 @@ pub fn uphill<'a>(dh: &'a [isize], posi: usize) -> impl Clone + Iterator<Item =
 /// Compute the neighbor "most downhill" from all chunks.
 ///
 /// TODO: See if allocating in advance is worthwhile.
-pub fn downhill<F: Float>(h: impl Fn(usize) -> F + Sync, is_ocean: impl Fn(usize) -> bool + Sync) -> Box<[isize]> {
+pub fn downhill<F: Float>(
+    h: impl Fn(usize) -> F + Sync,
+    is_ocean: impl Fn(usize) -> bool + Sync,
+) -> Box<[isize]> {
     // Constructs not only the list of downhill nodes, but also computes an ordering (visiting
     // nodes in order from roots to leaves).
-    (0..WORLD_SIZE.x * WORLD_SIZE.y).into_par_iter()
+    (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
         // .enumerate()
         .map(|(posi/*, &nh*/)| {
             let nh = h(posi);
@@ -702,5 +706,3 @@ impl<'a, F: NoiseFn<T> + 'a, T> NoiseFn<T> for ScaleBias<'a, F> {
         (self.source.get(point) * self.scale) + self.bias
     }
 }
-
-