From 14d1dc8e4ca03d5c63d06a0209011b2376692847 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Thu, 5 Dec 2024 11:34:35 +0100
Subject: [PATCH 01/14] Ignore profiles for salinity, simply add by volumes.

---
 src/sealock.c | 19 +++----------------
 1 file changed, 3 insertions(+), 16 deletions(-)

diff --git a/src/sealock.c b/src/sealock.c
index 8def704..33b3b71 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -434,29 +434,16 @@ static double sealock_collect(profile_t *profile, dfm_volumes_t *volumes, double
 
   int first = volumes->first_active_cell;
   int last = first + volumes->num_active_cells - 1;
-  double previous_volume = 0.0;
-  double next_volume = 0.0;
   double used_volume = 0.0;
-  double used_fraction = 0.0;
 
   for (int i = first; i <= last; i++) {
-    next_volume += volumes->normalized[i];
-    // Determine the faction in the direction of the lock according to the profile.
-    // Since this is the negative part of the profile, we need to add a minus sign.
-    double fraction = -1.0 * integrate_piecewise_linear_profile(profile, previous_volume, next_volume);
-    log_debug("Layer %d, value = %g, volume = %g, profile fraction = %g\n", i, buffer_ptr[i], volumes->normalized[i], fraction);
-    // Only aggregate fractions that flow into to the lock.
-    if (fraction > 0) {
-      aggregate += buffer_ptr[i] * volumes->normalized[i] * fraction;
-      used_volume += volumes->normalized[i];
-      used_fraction += fraction;
-    }
-    previous_volume = next_volume;
+    log_debug("Layer %d, value = %g, volume = %g, profile fraction = %g\n", i, buffer_ptr[i], volumes->normalized[i]);
+    aggregate += buffer_ptr[i] * volumes->normalized[i];
+    used_volume += volumes->normalized[i];
   }
 
   log_debug("Aggregate = %g\n", aggregate);
   log_debug("Used volume = %g\n", used_volume);
-  log_debug("Used fraction = %g\n", used_fraction);
 
   // Adjust aggregate to normalized volume.
   // Note: if used_volume is zero, so is aggregate.

From 3268d473c99b867183d2edac9485c9c39adbbdda Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Thu, 5 Dec 2024 16:45:20 +0100
Subject: [PATCH 02/14] Rewrite of layer distribution, taking into account the
 sign of the quantity in relation to the profile.

---
 src/io_layer_distribution.c        | 71 +++++++++++++++++++-----------
 src/io_layer_distribution.h        |  9 +++-
 tests/test_io_layer_distribution.c |  6 +++
 3 files changed, 59 insertions(+), 27 deletions(-)

diff --git a/src/io_layer_distribution.c b/src/io_layer_distribution.c
index 559750f..8fe4ab3 100644
--- a/src/io_layer_distribution.c
+++ b/src/io_layer_distribution.c
@@ -30,7 +30,7 @@ void cleanup_layered_discharge(layered_discharge_t *layered_discharge) {
   free(layered_discharge->discharge_per_layer);
 }
 
-// Generate a default array linear z positions for supplied number_of_layers.
+// Generate a default array linear relative_z positions for supplied number_of_layers.
 // Returns pointer to allocated array of doubles.
 // Note: Caller is responsible to deallocating the created array.
 double *io_layer_linear_z_positions(const int number_of_layers) {
@@ -139,6 +139,9 @@ int io_normalize_profile(profile_t *profile) {
   for (int i = index_after_zero; i <= max_index; i++) {
     profile->relative_discharge_from_lock[i] *= correction_factor_after;
   }
+  // Keep zero information.
+  profile->relative_z_zero = z_zero;
+  profile->start_sign = start_sign;
 
   // Integrate top.
   const double check_integral_before =
@@ -203,16 +206,16 @@ double integrate_piecewise_linear_profile(const profile_t *profile, const double
   return integrated_profile;
 }
 
-// Distribute the total_discharge over layers.
+// Distribute the total_quantity over layers.
 // The discharge that each layer receives is given by a relative profile
-int distribute_discharge_over_layers(double total_discharge, const profile_t *profile,
+int distribute_discharge_over_layers(double total_quantity, const profile_t *profile,
                                      const layers_t *layers,
-                                     layered_discharge_t *layered_discharge_result) {
+                                     layered_discharge_t *layered_quantity_result) {
   assert(layers != NULL);
   assert(layers->number_of_layers > 0);
-  assert(layered_discharge_result != NULL);
-  assert(layered_discharge_result->number_of_layers == layers->number_of_layers);
-  assert(layered_discharge_result->discharge_per_layer != NULL);
+  assert(layered_quantity_result != NULL);
+  assert(layered_quantity_result->number_of_layers == layers->number_of_layers);
+  assert(layered_quantity_result->discharge_per_layer != NULL);
   assert(profile->relative_z_position != NULL);
   assert(profile->relative_discharge_from_lock != NULL);
   assert(profile->number_of_positions > 0);
@@ -222,30 +225,48 @@ int distribute_discharge_over_layers(double total_discharge, const profile_t *pr
     return -1;
   }
 
-  double previous_volume = 0.0;
-  double next_volume = 0.0;
   double profile_integral = 0;
+  double relative_z_start = 0.0;
+  double relative_z_end = 1.0;
+  double relative_z_prev = 0.0;
+  double relative_z = 0.0;
 
-  log_debug("quantity total = %g\n", total_discharge);
+  if (fsign(total_quantity) == profile->start_sign) {
+    relative_z_end = profile->relative_z_zero;
+  } else {
+    relative_z_start = profile->relative_z_zero;
+  }
+
+  double relative_discharge_layer = 0.0;
+
+  log_debug("quantity total = %g\n", total_quantity);
   log_debug("num_layers = %d\n", layers->number_of_layers);
-  for (int layer = 0; layer < layers->number_of_layers; ++layer) {
-    log_debug("normalized_layer_volume[%d] = %g\n", layer,
-              layers->normalized_target_volumes[layer]);
-    next_volume += layers->normalized_target_volumes[layer];
-    const double relative_discharge_layer =
-        integrate_piecewise_linear_profile(profile, previous_volume, next_volume);
-    log_debug("profile_layer_volume   [%d] = %g\n", layer, relative_discharge_layer);
-    double layer_discharge = 0;
-    if (relative_discharge_layer * total_discharge > 0 || layers->number_of_layers == 1) {
-      layer_discharge = fabs(relative_discharge_layer) * total_discharge;
+  if (layers->number_of_layers > 1) {
+    for (int layer = 0; layer < layers->number_of_layers; ++layer) {
+      relative_z += layers->normalized_target_volumes[layer];
+      if (relative_z < relative_z_start || relative_z_prev > relative_z_end) {
+        relative_discharge_layer = 0.0;
+      } else if (relative_z >= relative_z_start && relative_z_prev < relative_z_start) {
+        relative_discharge_layer = integrate_piecewise_linear_profile(profile, relative_z_start, relative_z);
+      } else if (relative_z > relative_z_end && relative_z_prev >= relative_z_start) {
+        relative_discharge_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z_end);
+      } else {
+        relative_discharge_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z);
+      }
+      const double layer_discharge = fabs(relative_discharge_layer) * total_quantity;
       profile_integral += relative_discharge_layer;
+      layered_quantity_result->discharge_per_layer[layer] = layer_discharge;
+      log_debug("layer quantity         [%d] = %g (== %g * %g)\n", layer, layer_discharge,
+                fabs(relative_discharge_layer), total_quantity);
+      relative_z_prev = relative_z;
     }
-    layered_discharge_result->discharge_per_layer[layer] = layer_discharge;
-    log_debug("layer quantity         [%d] = %g (?= %g * %g)\n", layer, layer_discharge,
-              fabs(relative_discharge_layer), total_discharge);
-    previous_volume = next_volume;
+  } else {
+    layered_quantity_result->discharge_per_layer[0] = total_quantity;
+    profile_integral = 1.0;
+    log_debug("layer quantity         [0] = %g (== 1.0 * %g)\n", total_quantity, total_quantity);
   }
+
   log_debug("profile integral = %g (should be +/-1)\n\n", profile_integral);
-  
+
   return 0;
 }
diff --git a/src/io_layer_distribution.h b/src/io_layer_distribution.h
index b6395b5..a3d342d 100644
--- a/src/io_layer_distribution.h
+++ b/src/io_layer_distribution.h
@@ -13,8 +13,13 @@ static const struct profile_struct {
       relative_z_position; // Relative positions in the z-direction for which profile is provided (0 = bed level, 1 = water level)
   double *
       relative_discharge_from_lock; // Relative discharge at corresponding relative z-position (non-negative), need not be normalized
-} profile_default = {
-    .number_of_positions = 0, .relative_z_position = NULL, .relative_discharge_from_lock = NULL};
+  double relative_z_zero;
+  int start_sign;
+} profile_default = {.number_of_positions = 0,
+                     .relative_z_position = NULL,
+                     .relative_discharge_from_lock = NULL,
+                     .relative_z_zero = 1.0,
+                     .start_sign = 0};
 
 typedef struct profile_struct profile_t;
 
diff --git a/tests/test_io_layer_distribution.c b/tests/test_io_layer_distribution.c
index b10ff88..baaf6ad 100644
--- a/tests/test_io_layer_distribution.c
+++ b/tests/test_io_layer_distribution.c
@@ -14,6 +14,12 @@ static int create_linear_profile(profile_t *profile, const int number_of_positio
   profile->number_of_positions = number_of_positions;
   profile->relative_z_position = malloc(number_of_positions * sizeof(double));
   profile->relative_discharge_from_lock = malloc(number_of_positions * sizeof(double));
+  profile->start_sign = discharge_0 < 0 ? -1 : 1;
+  if (discharge_0 * discharge_1 < 0) {
+    profile->relative_z_zero = discharge_0 / (discharge_0 - discharge_1);
+  } else {
+    profile->relative_z_zero = discharge_0 < 0 ? 0.0 : 1.0;
+  }
 
   if (profile->relative_z_position == NULL || profile->relative_discharge_from_lock == NULL) {
     return -1;

From a1f18a5461297f4609ca9ac2a3d8a941355a003e Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Thu, 5 Dec 2024 17:30:01 +0100
Subject: [PATCH 03/14] Removed usage of profiles for concentration quantities
 (salinity).

---
 src/sealock.c | 35 +++++++++++++++++++++++------------
 1 file changed, 23 insertions(+), 12 deletions(-)

diff --git a/src/sealock.c b/src/sealock.c
index 33b3b71..9160f56 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -426,7 +426,7 @@ void sealock_get_active_layers(sealock_state_t *lock) {
 }
 
 // Collect aggregate of quantity from dfm layer buffer into scalar.
-static double sealock_collect(profile_t *profile, dfm_volumes_t *volumes, double* buffer_ptr) { 
+static double sealock_collect(dfm_volumes_t *volumes, double* buffer_ptr) {
   double aggregate = 0.0;
 
   assert(volumes);
@@ -469,7 +469,7 @@ static int sealock_collect_layers(sealock_state_t *lock) {
   } else {
     log_debug("Collecting salinity_lake from layers:\n");
     lock->parameters.salinity_lake =
-        sealock_collect(&lock->flow_profile, lake_volumes, lock->parameters3d.salinity_lake);
+        sealock_collect(lake_volumes, lock->parameters3d.salinity_lake);
   }
 
   sea_volumes = &lock->sea_volumes;
@@ -478,14 +478,16 @@ static int sealock_collect_layers(sealock_state_t *lock) {
   } else {
     log_debug("Collecting salinity_sea from layers:\n");
     lock->parameters.salinity_sea =
-        sealock_collect(&lock->flow_profile, sea_volumes, lock->parameters3d.salinity_sea);
+        sealock_collect(sea_volumes, lock->parameters3d.salinity_sea);
   }
 
   return SEALOCK_OK;
 }
 
-// Distribute scalar quantity over layers into layer buffer.
-static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double quantity, double* buffer_ptr) {
+// Distribute scalar quantity over layers into layer buffer using a provided profile.
+// If no profile is provided (NULL), the quanity is simply duplicated.
+static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double quantity,
+                              double *buffer_ptr) {
   layers_t layers;
   layered_discharge_t result;
   unsigned first_active;
@@ -497,12 +499,21 @@ static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double
   memset(buffer_ptr, 0, MAX_NUM_VOLUMES * sizeof(double));
 
   first_active = volumes->first_active_cell;
-  layers.number_of_layers = volumes->num_active_cells;
-  layers.normalized_target_volumes = volumes->normalized;
-  result.number_of_layers = layers.number_of_layers;
-  result.discharge_per_layer = &buffer_ptr[first_active];
 
-  return distribute_discharge_over_layers(quantity, profile, &layers, &result);
+  if (profile) {
+    layers.number_of_layers = volumes->num_active_cells;
+    layers.normalized_target_volumes = volumes->normalized;
+    result.number_of_layers = layers.number_of_layers;
+    result.discharge_per_layer = &buffer_ptr[first_active];
+
+    return distribute_discharge_over_layers(quantity, profile, &layers, &result);
+  }
+
+  for (int i = 0; i < volumes->num_active_cells; i++) {
+    buffer_ptr[first_active + i] = quantity;
+  }
+
+  return SEALOCK_OK;
 }
 
 static int sealock_distribute_results(sealock_state_t *lock) {
@@ -536,7 +547,7 @@ static int sealock_distribute_results(sealock_state_t *lock) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salinity_to_lake\n");
-  if (sealock_distribute(lake_volumes, lake_profile, lock->results.salinity_to_lake,
+  if (sealock_distribute(lake_volumes, NULL, lock->results.salinity_to_lake,
                          lock->results3d.salinity_to_lake) != 0) {
     return SEALOCK_ERROR;
   }
@@ -561,7 +572,7 @@ static int sealock_distribute_results(sealock_state_t *lock) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salinity_to_sea\n");
-  if (sealock_distribute(sea_volumes, sea_profile, lock->results.salinity_to_sea,
+  if (sealock_distribute(sea_volumes, NULL, lock->results.salinity_to_sea,
                          lock->results3d.salinity_to_sea) != 0) {
     return SEALOCK_ERROR;
   }

From 189dc3bd15dc33a31ef3fa4ebf6566b81cf90a4a Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Thu, 5 Dec 2024 17:46:44 +0100
Subject: [PATCH 04/14] Improved logging of salinity layers.

---
 src/sealock.c | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/sealock.c b/src/sealock.c
index 9160f56..795e220 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -437,7 +437,7 @@ static double sealock_collect(dfm_volumes_t *volumes, double* buffer_ptr) {
   double used_volume = 0.0;
 
   for (int i = first; i <= last; i++) {
-    log_debug("Layer %d, value = %g, volume = %g, profile fraction = %g\n", i, buffer_ptr[i], volumes->normalized[i]);
+    log_debug("Layer %d, value = %g, volume = %g\n", i, buffer_ptr[i], volumes->normalized[i]);
     aggregate += buffer_ptr[i] * volumes->normalized[i];
     used_volume += volumes->normalized[i];
   }
@@ -509,9 +509,11 @@ static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double
     return distribute_discharge_over_layers(quantity, profile, &layers, &result);
   }
 
-  for (int i = 0; i < volumes->num_active_cells; i++) {
-    buffer_ptr[first_active + i] = quantity;
+  for (int layer = 0; layer < volumes->num_active_cells; layer++) {
+    buffer_ptr[first_active + layer] = quantity;
+    log_debug("layer quantity         [%d] = %g\n", layer, quantity);
   }
+  log_debug("(no profile was used)\n\n");
 
   return SEALOCK_OK;
 }

From 711a83eb697a1c0235e07337c830dabeaf489093 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 10 Dec 2024 11:23:13 +0100
Subject: [PATCH 05/14] Split volumes in 'to' and 'from' parts, we use proper
 volumes for the configured laterals.

---
 src/dimr_bmi.c | 74 +++++++++++++++++++++++++++-----------------------
 src/sealock.c  | 73 +++++++++++++++++++++++++++++++------------------
 src/sealock.h  |  6 ++--
 3 files changed, 90 insertions(+), 63 deletions(-)

diff --git a/src/dimr_bmi.c b/src/dimr_bmi.c
index 8e0e676..f8aa45f 100644
--- a/src/dimr_bmi.c
+++ b/src/dimr_bmi.c
@@ -120,20 +120,26 @@ int set_var(const char *key, void *src_ptr) {
   // Set dest_ptr for intended lock and quantity.
   if (match_key(quantity, "salinity_lake")) {
     dest_ptr = config.locks[lock_index].parameters3d.salinity_lake;
-    dest_len = config.locks[lock_index].lake_volumes.num_volumes;
+    dest_len = config.locks[lock_index].from_lake_volumes.num_volumes;
   } else if (match_key(quantity, "head_lake")) {
     dest_ptr = &config.locks[lock_index].parameters.head_lake;
   } else if (match_key(quantity, "salinity_sea")) {
     dest_ptr = config.locks[lock_index].parameters3d.salinity_sea;
-    dest_len = config.locks[lock_index].sea_volumes.num_volumes;
+    dest_len = config.locks[lock_index].from_sea_volumes.num_volumes;
   } else if (match_key(quantity, "head_sea")) {
     dest_ptr = &config.locks[lock_index].parameters.head_sea;
-  } else if (match_key(quantity, "water_volume_lake")) {
-    dest_ptr = config.locks[lock_index].lake_volumes.volumes;
-    dest_len = config.locks[lock_index].lake_volumes.num_volumes;
-  } else if (match_key(quantity, "water_volume_sea")) {
-    dest_ptr = config.locks[lock_index].sea_volumes.volumes;
-    dest_len = config.locks[lock_index].sea_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_from_lake")) {
+    dest_ptr = config.locks[lock_index].from_lake_volumes.volumes;
+    dest_len = config.locks[lock_index].from_lake_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_from_sea")) {
+    dest_ptr = config.locks[lock_index].from_sea_volumes.volumes;
+    dest_len = config.locks[lock_index].from_sea_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_to_lake")) {
+    dest_ptr = config.locks[lock_index].to_lake_volumes.volumes;
+    dest_len = config.locks[lock_index].to_lake_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_to_sea")) {
+    dest_ptr = config.locks[lock_index].to_sea_volumes.volumes;
+    dest_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "temperature_lake")) {
     dest_ptr = &config.locks[lock_index].parameters.temperature_lake;
   } else if (match_key(quantity, "temperature_sea")) {
@@ -185,40 +191,40 @@ int get_var(const char *key, void **dst_ptr) {
   // Set source based on key(s)...
   if (match_key(quantity, "mass_transport_lake")) {
     source_ptr = config.locks[lock_index].results3d.mass_transport_lake;
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "salt_load_lake")) {
     source_ptr = config.locks[lock_index].results3d.salt_load_lake;
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_from_lake")) {
     source_ptr = config.locks[lock_index].results3d.discharge_from_lake;
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_to_lake")) {
     source_ptr = config.locks[lock_index].results3d.discharge_to_lake;
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_to_lake")) {
     source_ptr = config.locks[lock_index].results3d.salinity_to_lake;
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "mass_transport_sea")) {
     source_ptr = config.locks[lock_index].results3d.mass_transport_sea;
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salt_load_sea")) {
     source_ptr = config.locks[lock_index].results3d.salt_load_sea;
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_from_sea")) {
     source_ptr = config.locks[lock_index].results3d.discharge_from_sea;
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_to_sea")) {
     source_ptr = config.locks[lock_index].results3d.discharge_to_sea;
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_to_sea")) {
     source_ptr = config.locks[lock_index].results3d.salinity_to_sea;
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "water_volume_lake")) {
-    source_ptr = config.locks[lock_index].lake_volumes.volumes;
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_ptr = config.locks[lock_index].to_lake_volumes.volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "water_volume_sea")) {
-    source_ptr = config.locks[lock_index].sea_volumes.volumes;
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_ptr = config.locks[lock_index].to_sea_volumes.volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_sea")) {
     // NOTE: This is really a GET_VALUE_PTR(), called before the update.
     source_ptr = config.locks[lock_index].parameters3d.salinity_sea;
@@ -308,29 +314,29 @@ int get_var_shape(char *key, int dims[DIMR_BMI_MAXDIMS]) { // dims -> int[6]
 
   // Set source based on key(s)...
   if (match_key(quantity, "mass_transport_lake")) {
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "salt_load_lake")) {
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_from_lake")) {
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_to_lake")) {
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_to_lake")) {
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
   } else if (match_key(quantity, "mass_transport_sea")) {
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salt_load_sea")) {
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_from_sea")) {
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
   } else if (match_key(quantity, "discharge_to_sea")) {
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_to_sea")) {
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "water_volume_lake")) {
-    source_len = config.locks[lock_index].lake_volumes.num_volumes;
+    source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
   } else if (match_key(quantity, "water_volume_sea")) {
-    source_len = config.locks[lock_index].sea_volumes.num_volumes;
+    source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
   } else {
     log_debug("Unhandled get_var('%s', @%p)\n", key, dims);
     source_len = 1;
diff --git a/src/sealock.c b/src/sealock.c
index 795e220..696edf5 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -15,14 +15,22 @@ int sealock_defaults(sealock_state_t* lock) {
   // Init calculation parameters with defaults.
   zsf_param_default(&lock->parameters);
   // Set up default volumes/profile for '2D' case.
-  lock->lake_volumes.num_volumes = 1;
-  lock->lake_volumes.volumes[0] = 1.0;
-  lock->lake_volumes.first_active_cell = 0;
-  lock->lake_volumes.num_active_cells = 1;
-  lock->sea_volumes.num_volumes = 1;
-  lock->sea_volumes.volumes[0] = 1.0;
-  lock->sea_volumes.first_active_cell = 0;
-  lock->sea_volumes.num_active_cells = 1;
+  lock->from_lake_volumes.num_volumes = 1;
+  lock->from_lake_volumes.volumes[0] = 1.0;
+  lock->from_lake_volumes.first_active_cell = 0;
+  lock->from_lake_volumes.num_active_cells = 1;
+  lock->from_sea_volumes.num_volumes = 1;
+  lock->from_sea_volumes.volumes[0] = 1.0;
+  lock->from_sea_volumes.first_active_cell = 0;
+  lock->from_sea_volumes.num_active_cells = 1;
+  lock->to_lake_volumes.num_volumes = 1;
+  lock->to_lake_volumes.volumes[0] = 1.0;
+  lock->to_lake_volumes.first_active_cell = 0;
+  lock->to_lake_volumes.num_active_cells = 1;
+  lock->to_sea_volumes.num_volumes = 1;
+  lock->to_sea_volumes.volumes[0] = 1.0;
+  lock->to_sea_volumes.first_active_cell = 0;
+  lock->to_sea_volumes.num_active_cells = 1;
   return io_layer_init_2d(&lock->flow_profile);
 }
 
@@ -31,8 +39,10 @@ int sealock_init(sealock_state_t* lock, time_t start_time, unsigned int max_num_
   int status = SEALOCK_OK;
 
   // Set number of DFM volumes.
-  lock->lake_volumes.num_volumes = max_num_z_layers;
-  lock->sea_volumes.num_volumes = max_num_z_layers;
+  lock->from_lake_volumes.num_volumes = max_num_z_layers;
+  lock->from_sea_volumes.num_volumes = max_num_z_layers;
+  lock->to_lake_volumes.num_volumes = max_num_z_layers;
+  lock->to_sea_volumes.num_volumes = max_num_z_layers;
 
   // Load timeseries data when required.
   if (status == SEALOCK_OK && lock->operational_parameters_file) {
@@ -421,8 +431,14 @@ static void sealock_get_active_cells(dfm_volumes_t* volumes) {
 // Underlying assumption is that the volumes have one
 // contiguous run of non-zero cells.
 void sealock_get_active_layers(sealock_state_t *lock) {
-  sealock_get_active_cells(&lock->lake_volumes);
-  sealock_get_active_cells(&lock->sea_volumes);
+  log_debug("from_lake_volumes:\n");
+  sealock_get_active_cells(&lock->from_lake_volumes);
+  log_debug("from_sea _volumes:\n");
+  sealock_get_active_cells(&lock->from_sea_volumes);
+  log_debug("to_lake_volumes:\n");
+  sealock_get_active_cells(&lock->to_lake_volumes);
+  log_debug("to_sea_volumes:\n");
+  sealock_get_active_cells(&lock->to_sea_volumes);
 }
 
 // Collect aggregate of quantity from dfm layer buffer into scalar.
@@ -463,7 +479,7 @@ static int sealock_collect_layers(sealock_state_t *lock) {
   assert(lock);
   sealock_get_active_layers(lock);
 
-  lake_volumes = &lock->lake_volumes;
+  lake_volumes = &lock->from_lake_volumes;
   if (lake_volumes->num_volumes == 1) {
     lock->parameters.salinity_lake = lock->parameters3d.salinity_lake[lake_volumes->first_active_cell]; 
   } else {
@@ -472,7 +488,7 @@ static int sealock_collect_layers(sealock_state_t *lock) {
         sealock_collect(lake_volumes, lock->parameters3d.salinity_lake);
   }
 
-  sea_volumes = &lock->sea_volumes;
+  sea_volumes = &lock->from_sea_volumes;
   if (sea_volumes->num_volumes == 1) {
     lock->parameters.salinity_sea = lock->parameters3d.salinity_sea[sea_volumes->first_active_cell];
   } else {
@@ -520,61 +536,64 @@ static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double
 
 static int sealock_distribute_results(sealock_state_t *lock) {
   profile_t *lake_profile, *sea_profile;
-  dfm_volumes_t *lake_volumes, *sea_volumes;
+  dfm_volumes_t *from_lake_volumes, *from_sea_volumes;
+  dfm_volumes_t *to_lake_volumes, *to_sea_volumes;
 
   assert(lock);
 
   lake_profile = &lock->flow_profile;
   sea_profile = &lock->flow_profile;
-  lake_volumes = &lock->lake_volumes;
-  sea_volumes = &lock->sea_volumes;
+  from_lake_volumes = &lock->from_lake_volumes;
+  from_sea_volumes = &lock->from_sea_volumes;
+  to_lake_volumes = &lock->to_lake_volumes;
+  to_sea_volumes = &lock->to_sea_volumes;
 
   log_debug("quantity = mass_transport_lake\n");
-  if (sealock_distribute(lake_volumes, lake_profile, lock->results.mass_transport_lake, lock->results3d.mass_transport_lake) != 0) {
+  if (sealock_distribute(to_lake_volumes, lake_profile, lock->results.mass_transport_lake, lock->results3d.mass_transport_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salt_load_lake\n");
-  if (sealock_distribute(lake_volumes, lake_profile, lock->results.salt_load_lake,
+  if (sealock_distribute(to_lake_volumes, lake_profile, lock->results.salt_load_lake,
                          lock->results3d.salt_load_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_from_lake\n");
-  if (sealock_distribute(lake_volumes, lake_profile, lock->results.discharge_from_lake,
+  if (sealock_distribute(from_lake_volumes, lake_profile, lock->results.discharge_from_lake,
                              lock->results3d.discharge_from_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_to_lake\n");
-  if (sealock_distribute(lake_volumes, lake_profile, lock->results.discharge_to_lake,
+  if (sealock_distribute(from_lake_volumes, lake_profile, lock->results.discharge_to_lake,
                          lock->results3d.discharge_to_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salinity_to_lake\n");
-  if (sealock_distribute(lake_volumes, NULL, lock->results.salinity_to_lake,
+  if (sealock_distribute(to_lake_volumes, NULL, lock->results.salinity_to_lake,
                          lock->results3d.salinity_to_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = mass_transport_sea\n");
-  if (sealock_distribute(sea_volumes, sea_profile, lock->results.mass_transport_sea,
+  if (sealock_distribute(to_sea_volumes, sea_profile, lock->results.mass_transport_sea,
                          lock->results3d.mass_transport_sea) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salt_load_sea\n");
-  if (sealock_distribute(sea_volumes, sea_profile, lock->results.salt_load_sea,
+  if (sealock_distribute(to_sea_volumes, sea_profile, lock->results.salt_load_sea,
                          lock->results3d.salt_load_sea) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_from_sea\n");
-  if (sealock_distribute(sea_volumes, sea_profile, lock->results.discharge_from_sea,
+  if (sealock_distribute(from_sea_volumes, sea_profile, lock->results.discharge_from_sea,
                          lock->results3d.discharge_from_sea) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_to_sea\n");
-  if (sealock_distribute(sea_volumes, sea_profile, lock->results.discharge_to_sea,
+  if (sealock_distribute(to_sea_volumes, sea_profile, lock->results.discharge_to_sea,
                          lock->results3d.discharge_to_sea) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salinity_to_sea\n");
-  if (sealock_distribute(sea_volumes, NULL, lock->results.salinity_to_sea,
+  if (sealock_distribute(to_sea_volumes, NULL, lock->results.salinity_to_sea,
                          lock->results3d.salinity_to_sea) != 0) {
     return SEALOCK_ERROR;
   }
diff --git a/src/sealock.h b/src/sealock.h
index ecf0984..1a73b0a 100644
--- a/src/sealock.h
+++ b/src/sealock.h
@@ -89,8 +89,10 @@ typedef struct sealock_state_struct {
   time_t *times;
   size_t times_len;
   profile_t flow_profile;
-  dfm_volumes_t lake_volumes;
-  dfm_volumes_t sea_volumes;
+  dfm_volumes_t from_lake_volumes;
+  dfm_volumes_t from_sea_volumes;
+  dfm_volumes_t to_lake_volumes;
+  dfm_volumes_t to_sea_volumes;
   dfm_parameters3d_t parameters3d;
   dfm_results3d_t results3d;
 } sealock_state_t;

From e41cb845ed7e545d0ad935ba085a9f69a69aff26 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 10 Dec 2024 11:38:14 +0100
Subject: [PATCH 06/14] Added masking for salinity layers based on discharge
 result.

---
 src/sealock.c | 52 ++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 41 insertions(+), 11 deletions(-)

diff --git a/src/sealock.c b/src/sealock.c
index 696edf5..79359c4 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -10,6 +10,7 @@
 #include <string.h>
 #include <stdio.h>
 #include <float.h>
+#include <math.h>
 
 int sealock_defaults(sealock_state_t* lock) {
   // Init calculation parameters with defaults.
@@ -494,7 +495,7 @@ static int sealock_collect_layers(sealock_state_t *lock) {
   } else {
     log_debug("Collecting salinity_sea from layers:\n");
     lock->parameters.salinity_sea =
-        sealock_collect(sea_volumes, lock->parameters3d.salinity_sea);
+sealock_collect(sea_volumes, lock->parameters3d.salinity_sea);
   }
 
   return SEALOCK_OK;
@@ -502,8 +503,8 @@ static int sealock_collect_layers(sealock_state_t *lock) {
 
 // Distribute scalar quantity over layers into layer buffer using a provided profile.
 // If no profile is provided (NULL), the quanity is simply duplicated.
-static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double quantity,
-                              double *buffer_ptr) {
+static int sealock_distribute(dfm_volumes_t* volumes, profile_t* profile, double quantity,
+  double* buffer_ptr) {
   layers_t layers;
   layered_discharge_t result;
   unsigned first_active;
@@ -534,10 +535,31 @@ static int sealock_distribute(dfm_volumes_t *volumes, profile_t *profile, double
   return SEALOCK_OK;
 }
 
-static int sealock_distribute_results(sealock_state_t *lock) {
-  profile_t *lake_profile, *sea_profile;
-  dfm_volumes_t *from_lake_volumes, *from_sea_volumes;
-  dfm_volumes_t *to_lake_volumes, *to_sea_volumes;
+// Set layer entries in buffer to zero if mask buffer is zero.
+static int sealock_mask_layers(dfm_volumes_t* volumes, double* mask_ptr, double* buffer_ptr) {
+  layers_t layers;
+  layered_discharge_t result;
+  unsigned first_active;
+
+  assert(volumes);
+  assert(buffer_ptr);
+
+  first_active = volumes->first_active_cell;
+
+  for (int layer = 0; layer < volumes->num_active_cells; layer++) {
+    if (fabs(mask_ptr[layer]) < DBL_EPSILON) {
+      buffer_ptr[first_active + layer] = 0.0;
+    }
+    log_debug("masked layer quantity         [%d] = %g\n", layer, buffer_ptr[first_active + layer]);
+  }
+
+  return SEALOCK_OK;
+}
+
+static int sealock_distribute_results(sealock_state_t* lock) {
+  profile_t* lake_profile, * sea_profile;
+  dfm_volumes_t* from_lake_volumes, * from_sea_volumes;
+  dfm_volumes_t* to_lake_volumes, * to_sea_volumes;
 
   assert(lock);
 
@@ -554,22 +576,26 @@ static int sealock_distribute_results(sealock_state_t *lock) {
   }
   log_debug("quantity = salt_load_lake\n");
   if (sealock_distribute(to_lake_volumes, lake_profile, lock->results.salt_load_lake,
-                         lock->results3d.salt_load_lake) != 0) {
+    lock->results3d.salt_load_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_from_lake\n");
   if (sealock_distribute(from_lake_volumes, lake_profile, lock->results.discharge_from_lake,
-                             lock->results3d.discharge_from_lake) != 0) {
+    lock->results3d.discharge_from_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_to_lake\n");
   if (sealock_distribute(from_lake_volumes, lake_profile, lock->results.discharge_to_lake,
-                         lock->results3d.discharge_to_lake) != 0) {
+    lock->results3d.discharge_to_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = salinity_to_lake\n");
   if (sealock_distribute(to_lake_volumes, NULL, lock->results.salinity_to_lake,
-                         lock->results3d.salinity_to_lake) != 0) {
+    lock->results3d.salinity_to_lake) != 0) {
+    return SEALOCK_ERROR;
+  }
+  if (sealock_mask_layers(to_lake_volumes, lock->results3d.discharge_to_lake,
+    lock->results3d.salinity_to_lake) != 0) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = mass_transport_sea\n");
@@ -597,6 +623,10 @@ static int sealock_distribute_results(sealock_state_t *lock) {
                          lock->results3d.salinity_to_sea) != 0) {
     return SEALOCK_ERROR;
   }
+  if (sealock_mask_layers(to_sea_volumes, lock->results3d.discharge_to_sea,
+                          lock->results3d.salinity_to_sea) != 0) {
+    return SEALOCK_ERROR;
+  }
   return SEALOCK_OK;
 }
 

From c7c7f360f2588d960836dad02694eb0b5c811f23 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 10 Dec 2024 13:21:56 +0100
Subject: [PATCH 07/14] Fixed typo from->to.

---
 src/sealock.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/sealock.c b/src/sealock.c
index 79359c4..76fa3ed 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -585,7 +585,7 @@ static int sealock_distribute_results(sealock_state_t* lock) {
     return SEALOCK_ERROR;
   }
   log_debug("quantity = discharge_to_lake\n");
-  if (sealock_distribute(from_lake_volumes, lake_profile, lock->results.discharge_to_lake,
+  if (sealock_distribute(to_lake_volumes, lake_profile, lock->results.discharge_to_lake,
     lock->results3d.discharge_to_lake) != 0) {
     return SEALOCK_ERROR;
   }

From 826bac0699f9d6cee89f74a0b34081e2ee99f9e7 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 10 Dec 2024 17:49:03 +0100
Subject: [PATCH 08/14] Fixed quantity names for get

---
 src/dimr_bmi.c | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/dimr_bmi.c b/src/dimr_bmi.c
index f8aa45f..77d8709 100644
--- a/src/dimr_bmi.c
+++ b/src/dimr_bmi.c
@@ -219,10 +219,16 @@ int get_var(const char *key, void **dst_ptr) {
   } else if (match_key(quantity, "salinity_to_sea")) {
     source_ptr = config.locks[lock_index].results3d.salinity_to_sea;
     source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
-  } else if (match_key(quantity, "water_volume_lake")) {
+  } else if (match_key(quantity, "water_volume_from_lake")) {
+    source_ptr = config.locks[lock_index].from_lake_volumes.volumes;
+    source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_from_sea")) {
+    source_ptr = config.locks[lock_index].from_sea_volumes.volumes;
+    source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_to_lake")) {
     source_ptr = config.locks[lock_index].to_lake_volumes.volumes;
     source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
-  } else if (match_key(quantity, "water_volume_sea")) {
+  } else if (match_key(quantity, "water_volume_to_sea")) {
     source_ptr = config.locks[lock_index].to_sea_volumes.volumes;
     source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_sea")) {
@@ -333,15 +339,20 @@ int get_var_shape(char *key, int dims[DIMR_BMI_MAXDIMS]) { // dims -> int[6]
     source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else if (match_key(quantity, "salinity_to_sea")) {
     source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
-  } else if (match_key(quantity, "water_volume_lake")) {
+  } else if (match_key(quantity, "water_volume_from_lake")) {
     source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
-  } else if (match_key(quantity, "water_volume_sea")) {
+  } else if (match_key(quantity, "water_volume_from_sea")) {
     source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_to_lake")) {
+    source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
+  } else if (match_key(quantity, "water_volume_to_sea")) {
+    source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
   } else {
     log_debug("Unhandled get_var('%s', @%p)\n", key, dims);
     source_len = 1;
   }
 
+  memset(dims, 0, DIMR_BMI_MAXDIMS * sizeof(int));
   dims[0] = source_len;
   log_info("%s yielded %d for quantity '%s' of lock %d.\n", __func__, source_len,
            quantity, lock_index);

From 9a3e5c35784f7f28a1eca5df0d9ebd27bfb35ab2 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 17 Dec 2024 13:43:59 +0100
Subject: [PATCH 09/14] Renamed variables to adhere to new terminology.

---
 src/io_layer_distribution.c | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/io_layer_distribution.c b/src/io_layer_distribution.c
index 8fe4ab3..d3520c7 100644
--- a/src/io_layer_distribution.c
+++ b/src/io_layer_distribution.c
@@ -30,7 +30,7 @@ void cleanup_layered_discharge(layered_discharge_t *layered_discharge) {
   free(layered_discharge->discharge_per_layer);
 }
 
-// Generate a default array linear relative_z positions for supplied number_of_layers.
+// Generate a default array linear z positions for supplied number_of_layers.
 // Returns pointer to allocated array of doubles.
 // Note: Caller is responsible to deallocating the created array.
 double *io_layer_linear_z_positions(const int number_of_layers) {
@@ -237,7 +237,7 @@ int distribute_discharge_over_layers(double total_quantity, const profile_t *pro
     relative_z_start = profile->relative_z_zero;
   }
 
-  double relative_discharge_layer = 0.0;
+  double relative_quantity_layer = 0.0;
 
   log_debug("quantity total = %g\n", total_quantity);
   log_debug("num_layers = %d\n", layers->number_of_layers);
@@ -245,19 +245,19 @@ int distribute_discharge_over_layers(double total_quantity, const profile_t *pro
     for (int layer = 0; layer < layers->number_of_layers; ++layer) {
       relative_z += layers->normalized_target_volumes[layer];
       if (relative_z < relative_z_start || relative_z_prev > relative_z_end) {
-        relative_discharge_layer = 0.0;
+        relative_quantity_layer = 0.0;
       } else if (relative_z >= relative_z_start && relative_z_prev < relative_z_start) {
-        relative_discharge_layer = integrate_piecewise_linear_profile(profile, relative_z_start, relative_z);
+        relative_quantity_layer = integrate_piecewise_linear_profile(profile, relative_z_start, relative_z);
       } else if (relative_z > relative_z_end && relative_z_prev >= relative_z_start) {
-        relative_discharge_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z_end);
+        relative_quantity_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z_end);
       } else {
-        relative_discharge_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z);
+        relative_quantity_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z);
       }
-      const double layer_discharge = fabs(relative_discharge_layer) * total_quantity;
-      profile_integral += relative_discharge_layer;
-      layered_quantity_result->discharge_per_layer[layer] = layer_discharge;
-      log_debug("layer quantity         [%d] = %g (== %g * %g)\n", layer, layer_discharge,
-                fabs(relative_discharge_layer), total_quantity);
+      const double layer_quantity = fabs(relative_quantity_layer) * total_quantity;
+      profile_integral += relative_quantity_layer;
+      layered_quantity_result->discharge_per_layer[layer] = layer_quantity;
+      log_debug("layer quantity         [%d] = %g (== %g * %g)\n", layer, layer_quantity,
+                fabs(relative_quantity_layer), total_quantity);
       relative_z_prev = relative_z;
     }
   } else {

From c4764520f5daae4ffa400d65cc5c84fceb6077be Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Wed, 18 Dec 2024 13:38:33 +0100
Subject: [PATCH 10/14] Fixed log message for phase wise update step.

---
 src/sealock.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/sealock.c b/src/sealock.c
index 76fa3ed..58188a4 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -324,9 +324,9 @@ static int sealock_update_phase_wise_parameters(sealock_state_t *lock, time_t ti
 static int sealock_phase_wise_step(sealock_state_t *lock, time_t time) {
   int status = SEALOCK_OK;
   time_t duration = 0;
-  log_debug("%s: Handling '%d' (run_update = %d)\n", lock->id, lock->phase_args.run_update);
+  log_debug("%s: Handling '%s' (run_update = %d)\n", __func__, lock->id, lock->phase_args.run_update);
   if (lock->phase_args.run_update) {
-    log_info("%s: Updating '%d' to phase %d.\n", lock->id, lock->phase_args.routine);
+    log_info("%s: Updating '%s' to phase %d.\n", __func__, lock->id, lock->phase_args.routine);
     zsf_results_t previous_step_results = lock->results;
     switch (lock->phase_args.routine) {
     case 1:

From b603f91da35d95f9298d31f00a3f09a0d37628f0 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Thu, 19 Dec 2024 16:20:33 +0100
Subject: [PATCH 11/14] Made open phasse a a bit more lenient towards inner and
 outer head differences (max 1cm), also made the delta_time check towards time
 steps and phase-wise phase lengths a bit more strict to prevent having
 corrected phase lenghts of 0, wich doesn't allow for any correction in a new
 phase.

---
 src/dimr_bmi.c | 1 +
 src/sealock.c  | 2 +-
 src/zsf.c      | 4 ++--
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/dimr_bmi.c b/src/dimr_bmi.c
index 77d8709..c164dce 100644
--- a/src/dimr_bmi.c
+++ b/src/dimr_bmi.c
@@ -279,6 +279,7 @@ int update(double dt) {
   if (config.current_time == config.start_time && delta_time > 0) {
     // Check if timestep is compatible with all loaded timeseries.
     if (!sealock_delta_time_ok(&config.locks[lock_index], delta_time)) {
+      log_error("delta_time (%d) >= shortest phase in the timeseries!\n", delta_time);
       return DIMR_BMI_FAILURE;
     }
   }
diff --git a/src/sealock.c b/src/sealock.c
index 58188a4..0a46674 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -658,7 +658,7 @@ int sealock_update(sealock_state_t *lock, time_t time) {
 // Check if none of the time steps in the timeseries is shorter than delta_time
 int sealock_delta_time_ok(sealock_state_t* lock, time_t delta_time) {
   for (int i = 0; i < lock->times_len-1; i++) {
-    if (lock->times[i + 1] - lock->times[i] < delta_time) {
+    if (lock->times[i + 1] - lock->times[i] <= delta_time) {
       return 0;
     }
   }
diff --git a/src/zsf.c b/src/zsf.c
index be91cb3..4423c59 100644
--- a/src/zsf.c
+++ b/src/zsf.c
@@ -736,7 +736,7 @@ int ZSF_CALLCONV zsf_step_phase_2(const zsf_param_t *p, double t_open_lake,
   if (err) {
     return err;
   }
-  if (fabs(state->head_lock - p->head_lake) > 1E-8) {
+  if (fabs(state->head_lock - p->head_lake) > 0.01 ) { // was: 1E-8
     return ZSF_ERR_REMAINING_HEAD_DIFF;
   }
 
@@ -788,7 +788,7 @@ int ZSF_CALLCONV zsf_step_phase_4(const zsf_param_t *p, double t_open_sea, zsf_p
   if (err) {
     return err;
   }
-  if (fabs(state->head_lock - p->head_sea) > 1E-8) {
+  if (fabs(state->head_lock - p->head_sea) > 0.01) { // Was: 1E-8
     return ZSF_ERR_REMAINING_HEAD_DIFF;
   }
 

From 65aba4705c3e2330388a13c81faef53b378e83f1 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 7 Jan 2025 17:34:55 +0100
Subject: [PATCH 12/14] Fixed division by zero while correcting phase values in
 case lake or sea volumes are zero in phases where that door is closed.

---
 src/sealock.c | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/sealock.c b/src/sealock.c
index 0a46674..1434337 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -235,12 +235,16 @@ static int sealock_apply_phase_wise_result_correction(sealock_state_t *lock, tim
     log_info("Applying correction to discharge_to_sea  : %g -> %g\n",
              lock->results.discharge_to_sea, total_volume_to_sea / new_phase_len);
     lock->results.discharge_to_sea = total_volume_to_sea / new_phase_len;
-    log_info("Applying correction to salinity_to_lake  : %g -> %g\n",
-             lock->results.salinity_to_lake, total_salt_to_lake / total_volume_to_lake);
-    lock->results.salinity_to_lake = total_salt_to_lake / total_volume_to_lake;
-    log_info("Applying correction to salinity_to_sea  : %g -> %g\n",
-             lock->results.salinity_to_sea, total_salt_to_sea / total_volume_to_sea);
-    lock->results.salinity_to_sea = total_salt_to_sea / total_volume_to_sea;
+    if (total_volume_to_lake > 0) {
+      log_info("Applying correction to salinity_to_lake  : %g -> %g\n",
+               lock->results.salinity_to_lake, total_salt_to_lake / total_volume_to_lake);
+      lock->results.salinity_to_lake = total_salt_to_lake / total_volume_to_lake;
+    }
+    if (total_volume_to_sea > 0) {
+      log_info("Applying correction to salinity_to_sea  : %g -> %g\n",
+               lock->results.salinity_to_sea, total_salt_to_sea / total_volume_to_sea);
+      lock->results.salinity_to_sea = total_salt_to_sea / total_volume_to_sea;
+    }
     log_info("Applying correction to discharge_from_lake : %g -> %g\n",
              lock->results.discharge_from_lake, total_volume_from_lake / new_phase_len);
     lock->results.discharge_from_lake = total_volume_from_lake / new_phase_len;

From 5c4e34f63be60ef3e9a6bb2bc3caef06cf50f773 Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 21 Jan 2025 16:07:15 +0100
Subject: [PATCH 13/14] Addressed PR 58 review comments.

---
 src/io_layer_distribution.c |  2 +-
 src/sealock.c               | 27 +++++++++++----------------
 src/sealock.h               |  6 ++++++
 3 files changed, 18 insertions(+), 17 deletions(-)

diff --git a/src/io_layer_distribution.c b/src/io_layer_distribution.c
index d3520c7..4c74e91 100644
--- a/src/io_layer_distribution.c
+++ b/src/io_layer_distribution.c
@@ -207,7 +207,7 @@ double integrate_piecewise_linear_profile(const profile_t *profile, const double
 }
 
 // Distribute the total_quantity over layers.
-// The discharge that each layer receives is given by a relative profile
+// The quantity that each layer receives is given by a relative profile
 int distribute_discharge_over_layers(double total_quantity, const profile_t *profile,
                                      const layers_t *layers,
                                      layered_discharge_t *layered_quantity_result) {
diff --git a/src/sealock.c b/src/sealock.c
index 1434337..41d14b0 100644
--- a/src/sealock.c
+++ b/src/sealock.c
@@ -12,26 +12,21 @@
 #include <float.h>
 #include <math.h>
 
+static void sealock_set_2d_defaults(dfm_volumes_t *volumes) {
+  volumes->num_volumes = 1;
+  volumes->volumes[0] = 1.0;
+  volumes->first_active_cell = 0;
+  volumes->num_active_cells = 1;
+}
+
 int sealock_defaults(sealock_state_t* lock) {
   // Init calculation parameters with defaults.
   zsf_param_default(&lock->parameters);
   // Set up default volumes/profile for '2D' case.
-  lock->from_lake_volumes.num_volumes = 1;
-  lock->from_lake_volumes.volumes[0] = 1.0;
-  lock->from_lake_volumes.first_active_cell = 0;
-  lock->from_lake_volumes.num_active_cells = 1;
-  lock->from_sea_volumes.num_volumes = 1;
-  lock->from_sea_volumes.volumes[0] = 1.0;
-  lock->from_sea_volumes.first_active_cell = 0;
-  lock->from_sea_volumes.num_active_cells = 1;
-  lock->to_lake_volumes.num_volumes = 1;
-  lock->to_lake_volumes.volumes[0] = 1.0;
-  lock->to_lake_volumes.first_active_cell = 0;
-  lock->to_lake_volumes.num_active_cells = 1;
-  lock->to_sea_volumes.num_volumes = 1;
-  lock->to_sea_volumes.volumes[0] = 1.0;
-  lock->to_sea_volumes.first_active_cell = 0;
-  lock->to_sea_volumes.num_active_cells = 1;
+  sealock_set_2d_defaults(&lock->from_lake_volumes);
+  sealock_set_2d_defaults(&lock->from_sea_volumes);
+  sealock_set_2d_defaults(&lock->to_lake_volumes);
+  sealock_set_2d_defaults(&lock->to_sea_volumes);
   return io_layer_init_2d(&lock->flow_profile);
 }
 
diff --git a/src/sealock.h b/src/sealock.h
index 1a73b0a..0d57366 100644
--- a/src/sealock.h
+++ b/src/sealock.h
@@ -89,11 +89,17 @@ typedef struct sealock_state_struct {
   time_t *times;
   size_t times_len;
   profile_t flow_profile;
+  // Volumes
+  // Note that the 'to' and 'from' locations may differ
+  // so their actual water volumes may be different and
+  // should be trated as such.
   dfm_volumes_t from_lake_volumes;
   dfm_volumes_t from_sea_volumes;
   dfm_volumes_t to_lake_volumes;
   dfm_volumes_t to_sea_volumes;
+  // 3D parameters
   dfm_parameters3d_t parameters3d;
+  // Results split into layers
   dfm_results3d_t results3d;
 } sealock_state_t;
 

From 79b2b1db83b38ab7c591ddbb4a44299b7d4d248c Mon Sep 17 00:00:00 2001
From: Erwin Mulder <160875449+erwin-mulder@users.noreply.github.com>
Date: Tue, 21 Jan 2025 16:08:45 +0100
Subject: [PATCH 14/14] Typo in comment

---
 src/sealock.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/sealock.h b/src/sealock.h
index 0d57366..e9c8b8c 100644
--- a/src/sealock.h
+++ b/src/sealock.h
@@ -92,7 +92,7 @@ typedef struct sealock_state_struct {
   // Volumes
   // Note that the 'to' and 'from' locations may differ
   // so their actual water volumes may be different and
-  // should be trated as such.
+  // should be treated as such.
   dfm_volumes_t from_lake_volumes;
   dfm_volumes_t from_sea_volumes;
   dfm_volumes_t to_lake_volumes;