From 314134459b762f438ad0189f5fac4334776de4ff Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 21 Aug 2026 15:04:22 +0200 Subject: [PATCH 1/5] Add regression test for the LOOP_SNOW column shifts The two vertical grid shifts in LOOP_SNOW have never been covered by a test. Add one that pins down both: snowfall_and_deposition must move every interior layer down one index and drop the deepest layer, and melt_sublimation must move every interior layer up one index and add a fresh layer at the base. Also assert the two budgets that follow from that: column mass changes by exactly the snowfall in minus the deepest layer out, and the water credited to runoff_irr_deep equals the water that really left the column. The test column deliberately gives every layer a distinct thickness, density and water content. This matters: the model's own initial state gives all deep layers the same values (INIT sets subD = C["Dice"] throughout and subZ[34:] = 0.8), and duplicating or dropping a layer that is identical to its neighbour is numerically invisible. A test built on a uniform deep column passes whether the shifts are right or wrong. subW is used as the tracer because nothing but the shift can touch it here: the column sits at the melting point so it has no cold content and no water can refreeze, no rain or melt adds water, and every layer's water stays below the mliqmax cap applied by compaction. That makes the permutation assertions exact rather than approximate. The upward shift is driven by sublimation rather than melt for the same reason - both feed the same shift_tot, but sublimation produces no liquid water. These tests fail on the current implementation, which leaves one interior index unwritten in each shift. The fix follows separately. Assisted-by: Claude Code:claude-opus-5. --- tests/core/test_loop_snow_layer_shift.py | 253 +++++++++++++++++++++++ 1 file changed, 253 insertions(+) create mode 100644 tests/core/test_loop_snow_layer_shift.py diff --git a/tests/core/test_loop_snow_layer_shift.py b/tests/core/test_loop_snow_layer_shift.py new file mode 100644 index 00000000..6e7f6eeb --- /dev/null +++ b/tests/core/test_loop_snow_layer_shift.py @@ -0,0 +1,253 @@ +# SPDX-FileCopyrightText: 2026 EBFM Authors +# +# SPDX-License-Identifier: BSD-3-Clause +# +# Content partially generated with the assistance of AI tools. +# Claude Code: Opus 5 + +"""The vertical grid shifts in LOOP_SNOW. + +Two blocks move the subsurface column by one layer index: snowfall_and_deposition +pushes it down when a new top layer is started, and melt_sublimation pulls it up +when the top layer is gone. Both must move *every* interior layer, and the layer +that drops off the bottom is the one snowfall_and_deposition credits to +runoff_irr_deep as water leaving the domain. + +These tests use a column whose per-layer values are all distinct, and in +particular whose two deepest layers differ. That is essential: the model's own +initial state gives every deep layer the same thickness, temperature, density and +water content (INIT sets subD = C["Dice"] throughout and subZ[34:] = 0.8), and +duplicating or dropping a layer identical to its neighbour is numerically +invisible. A test built on a uniform deep column passes whether the shifts are +right or wrong. + +subW is the tracer of choice, because in this setup nothing but the shift can +touch it: the column sits exactly at the melting point so it has no cold content +and no water can refreeze, no rain or melt adds water, and every layer's water +stays below the mliqmax cap that compaction applies. That makes the assertions on +subW exact rather than approximate. +""" + +import unittest + +import numpy as np + +GPSUM, NL = 2, 50 +# Layer thickness (m). The top layer is a full max_subZ, so any snowfall starts a +# new one and forces a downward shift. +MAX_SUBZ = 0.1 +THICKNESS = 0.5 +# Distinct per layer, and inside [C["Dfirn"], C["Dice"]) so only firn compaction +# applies. The 3.0 kg m-3 spacing is what makes a misplaced layer detectable. +DENSITY_STEP = 3.0 +DENSITY_BASE = 600.0 +# Distinct per layer, and far below the mliqmax of these layers (9.6 kg m-2 at the +# base rising to 13.5 at the top), so compaction's cap leaves them untouched. +WATER_BASE, WATER_STEP = 1.0, 0.1 + + +def _make_case(C): + """A column set at the melting point whose every layer differs from its neighbours.""" + subZ = np.full(NL, THICKNESS) + subZ[0] = MAX_SUBZ + subD = DENSITY_BASE + DENSITY_STEP * np.arange(NL) + subW = WATER_BASE + WATER_STEP * np.arange(NL) + + grid = { + "gpsum": GPSUM, + "nl": NL, + "max_subZ": MAX_SUBZ, + "doubledepth": False, + "split": np.array([15]), + # Layer merging/splitting does its own shifts and its own runoff_irr_deep + # credit; switch it off so these tests see the snowfall/melt shifts alone. + "mask": np.zeros(GPSUM, dtype=int), + } + OUT = { + "subZ": np.tile(subZ, (GPSUM, 1)), + "subD": np.tile(subD, (GPSUM, 1)), + "subW": np.tile(subW, (GPSUM, 1)), + "subT": np.full((GPSUM, NL), C["T0"]), + "subTmean": np.full((GPSUM, NL), C["T0"]), + "subS": np.zeros((GPSUM, NL)), + "surfH": np.zeros(GPSUM), + "Tsurf": np.full(GPSUM, C["T0"]), + "melt": np.zeros(GPSUM), + "moist_deposition": np.zeros(GPSUM), + "moist_sublimation": np.zeros(GPSUM), + "moist_condensation": np.zeros(GPSUM), + "moist_evaporation": np.zeros(GPSUM), + "runoff_irr_deep_mean": np.zeros(GPSUM), + } + yearsnow = np.full(GPSUM, 500.0) + IN = { + "T": np.full(GPSUM, C["T0"]), + "WS": np.full(GPSUM, 4.0), + "snow": np.zeros(GPSUM), + "rain": np.zeros(GPSUM), + "yearsnow": np.tile(yearsnow[:, None], (1, NL)), + "logyearsnow": np.tile(np.log(yearsnow)[:, None], (1, NL)), + } + return grid, OUT, IN + + +class TestColumnShift(unittest.TestCase): + def setUp(self): + from ebfm.core import INIT, LOOP_SNOW + + self.LOOP_SNOW = LOOP_SNOW + self.C = INIT.init_constants() + self.phys = {"snow_compaction": "firn+snow", "percolation": "normal"} + + def _step(self, grid, OUT, IN): + self.LOOP_SNOW.main(self.C, OUT, IN, 1.0 / 24.0, grid, self.phys) + + def _snowfall_case(self): + """A case that performs exactly one downward shift. + + The fresh snow is sized from the fresh-snow density the model will compute, so + that it overfills the top layer without filling a second one: one pass of the + shift loop, whose result is a plain permutation of the column. + """ + grid, OUT, IN = _make_case(self.C) + Dfresh = 50 + 1.7 * 17 ** (3 / 2) + 266.86 * (0.5 * (1 + np.tanh(4.0 / 5))) ** 8.8 + IN["snow"][:] = 0.9 * MAX_SUBZ * Dfresh / self.C["Dwater"] + return grid, OUT, IN + + def _sublimation_case(self): + """A case that performs exactly one upward shift. + + Sublimation, not melt, removes the top layer: it takes away exactly as much mass + as the top layer holds, and unlike melt it produces no liquid water, so subW + stays an exact tracer. Both drive the same block through the same `shift_tot`. + """ + grid, OUT, IN = _make_case(self.C) + OUT["moist_sublimation"][:] = OUT["subD"][:, 0] * OUT["subZ"][:, 0] * 1e-3 + return grid, OUT, IN + + def test_snowfall_shifts_every_layer_down(self): + """Every layer moves down one index; the deepest layer leaves the domain. + + The layer that must not survive is the deepest one: if it is still there + afterwards, the column kept a layer that the runoff bookkeeping already + counted as gone. + """ + grid, OUT, IN = self._snowfall_case() + subW_old = OUT["subW"].copy() + + self._step(grid, OUT, IN) + + want = np.empty((GPSUM, NL)) + want[:, 0] = 0.0 # new top layer is fresh snow, no water + want[:, 1] = subW_old[:, 0] # former top layer, handed down whole + want[:, 2:NL] = subW_old[:, 1 : NL - 1] # everything else moves down one index + np.testing.assert_array_equal( + OUT["subW"], + want, + err_msg="subW is not the column shifted down by one index", + ) + + def test_sublimation_shifts_every_layer_up(self): + """Every layer moves up one index; a fresh, dry layer appears at the base.""" + grid, OUT, IN = self._sublimation_case() + subW_old = OUT["subW"].copy() + + self._step(grid, OUT, IN) + + want = np.empty((GPSUM, NL)) + want[:, 0] = subW_old[:, 1] # former layer 1 becomes the whole top layer + want[:, 1 : NL - 1] = subW_old[:, 2:NL] # everything else moves up one index + want[:, NL - 1] = 0.0 # fresh layer added at the base + np.testing.assert_allclose( + OUT["subW"], + want, + rtol=1e-12, + atol=0.0, + err_msg="subW is not the column shifted up by one index", + ) + + def test_snowfall_shifts_density_down(self): + """The density profile rides the same shift as the water profile. + + Firn compaction perturbs subD by ~1e-3 kg m-3 over one hourly step, three + orders of magnitude below the layer-to-layer spacing, so a misplaced layer is + unambiguous at this tolerance. + """ + grid, OUT, IN = self._snowfall_case() + subD_old = OUT["subD"].copy() + + self._step(grid, OUT, IN) + + np.testing.assert_allclose( + OUT["subD"][:, 2:NL], + subD_old[:, 1 : NL - 1], + rtol=1e-4, + err_msg="subD is not the column shifted down by one index", + ) + + def test_sublimation_shifts_density_up(self): + """The density profile rides the same shift as the water profile.""" + grid, OUT, IN = self._sublimation_case() + subD_old = OUT["subD"].copy() + + self._step(grid, OUT, IN) + + np.testing.assert_allclose( + OUT["subD"][:, 1 : NL - 1], + subD_old[:, 2:NL], + rtol=1e-4, + err_msg="subD is not the column shifted up by one index", + ) + + def test_snowfall_shift_conserves_column_mass(self): + """Column mass changes by exactly the snow added minus the layer that left. + + Nothing else in this configuration moves mass: compaction preserves subD*subZ + by construction, heat conduction touches neither, and with no rain, no melt and + no cold content there is no refreezing. So the whole residual is the shift's. + """ + grid, OUT, IN = self._snowfall_case() + mass_before = np.sum(OUT["subD"] * OUT["subZ"], axis=1) + # The fresh snow mass is the snowfall in water equivalent, whatever density the + # model gives it: shift_snowfall * Dfreshsnow == snow * Dwater. + mass_added = IN["snow"] * self.C["Dwater"] + mass_left = OUT["subD"][:, NL - 1] * OUT["subZ"][:, NL - 1] + + self._step(grid, OUT, IN) + + mass_after = np.sum(OUT["subD"] * OUT["subZ"], axis=1) + np.testing.assert_allclose( + mass_after - mass_before, + mass_added - mass_left, + atol=1e-6, + err_msg="column mass across the shift does not balance snowfall in minus deepest layer out", + ) + + def test_water_credited_to_deep_runoff_actually_leaves(self): + """runoff_irr_deep must be credited the water that really left the column. + + This is the bookkeeping the shift has to agree with: the block credits the + deepest layer's water as having left the domain, so the deepest layer's water + is what the shift has to discard. + """ + grid, OUT, IN = self._snowfall_case() + water_before = np.sum(OUT["subW"], axis=1) + + self._step(grid, OUT, IN) + + water_after = np.sum(OUT["subW"], axis=1) + # runoff() has already rescaled runoff_irr_deep into the running mean, so read + # the credit back off the mean, which started at zero. + weight = (1.0 / 24.0) / self.C["yeardays"] + credited = OUT["runoff_irr_deep_mean"] / weight + + np.testing.assert_allclose( + credited, + water_before - water_after, + rtol=1e-12, + err_msg="water credited to runoff_irr_deep differs from the water that left the column", + ) + + +if __name__ == "__main__": + unittest.main() From 5d97a9d36be248245869b3733662b11f2aeb82b9 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 21 Aug 2026 16:10:41 +0200 Subject: [PATCH 2/5] Fix off-by-one in the LOOP_SNOW column shifts Both blocks that shift the subsurface column by one index left exactly one interior index unwritten, so it silently kept its pre-shift value (indices zero-based). The downward shift in snowfall_and_deposition never assigned nl-1, dropping the layer at nl-2 and keeping the one at nl-1 that the runoff_irr_deep credit had already counted as gone. The upward shift in melt_sublimation never assigned nl-2, duplicating that layer and losing the one at nl-1. Widen both slices by one so every interior layer moves. The credit against subW_old[:, nl-1] is deliberately left alone: it was already right for the intended shift, and it was the shift that disagreed with it. layer_merging_and_splitting performs the same two shifts and was audited as correct. Per shift event the mass budget was off by subZ[nl-1] * (subD[nl-1] - subD[nl-2]), gaining mass on snowfall and losing it on melt, and the deep runoff credit by subW[nl-1] - subW[nl-2]. A 10-year single-column run at 500 mm w.e. per year drifted 8 % in column mass and was still diverging. The error vanishes when the two deepest layers are identical, which is how it survived undetected: the model's initial state makes them so. Both blocks are fixed together because their errors pushed thickness in opposite directions, so fixing one alone could expose a partly cancelling error. Assisted-by: Claude Code:claude-opus-5. --- src/ebfm/core/LOOP_SNOW.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/ebfm/core/LOOP_SNOW.py b/src/ebfm/core/LOOP_SNOW.py index d6cee1cd..3638e476 100644 --- a/src/ebfm/core/LOOP_SNOW.py +++ b/src/ebfm/core/LOOP_SNOW.py @@ -142,10 +142,12 @@ def snowfall_and_deposition(): subD_old = OUT["subD"][idx_shift] subW_old = OUT["subW"][idx_shift] - OUT["subZ"][idx_shift, 2 : nl - 1] = subZ_old[:, 1 : nl - 2] - OUT["subT"][idx_shift, 2 : nl - 1] = subT_old[:, 1 : nl - 2] - OUT["subD"][idx_shift, 2 : nl - 1] = subD_old[:, 1 : nl - 2] - OUT["subW"][idx_shift, 2 : nl - 1] = subW_old[:, 1 : nl - 2] + # Shift the column down by one: the layers at indices 1..nl-2 move to + # indices 2..nl-1. The layer at index nl-1 falls out of the domain. + OUT["subZ"][idx_shift, 2:nl] = subZ_old[:, 1 : nl - 1] + OUT["subT"][idx_shift, 2:nl] = subT_old[:, 1 : nl - 1] + OUT["subD"][idx_shift, 2:nl] = subD_old[:, 1 : nl - 1] + OUT["subW"][idx_shift, 2:nl] = subW_old[:, 1 : nl - 1] OUT["subZ"][idx_shift, 1] = max_subZ OUT["subZ"][idx_shift, 0] = (subZ_old[:, 0] + shift_amount[idx_shift]) - max_subZ @@ -163,7 +165,8 @@ def snowfall_and_deposition(): OUT["subW"][idx_shift, 1] = subW_old[:, 0] OUT["subW"][idx_shift, 0] = 0.0 - # Update runoff for shifted layers + # The layer that just fell out of the domain, index nl-1 before the + # shift, takes its irreducible water with it. OUT["runoff_irr_deep"][idx_shift] += subW_old[:, nl - 1] return _SUCCESS @@ -204,6 +207,8 @@ def melt_sublimation(): OUT["surfH"] += shift_amount + nl = grid["nl"] + # No-shift branch only touches top layer, only column 0 copied (instead of whole grid). subZ_top_old = OUT["subZ"][:, 0].copy() subT_top_old = OUT["subT"][:, 0].copy() @@ -233,10 +238,12 @@ def melt_sublimation(): subD_old = OUT["subD"][idx_shift] subW_old = OUT["subW"][idx_shift] - OUT["subZ"][idx_shift, 1 : nl - 2] = subZ_old[:, 2 : nl - 1] - OUT["subT"][idx_shift, 1 : nl - 2] = subT_old[:, 2 : nl - 1] - OUT["subD"][idx_shift, 1 : nl - 2] = subD_old[:, 2 : nl - 1] - OUT["subW"][idx_shift, 1 : nl - 2] = subW_old[:, 2 : nl - 1] + # Shift the column up by one: the layers at indices 2..nl-1 move to + # indices 1..nl-2. + OUT["subZ"][idx_shift, 1 : nl - 1] = subZ_old[:, 2:nl] + OUT["subT"][idx_shift, 1 : nl - 1] = subT_old[:, 2:nl] + OUT["subD"][idx_shift, 1 : nl - 1] = subD_old[:, 2:nl] + OUT["subW"][idx_shift, 1 : nl - 1] = subW_old[:, 2:nl] subZ_top_new = subZ_old[:, 0] + subZ_old[:, 1] + shift_amount[idx_shift] OUT["subZ"][idx_shift, 0] = subZ_top_new @@ -244,6 +251,9 @@ def melt_sublimation(): OUT["subD"][idx_shift, 0] = subD_old[:, 1] # Liquid water scales with the remaining fraction of the former layer 1 OUT["subW"][idx_shift, 0] = subW_old[:, 1] * (subZ_top_new / subZ_old[:, 1]) + # Index nl-1 becomes a fresh layer at the base (vectorized over the + # shift rows): it inherits the temperature and density of the pre-shift + # deepest layer, starts dry, and gets the full deep-layer thickness. OUT["subT"][idx_shift, nl - 1] = subT_old[:, nl - 1] OUT["subW"][idx_shift, nl - 1] = 0.0 @@ -262,6 +272,7 @@ def compaction(gpu=None): Calculate snow and firn compaction and update density and layer thickness """ + gpsum, nl = grid["gpsum"], grid["nl"] Dice, Dfirn = C["Dice"], C["Dfirn"] dt_yearfrac = dt / C["yeardays"] @@ -486,6 +497,7 @@ def heat_conduction(gpu=None): if get_backend() == ComputeBackend.NUMBA: # Numba parallel path # per-column precompute + CFL solve is done in _heat_conduction_prep_kernel + gpsum, nl = OUT["subT"].shape kk = np.empty((gpsum, nl)) c_eff = np.empty((gpsum, nl)) kk_sz_top = np.empty(gpsum) @@ -633,6 +645,8 @@ def percolation_refreezing_and_storage(gpu=None): ######################################################### # Percolation, refreezing and irreducible water storage ######################################################### + gpsum, nl = OUT["subT"].shape + backend = get_backend() if backend in (ComputeBackend.NUMBA, ComputeBackend.GPU): # Numba / GPU kernel path: From 704299fc40555ead73c18196e95290387120336f Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 21 Aug 2026 18:31:21 +0200 Subject: [PATCH 3/5] Add changelog. --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 13065fc9..f2743429 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ SPDX-FileCopyrightText: 2025 EBFM Authors SPDX-License-Identifier: CC-BY-4.0 --> +# develop + +* Fixed off-by-one error in grid shifts that can lead to inaccuracies for long-running simulations. https://github.com/EBFMorg/EBFM/pull/164 + # v0.7.0 * Moved performance and profiling sections out of `README.md` into `docs/Performance.md`, together with the GPU setup and CUDA runtime troubleshooting entry. https://github.com/EBFMorg/EBFM/pull/161 From f905d5b49090f8c2d48b6deb119cdfddcb0d964b Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 24 Aug 2026 09:09:41 +0200 Subject: [PATCH 4/5] Fix merge. --- src/ebfm/core/LOOP_SNOW.py | 6 ------ tests/core/test_loop_snow_layer_shift.py | 5 ++++- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/ebfm/core/LOOP_SNOW.py b/src/ebfm/core/LOOP_SNOW.py index 3638e476..82377d3e 100644 --- a/src/ebfm/core/LOOP_SNOW.py +++ b/src/ebfm/core/LOOP_SNOW.py @@ -207,8 +207,6 @@ def melt_sublimation(): OUT["surfH"] += shift_amount - nl = grid["nl"] - # No-shift branch only touches top layer, only column 0 copied (instead of whole grid). subZ_top_old = OUT["subZ"][:, 0].copy() subT_top_old = OUT["subT"][:, 0].copy() @@ -272,7 +270,6 @@ def compaction(gpu=None): Calculate snow and firn compaction and update density and layer thickness """ - gpsum, nl = grid["gpsum"], grid["nl"] Dice, Dfirn = C["Dice"], C["Dfirn"] dt_yearfrac = dt / C["yeardays"] @@ -497,7 +494,6 @@ def heat_conduction(gpu=None): if get_backend() == ComputeBackend.NUMBA: # Numba parallel path # per-column precompute + CFL solve is done in _heat_conduction_prep_kernel - gpsum, nl = OUT["subT"].shape kk = np.empty((gpsum, nl)) c_eff = np.empty((gpsum, nl)) kk_sz_top = np.empty(gpsum) @@ -645,8 +641,6 @@ def percolation_refreezing_and_storage(gpu=None): ######################################################### # Percolation, refreezing and irreducible water storage ######################################################### - gpsum, nl = OUT["subT"].shape - backend = get_backend() if backend in (ComputeBackend.NUMBA, ComputeBackend.GPU): # Numba / GPU kernel path: diff --git a/tests/core/test_loop_snow_layer_shift.py b/tests/core/test_loop_snow_layer_shift.py index 6e7f6eeb..88066053 100644 --- a/tests/core/test_loop_snow_layer_shift.py +++ b/tests/core/test_loop_snow_layer_shift.py @@ -32,6 +32,8 @@ import numpy as np +from ebfm.core.config import ColumnDiscretizationConfig + GPSUM, NL = 2, 50 # Layer thickness (m). The top layer is a full max_subZ, so any snowfall starts a # new one and forces a downward shift. @@ -98,9 +100,10 @@ def setUp(self): self.LOOP_SNOW = LOOP_SNOW self.C = INIT.init_constants() self.phys = {"snow_compaction": "firn+snow", "percolation": "normal"} + self.column = ColumnDiscretizationConfig() def _step(self, grid, OUT, IN): - self.LOOP_SNOW.main(self.C, OUT, IN, 1.0 / 24.0, grid, self.phys) + self.LOOP_SNOW.main(self.C, OUT, IN, 1.0 / 24.0, grid, self.phys, self.column) def _snowfall_case(self): """A case that performs exactly one downward shift. From eda14a684eaf6b57dd0cec43a1033ea4b3b1f9ee Mon Sep 17 00:00:00 2001 From: wardvp <63638856+wardvp@users.noreply.github.com> Date: Fri, 30 Oct 2026 09:35:04 +0100 Subject: [PATCH 5/5] New inline comments describing the layer shifting --- src/ebfm/core/LOOP_SNOW.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/ebfm/core/LOOP_SNOW.py b/src/ebfm/core/LOOP_SNOW.py index d443983a..1b34e80f 100644 --- a/src/ebfm/core/LOOP_SNOW.py +++ b/src/ebfm/core/LOOP_SNOW.py @@ -170,8 +170,8 @@ def snowfall_and_deposition(): OUT["subW"][idx_shift, 1] = subW_old[:, 0] OUT["subW"][idx_shift, 0] = 0.0 - # The layer that just fell out of the domain, index nl-1 before the - # shift, takes its irreducible water with it. + # The irreducible water in the previously deepest layer (now removed) is added to runoff_irr_deep for + # mass conservation. OUT["runoff_irr_deep"][idx_shift] += subW_old[:, nl - 1] return _SUCCESS @@ -254,13 +254,11 @@ def melt_sublimation(): OUT["subD"][idx_shift, 0] = subD_old[:, 1] # Liquid water scales with the remaining fraction of the former layer 1 OUT["subW"][idx_shift, 0] = subW_old[:, 1] * (subZ_top_new / subZ_old[:, 1]) - # Index nl-1 becomes a fresh layer at the base (vectorized over the - # shift rows): it inherits the temperature and density of the pre-shift - # deepest layer, starts dry, and gets the full deep-layer thickness. + + # Index nl-1 becomes a fresh layer at the base: it inherits the temperature and density of the + # pre-shift deepest layer, starts dry, and gets the full deep-layer thickness. OUT["subT"][idx_shift, nl - 1] = subT_old[:, nl - 1] OUT["subW"][idx_shift, nl - 1] = 0.0 - - # Update the deepest layer properties (vectorized over shift rows) if column.doubledepth: OUT["subZ"][idx_shift, nl - 1] = 2.0 ** len(column.split) * column.max_subZ else: