Enable Option for High Res Mesh in Ocean

compass/landice/mesh.py

@@ -265,12 +265,10 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
     low_bed = section.getfloat('low_bed')
     high_bed = section.getfloat('high_bed')
 
-    # convert km to m
-    cull_distance = section.getfloat('cull_distance') * 1.e3
-
     # Cell spacing function based on union of masks
     if section.get('use_bed') == 'True':
         logger.info('Using bed elevation for spacing.')
+
         if flood_fill_iStart is not None and flood_fill_jStart is not None:
             logger.info('calling gridded_flood_fill to find \
                         bedTopography <= low_bed connected to the ocean.')
@@ -290,9 +288,11 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
         k = 0.05  # This works well, but could try other values
         spacing_bed = min_spac + (max_spac - min_spac) / (1.0 + np.exp(
             -k * (bed - np.mean([high_bed, low_bed]))))
-        # We only want bed topography to influence spacing within high_dist_bed
-        # from the ice margin. In the region between high_dist_bed and
-        # low_dist_bed, use a linear ramp to damp influence of bed topo.
+
+        # We only want bed topography to influence spacing within
+        # high_dist_bed from the ice margin. In the region
+        # between high_dist_bed and low_dist_bed, use a linear
+        # ramp to damp influence of bed topo.
         spacing_bed[dist_to_grounding_line >= low_dist_bed] = (
             (1.0 - (dist_to_grounding_line[
                 dist_to_grounding_line >= low_dist_bed] -
@@ -302,6 +302,17 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
                                     low_dist_bed] - low_dist_bed) /
             (high_dist_bed - low_dist_bed) * max_spac)
         spacing_bed[dist_to_grounding_line >= high_dist_bed] = max_spac
+
+        # If max_res_in_ocn is true, use the minimum cell spacing for
+        # all ocean cells. Important for ocean-coupled runs where
+        # resolving fjords and ocean bathymetry is necessary for
+        # thermal forcing parameterizations. Otherwize, telescope out
+        # to coarse resolution with distance from the ice front.
+        if section.get('max_res_in_ocn') == 'True':
+            spacing_bed[np.logical_and(bed < 0, thk == 0)] = min_spac
+            print("Maximizing resolution in ocean. {} ocean cells \
+                  ".format(np.sum(np.logical_and(bed < 0, thk == 0))))
+
         if flood_fill_iStart is not None and flood_fill_jStart is not None:
             spacing_bed[low_bed_mask == 0] = max_spac
             # Do one more flood fill to eliminate isolated pockets
@@ -366,16 +377,30 @@ def set_cell_width(self, section_name, thk, bed=None, vx=None, vy=None,
     for width in [spacing_bed, spacing_speed, spacing_edge, spacing_gl]:
         cell_width = np.minimum(cell_width, width)
 
-    # Set large cell_width in areas we are going to cull anyway (speeds up
-    # whole process). Use 3x the cull_distance to avoid this affecting
-    # cell size in the final mesh. There may be a more rigorous way to set
-    # that distance.
-    if dist_to_edge is not None:
-        mask = np.logical_and(
-            thk == 0.0, dist_to_edge > (3. * cull_distance))
-        logger.info('Setting cell_width in outer regions to max_spac '
-                    f'for {mask.sum()} cells')
-        cell_width[mask] = max_spac
+    if section.get('define_bnds_by_geojson') == 'False':
+        # Set large cell_width in areas we are going to cull anyway (speeds
+        # up whole process). If max_res_in_ocn is True, then use a larger
+        # multiplier to ensure ocean cells are not accidentally coarsened
+        # within the final domain. Otherwise, we can get away with
+        # something smaller, like 3x the cull_distance, to avoid this
+        # affecting the cell size in the final mesh. There may eventually
+        # be a more rigorous way to set this distance. Skip this step if
+        # using geojson to pre-cull domain
+
+        # convert km to m
+        cull_distance = section.getfloat('cull_distance') * 1.e3
+
+        if dist_to_edge is not None:
+            if section.get('max_res_in_ocn') == 'True':
+                mask = np.logical_and(
+                    # Try 20x cull_distance for now
+                    thk == 0.0, dist_to_edge > (20. * cull_distance))
+            else:
+                mask = np.logical_and(
+                    thk == 0.0, dist_to_edge > (3. * cull_distance))
+            logger.info('Setting cell_width in outer regions to max_spac '
+                        f'for {mask.sum()} cells')
+            cell_width[mask] = max_spac
 
     return cell_width
 
@@ -518,7 +543,7 @@ def get_dist_to_edge_and_gl(self, thk, topg, x, y,
 
 
 def build_cell_width(self, section_name, gridded_dataset,
-                     flood_fill_start=[None, None]):
+                     flood_fill_start=[None, None], calc_geom_bnds=True):
     """
     Determine MPAS mesh cell size based on user-defined density function.
 
@@ -543,6 +568,11 @@ def build_cell_width(self, section_name, gridded_dataset,
         fill. Most cases will use ``[None, None]``, which will just start the
         flood fill in the center of the gridded dataset.
 
+    calc_geom_bnds : logical
+        Option to calculate geom_points and geom_edges needed for jigsaw within
+        build_cell_width. Default is to perform calculation, but the user may
+        opt out if these are determined elsewhere (e.g., using a geojson file)
+
     Returns
     -------
     cell_width : numpy.ndarray
@@ -581,15 +611,16 @@ def build_cell_width(self, section_name, gridded_dataset,
 
     f.close()
 
-    # Get bounds defined by user, or use bound of gridded dataset
-    bnds = [np.min(x1), np.max(x1), np.min(y1), np.max(y1)]
-    bnds_options = ['x_min', 'x_max', 'y_min', 'y_max']
-    for index, option in enumerate(bnds_options):
-        bnd = section.get(option)
-        if bnd != 'None':
-            bnds[index] = float(bnd)
+    # If necessary, get bounds defined by user or use bound of gridded dataset
+    if calc_geom_bnds:
+        bnds = [np.min(x1), np.max(x1), np.min(y1), np.max(y1)]
+        bnds_options = ['x_min', 'x_max', 'y_min', 'y_max']
+        for index, option in enumerate(bnds_options):
+            bnd = section.get(option)
+            if bnd != 'None':
+                bnds[index] = float(bnd)
 
-    geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds)
+        geom_points, geom_edges = set_rectangular_geom_points_and_edges(*bnds)
 
     # Remove ice not connected to the ice sheet.
     flood_mask = gridded_flood_fill(thk)
@@ -611,8 +642,12 @@ def build_cell_width(self, section_name, gridded_dataset,
                                 flood_fill_iStart=flood_fill_start[0],
                                 flood_fill_jStart=flood_fill_start[1])
 
-    return (cell_width.astype('float64'), x1.astype('float64'),
-            y1.astype('float64'), geom_points, geom_edges, flood_mask)
+    if not calc_geom_bnds:
+        return (cell_width.astype('float64'), x1.astype('float64'),
+                y1.astype('float64'), flood_mask)
+    else:
+        return (cell_width.astype('float64'), x1.astype('float64'),
+                y1.astype('float64'), geom_points, geom_edges, flood_mask)
 
 
 def build_mali_mesh(self, cell_width, x1, y1, geom_points,