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