Make an ISOMIP+ base mesh

isomip_plus_base_mesh.py

#!/usr/bin/env python3

import numpy as np
from mpas_tools.io import write_netcdf
from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.translate import translate


def main():
    """
    Make the base mesh
    """
    lx = 800.0
    ly = 80.0
    buffer = 80.0

    resolution = 2.0

    nx, ny = compute_planar_hex_nx_ny(lx + 2 * buffer, ly + 2 * buffer,
                                      resolution)
    # distance between cells in meters
    dc = 1e3 * resolution

    ds_mesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True,
                                   nonperiodic_y=True)

    translate(mesh=ds_mesh, xOffset=-1e3 * buffer, yOffset=-1e3 * buffer)

    write_netcdf(ds_mesh, 'base_mesh.nc')

        
def compute_planar_hex_nx_ny(lx, ly, resolution):
    """
    Compute number of grid cells in each direction for the uniform, hexagonal
    planar mesh with the given physical sizes and resolution.  The resulting
    ``nx`` and ``ny`` account for the staggered nature of the hexagonal grid
    in the y direction, and are appropriate for passing to
    :py:func:`mpas_tools.planar_hex.make_planar_hex_mesh()`.

    Parameters
    ----------
    lx : float
        The size of the domain in km in the x direction

    ly : float
        The size of the domain in km in the y direction

    resolution : float
        The resolution of the mesh (distance between cell centers) in km

    Returns
    -------
    nx : int
        The number of grid cells in the x direction

    ny : int
        The number of grid cells in the y direction
    """
    # these could be hard-coded as functions of specific supported
    # resolutions but it is preferable to make them algorithmic like here
    # for greater flexibility
    nx = max(2 * int(0.5 * lx / resolution + 0.5), 4)
    # factor of 2/sqrt(3) because of hexagonal mesh
    ny = max(2 * int(0.5 * ly * (2. / np.sqrt(3)) / resolution + 0.5), 4)
    return nx, ny


if __name__ == '__main__':
    main()