From 5ad0f53073915a70caabc4d9c468b4af72f69eba Mon Sep 17 00:00:00 2001
From: Alice Barthel <abarthel@lanl.gov>
Date: Mon, 30 Oct 2023 18:29:40 -0500
Subject: [PATCH] added linear dz in vertical grid

---
 compass/ocean/vertical/grid_1d/__init__.py  |  8 +++++
 compass/ocean/vertical/grid_1d/linear_dz.py | 38 +++++++++++++++++++++
 2 files changed, 46 insertions(+)
 create mode 100644 compass/ocean/vertical/grid_1d/linear_dz.py

diff --git a/compass/ocean/vertical/grid_1d/__init__.py b/compass/ocean/vertical/grid_1d/__init__.py
index fab57b5f23..89064fb8ec 100644
--- a/compass/ocean/vertical/grid_1d/__init__.py
+++ b/compass/ocean/vertical/grid_1d/__init__.py
@@ -8,6 +8,7 @@
 from compass.ocean.vertical.grid_1d.index_tanh_dz import (
     create_index_tanh_dz_grid,
 )
+from compass.ocean.vertical.grid_1d.linear_dz import create_linear_dz_grid
 from compass.ocean.vertical.grid_1d.tanh_dz import create_tanh_dz_grid
 
 
@@ -32,6 +33,13 @@ def generate_1d_grid(config):
     if grid_type == 'uniform':
         vert_levels = section.getint('vert_levels')
         interfaces = _generate_uniform(vert_levels)
+    elif grid_type == 'linear_dz':
+        vert_levels = section.getint('vert_levels')
+        bottom_depth = section.getfloat('bottom_depth')
+        linear_dz_rate = section.getfloat('linear_dz_rate')
+        interfaces = create_linear_dz_grid(vert_levels,
+                                           bottom_depth,
+                                           linear_dz_rate)
     elif grid_type == 'tanh_dz':
         vert_levels = section.getint('vert_levels')
         min_layer_thickness = section.getfloat('min_layer_thickness')
diff --git a/compass/ocean/vertical/grid_1d/linear_dz.py b/compass/ocean/vertical/grid_1d/linear_dz.py
new file mode 100644
index 0000000000..465fd8bd9a
--- /dev/null
+++ b/compass/ocean/vertical/grid_1d/linear_dz.py
@@ -0,0 +1,38 @@
+import numpy
+import numpy as np
+
+
+def create_linear_dz_grid(num_vert_levels, bottom_depth,
+                          linear_dz_rate):
+    """
+    Creates the linear vertical grid for MPAS-Ocean and
+    writes it to a NetCDF file
+
+    Parameters
+    ----------
+    num_vert_levels : int
+        Number of vertical levels for the grid
+
+    bottom_depth : float
+        bottom depth for the chosen vertical coordinate [m]
+
+    linear_dz_rate : float
+        rate of layer thickness increase (for linear_dz) [m]
+
+    Returns
+    -------
+    interfaces : numpy.ndarray
+        A 1D array of positive depths for layer interfaces in meters
+    """
+
+    nz = num_vert_levels
+    layerThickness = [(bottom_depth / nz) - (np.floor(nz / 2) - k) *
+                      linear_dz_rate for k in np.arange(0, nz)]
+    min_layer_thickness = layerThickness[0]
+    max_layer_thickness = layerThickness[-1]
+    print('Linear dz vertical grid')
+    print(f'min layer thickness: {min_layer_thickness}; '
+          f'max layer thickness {max_layer_thickness} in m;')
+    interfaces = - np.append([0], np.cumsum(layerThickness))
+
+    return interfaces