Source code for polaris.tasks.ocean.overflow.init

import numpy as np
import xarray as xr
from mpas_tools.mesh.conversion import convert, cull
from mpas_tools.planar_hex import make_planar_hex_mesh

from polaris.mesh.planar import compute_planar_hex_nx_ny
from polaris.ocean.model import OceanIOStep
from polaris.ocean.model.eos import compute_density
from polaris.ocean.vertical import init_vertical_coord


[docs] class Init(OceanIOStep): """ A step for creating a mesh and initial condition for overflow test cases. """
[docs] def __init__(self, component, name='init', indir=None): """ Create the step Parameters ---------- component : polaris.Component The component the step belongs to name : str The name of the step indir : str The name of the directory the task will be set up in """ super().__init__(component=component, name=name, indir=indir)
[docs] def setup(self): super().setup() output_filenames = ['base_mesh.nc', 'culled_mesh.nc', 'init.nc'] model = self.config.get('ocean', 'model') if model == 'mpas-ocean': output_filenames.append('culled_graph.info') for filename in output_filenames: self.add_output_file(filename=filename)
[docs] def run(self): """ Run this step of the test case """ config = self.config logger = self.logger section = config['overflow'] lx = section.getfloat('lx') ly = section.getfloat('ly') resolution = section.getfloat('resolution') nx, ny = compute_planar_hex_nx_ny(lx, ly, resolution) dc = 1e3 * resolution ds_mesh = make_planar_hex_mesh( nx=nx, ny=ny, dc=dc, nonperiodic_x=True, nonperiodic_y=False ) self.write_model_dataset(ds_mesh, 'base_mesh.nc') ds_mesh = cull(ds_mesh, logger=logger) ds_mesh = convert( ds_mesh, graphInfoFileName='culled_graph.info', logger=logger ) self.write_model_dataset(ds_mesh, 'culled_mesh.nc') max_bottom_depth = section.getfloat('max_bottom_depth') shelf_depth = section.getfloat('shelf_depth') x_slope = section.getfloat('x_slope') l_slope = section.getfloat('l_slope') ds = ds_mesh.copy() # Form a continental shelf-like bathymetry ds['bottomDepth'] = shelf_depth + 0.5 * ( max_bottom_depth - shelf_depth ) * ( 1.0 + np.tanh( (ds.xCell - ds.xCell.min() - x_slope * 1.0e3) / (l_slope * 1.0e3) ) ) # ssh is zero ds['ssh'] = xr.zeros_like(ds.xCell) init_vertical_coord(config, ds) # initial temperature is constant except for a block of cold water on # the shelf x_dense = section.getfloat('x_dense') lower_temperature = section.getfloat('lower_temperature') higher_temperature = section.getfloat('higher_temperature') _, x = np.meshgrid(np.zeros(ds.sizes['nVertLevels']), ds.xCell) temperature = np.where( (x - ds.xCell.min().values) < x_dense * 1.0e3, lower_temperature, higher_temperature, ) ds['temperature'] = ( ( 'Time', 'nCells', 'nVertLevels', ), np.expand_dims(temperature, axis=0), ) # initial salinity is constant salinity = section.getfloat('salinity') * np.ones_like(temperature) ds['salinity'] = salinity * xr.ones_like(ds.temperature) density = compute_density(config, temperature, salinity) ds['density'] = ( ( 'Time', 'nCells', 'nVertLevels', ), np.expand_dims(density, axis=0), ) # initial velocity on edges is stationary ds['normalVelocity'] = ( ( 'Time', 'nEdges', 'nVertLevels', ), np.zeros([1, ds.sizes['nEdges'], ds.sizes['nVertLevels']]), ) # Coriolis parameter is zero ds['fCell'] = ( ( 'nCells', 'nVertLevels', ), np.zeros([ds.sizes['nCells'], ds.sizes['nVertLevels']]), ) ds['fEdge'] = ( ( 'nEdges', 'nVertLevels', ), np.zeros([ds.sizes['nEdges'], ds.sizes['nVertLevels']]), ) ds['fVertex'] = ( ( 'nVertices', 'nVertLevels', ), np.zeros([ds.sizes['nVertices'], ds.sizes['nVertLevels']]), ) # finalize and write file self.write_model_dataset(ds, 'init.nc')