Bug fixes

Author: samharrison7

CHANGELOG.md

@@ -4,6 +4,14 @@ All notable changes to the model will be documented in this file. Breaking chang
 
 ## [Unreleased]
 
+## [0.0.7] - 2026-08-19
+
+### Fixed
+
+- Fixed `V_w_b` being carried between iterations in `BedSediment%deposit` (`depositSediment1`). It is assigned only inside the `if (O%A_f(s) > 0 .or. O%A_w(s) > 0)` capacity check but consumed unconditionally when tallying `V_w_tot`, so a layer with no free capacity re-added the previous iteration's value, compounding once per sediment layer and inflating the water volume returned to the water column. It is now reset at the top of each layer iteration. `V_w_tot` only feeds the reporting-only depth adjustment in `depositToBed`, so this makes no difference to model results (verified bit-identical over the Thames scenario), but the previous behaviour was unbounded.
+- `water__depth` (and the `depth(m)` CSV column) now report the hydraulic depth implied by the reach volume, `volume / bedArea`, rather than the `me%depth` attribute. `depositToBed` subtracts the water displaced into the bed sediment from `me%depth`, which left the reported depth inconsistent with the reported `water__volume` and `sediment__bed_area` by ~1e-4 relative.
+- Fixed an error stop in `depositToBed` when `&sediment > include_bed_sediment = .false.`. `V_water_toDeposit = .dp. depositRslt` sat outside the `includeBedSediment` block, so `depositRslt` was never assigned and asking it for its data aborted the run on the first time step with "Error trying to return 0D data as REAL(DP)". The extraction is now inside the block and `V_water_toDeposit` defaults to zero.
+
 ## [0.0.6] - 2026-08-18
 
 ### Fixed
@@ -79,7 +87,8 @@ All notable changes to the model will be documented in this file. Breaking chang
 - Added option to aggregate CSV output for waterbodies at grid cell level, rather than breaking it down to waterbody level. Internal functions for aggregating to grid cell added (e.g. weighted means, fetching outflow reaches). This option can be used by specifying `&output > include_waterbody_breakdown = .false.` in the [model config file](./config.example/config.example.nml). Default is `.true.`.
 - This changelog.
 
-[unreleased]: https://github.com/nerc-ceh/nanofase/compare/0.0.6...HEAD
+[unreleased]: https://github.com/nerc-ceh/nanofase/compare/0.0.7...HEAD
+[0.0.7]: https://github.com/nerc-ceh/nanofase/releases/tag/0.0.7
 [0.0.6]: https://github.com/nerc-ceh/nanofase/releases/tag/0.0.6
 [0.0.5]: https://github.com/nerc-ceh/nanofase/releases/tag/0.0.5
 [0.0.4]: https://github.com/nerc-ceh/nanofase/releases/tag/0.0.4

src/BedSediment/BedSedimentModule.f90

@@ -521,6 +521,11 @@ function depositSediment1(me, FS_dep) result(r)
         do s = 1, me%nSizeClasses                                       ! now add in the depositing sediment, work by size class
             do l = C%nSedimentLayers, 1, -1                             ! start with the bottom layer and work upwards
                 if (FS_dep(s)%M_f() > 0.0_dp) then
+                    V_w_b = 0.0_dp                                      ! reset the receiving layer's water requirement for THIS layer.
+                                                                        ! V_w_b is only assigned inside the capacity check below but is
+                                                                        ! consumed unconditionally when tallying V_w_tot, so without this
+                                                                        ! reset a layer with no free capacity re-adds the previous
+                                                                        ! iteration's value, compounding once per layer
                     associate(O => me%colBedSedimentLayers(l)%item)     ! size class S in Layer L
                         if (O%A_f(s) > 0.0_dp .or. &
                             O%A_w(s) > 0.0_dp) then                ! if there is available capacity in this layer, add deposition here

src/Data/DataOutputModule.f90

@@ -207,7 +207,8 @@ subroutine updateWaterDataOutput(me, t, tInChunk, x, y, date, easts, norths)
                                 trim(str(sum(reach%j_spm%inflow))) // "," // trim(str(sum(reach%j_spm%outflow))) // "," // &
                                 trim(str(sum(reach%j_spm%bankErosion))) // ","
                         end if
