Floodfill landIceFraction, not landIceMask

During the `cull_mesh` step of creating a global ocean mesh,
if ice-shelf cavities are present, we need to flood fill to make
sure all land ice is connected (no islands of isolated land ice).
This merge changes the flood fill to be on all locations where
`landIceFraction > 0`, not all locations where `landIceFraction > 0.5`.
It then masks all `landIce*` fields to only be nonzero within the
flood-filled region.  This prevents islands of nonzero
`landIceFraction`, `landIcePressure`, `landIceDraft`, etc.

compass/ocean/mesh/cull.py

@@ -488,7 +488,7 @@ def _cull_topo(with_cavities, process_count, logger, latitude_threshold,
     write_netcdf(ds_map_culled_to_base, 'map_culled_to_base.nc')
 
     if with_cavities:
-        _add_land_ice_mask_and_mask_draft(ds_topo, ds_base,
+        _flood_fill_and_add_land_ice_mask(ds_topo, ds_base,
                                           ds_map_culled_to_base, ds_preserve,
                                           logger, latitude_threshold,
                                           sweep_count)
@@ -525,7 +525,7 @@ def _land_mask_from_topo(with_cavities, topo_filename, mask_filename):
     write_netcdf(ds_mask, mask_filename)
 
 
-def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
+def _flood_fill_and_add_land_ice_mask(ds_topo, ds_base_mesh,
                                       ds_map_culled_to_base, ds_perserve,
                                       logger, latitude_threshold, sweep_count):
 
@@ -534,7 +534,8 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
         not_preserve = ds.transectCellMasks.sum(dim='nTransects') == 0
 
         for var in ['landIceFracObserved', 'landIcePressureObserved',
-                    'landIceDraftObserved', 'landIceGroundedFracObserved',
+                    'landIceDraftObserved', 'ssh',
+                    'landIceGroundedFracObserved',
                     'landIceFloatingFracObserved', 'landIceThkObserved']:
             ds_topo[var] = ds_topo[var].where(not_preserve, 0.0)
 
@@ -543,22 +544,30 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
 
     land_ice_frac = ds_topo.landIceFracObserved
 
-    # we want the mask to be 1 where there's at least half land-ice
-    land_ice_mask = xr.where(land_ice_frac > 0.5, 1, 0)
+    land_ice_present = xr.where(land_ice_frac > 0.0, 1, 0)
 
     gf = GeometricFeatures()
     fc_ocean_seed = gf.read(componentName='ocean', objectType='point',
                             tags=['seed_point'])
 
     fc_south_pole_seed = read_feature_collection('south_pole.geojson')
 
-    # flood fill the ice portion to remove isolated land ice
+    # flood fill anywhere land ice is present to remove isolated land ice
 
     ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh,
-                                           daGrow=land_ice_mask,
+                                           daGrow=land_ice_present,
                                            fcSeed=fc_south_pole_seed,
                                            logger=logger)
-    land_ice_mask = ds_mask.cellSeedMask
+    land_ice_present = ds_mask.cellSeedMask
+
+    # update land-ice variables and ocean fraction accordingly
+    for var in ['landIceFracObserved', 'landIcePressureObserved',
+                'landIceDraftObserved', 'ssh', 'landIceGroundedFracObserved',
+                'landIceFloatingFracObserved', 'landIceThkObserved']:
+        ds_topo[var] = ds_topo[var].where(land_ice_present, 0.0)
+
+    land_ice_frac = ds_topo.landIceFracObserved
+    land_ice_mask = xr.where(land_ice_frac > 0.5, 1, 0)
 
     # now, remove land-locked or land-ice-locked cells
 
@@ -600,19 +609,6 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
     region_cell_mask = land_ice_mask.expand_dims(dim='nRegions', axis=1)
     ds_topo['regionCellMasks'] = region_cell_mask
 
-    # we also want to mask out the land-ice draft for cells detatched from the
-    # ice sheet, this time anywhere the land-ice fraction > 0
-    land_ice_draft_mask = xr.where(land_ice_frac > 0.0, 1, 0)
-    ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh,
-                                           daGrow=land_ice_draft_mask,
-                                           fcSeed=fc_south_pole_seed,
-                                           logger=logger)
-
-    land_ice_draft_mask = ds_mask.cellSeedMask
-    ds_topo['landIceDraftMask'] = land_ice_draft_mask
-    ds_topo['landIceDraftObserved'] = (
-        land_ice_draft_mask * ds_topo.landIceDraftObserved)
-
 
 def _land_mask_from_geojson(with_cavities, process_count, logger,
                             mesh_filename, geojson_filename, mask_filename):