-                        write(iouOutputWater, '(a)') trim(str(reach%volume)) // "," // trim(str(reach%depth)) // "," // &
+                        write(iouOutputWater, '(a)') trim(str(reach%volume)) &
+                            // "," // trim(str(divideCheckZero(reach%volume, reach%bedArea))) // "," // &
                             trim(str(reach%Q%outflow / C%timeStep))
                     end associate
                 end do

src/Data/NetCDFOutputModule.f90

@@ -205,7 +205,7 @@ subroutine updateWaterNetCDFOutput(me, t, tInChunk, x, y)
                         me%output_water__m_spm_bank_erosion(w,x,y,tInChunk) = sum(reach%j_spm%bankErosion)
                     end if
                     me%output_water__volume(w,x,y,tInChunk) = reach%volume
-                    me%output_water__depth(w,x,y,tInChunk) = reach%depth
+                    me%output_water__depth(w,x,y,tInChunk) = divideCheckZero(reach%volume, reach%bedArea)
                     me%output_water__flow(w,x,y,tInChunk) = reach%Q%outflow / C%timeStep
                 ! If we're in iterative write mode, then write straight to the NetCDF file, which is time-indexed
                 ! by the whole batch, not just this chunk
@@ -235,7 +235,7 @@ subroutine updateWaterNetCDFOutput(me, t, tInChunk, x, y)
                         call me%nc__water__m_spm_bank_erosion%setData(sum(reach%j_spm%bankErosion), start=[w,x,y,t])
                     end if 
                     call me%nc__water__volume%setData(reach%volume, start=[w,x,y,t])
-                    call me%nc__water__depth%setData(reach%depth, start=[w,x,y,t])
+                    call me%nc__water__depth%setData(divideCheckZero(reach%volume, reach%bedArea), start=[w,x,y,t])
                     call me%nc__water__flow%setData(reach%Q%outflow / C%timeStep, start=[w,x,y,t])
                 end if
             end associate

src/GridCell/GridCellModule.f90

@@ -1006,7 +1006,11 @@ function getWaterDepthGridCell(me) result(waterDepth)
         integer             :: i
         ! Loop over reaches and get their depths and lengths
         do i = 1, me%nReaches
-            waterDepth_i(i) = me%colRiverReaches(i)%item%depth
+            ! Report the hydraulic depth implied by the reach volume, rather than me%depth:
+            ! depositToBed subtracts the water displaced into the bed sediment from me%depth,
+            ! which leaves the reported depth inconsistent with the reported volume and bed area
+            waterDepth_i(i) = divideCheckZero(me%colRiverReaches(i)%item%volume, &
+                                             me%colRiverReaches(i)%item%bedArea)
             lengths(i) = me%colRiverReaches(i)%item%length
         end do
         waterDepth = weightedAverage(waterDepth_i, lengths)

src/VersionModule.f90

@@ -1,4 +1,4 @@
 module VersionModule
     implicit none
-    character(len=17), parameter :: MODEL_VERSION = "0.0.6" 
+    character(len=17), parameter :: MODEL_VERSION = "0.0.7"
 end module

src/WaterBody/ReachModule.f90

@@ -275,17 +275,20 @@ function depositToBedReach(me, spmDep) result(rslt)
             )
         end do
 
+        ! Default to no water being deposited, in case we're not modelling bed sediment
+        V_water_toDeposit = 0.0_dp
         if (C%includeBedSediment) then
             ! Deposit the fine sediment to the bed sediment
             depositRslt = Me%bedSediment%deposit(fineSed)
             call rslt%addErrors(.errors. depositRslt)
             if (rslt%hasCriticalError()) then
                 return
             end if
+            ! Retrieve the amount of water to be taken from the reach. This must be inside the
+            ! includeBedSediment block: depositRslt is only set by the deposit call above, and
+            ! asking an unset Result0D for its data is an error stop
+            V_water_toDeposit = .dp. depositRslt            ! [m3/m2]
         end if
-        ! TODO add error handling to line above as it causes a crash if there is a critical error in the called method
-        ! Retrieve the amount of water to be taken from the reach
-        V_water_toDeposit = .dp. depositRslt                ! [m3/m2]
         ! Subtract that volume for the reach (as a depth). This doesn't have any effect on
         ! the model calculations, as the model recalculates depth depth on hydrology at the
         ! start of every timestep. However, it is this updated depth that is saved to data.