From e8df6203b7fddfbbe5bf5c1c0489745c7fb00199 Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Thu, 15 Aug 2024 03:06:00 -0500
Subject: [PATCH] Set ssh to land-ice draft only in meshes with cavities

Without this fix, meshes without cavities were ending up with a
ring of depressed sea surface height around Antarctica.
---
 compass/ocean/mesh/cull.py             | 6 ++++--
 compass/ocean/mesh/remap_topography.py | 2 --
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py
index 425cfc7cd..d4ab6b004 100644
--- a/compass/ocean/mesh/cull.py
+++ b/compass/ocean/mesh/cull.py
@@ -492,6 +492,8 @@ def _cull_topo(with_cavities, process_count, logger, latitude_threshold,
                                           ds_map_culled_to_base, ds_preserve,
                                           logger, latitude_threshold,
                                           sweep_count)
+        ds_topo['ssh'] = ds_topo['landIceDraftObserved']
+
         write_netcdf(ds_topo, 'topography_with_land_ice_mask.nc')
 
     logger.info('Culling topography')
@@ -534,7 +536,7 @@ def _flood_fill_and_add_land_ice_mask(ds_topo, ds_base_mesh,
         not_preserve = ds.transectCellMasks.sum(dim='nTransects') == 0
 
         for var in ['landIceFracObserved', 'landIcePressureObserved',
-                    'landIceDraftObserved', 'ssh',
+                    'landIceDraftObserved',
                     'landIceGroundedFracObserved',
                     'landIceFloatingFracObserved', 'landIceThkObserved']:
             ds_topo[var] = ds_topo[var].where(not_preserve, 0.0)
@@ -562,7 +564,7 @@ def _flood_fill_and_add_land_ice_mask(ds_topo, ds_base_mesh,
 
     # update land-ice variables and ocean fraction accordingly
     for var in ['landIceFracObserved', 'landIcePressureObserved',
-                'landIceDraftObserved', 'ssh', 'landIceGroundedFracObserved',
+                'landIceDraftObserved', 'landIceGroundedFracObserved',
                 'landIceFloatingFracObserved', 'landIceThkObserved']:
         ds_topo[var] = ds_topo[var].where(land_ice_present, 0.0)
 
diff --git a/compass/ocean/mesh/remap_topography.py b/compass/ocean/mesh/remap_topography.py
index e1c164b93..0e0356220 100644
--- a/compass/ocean/mesh/remap_topography.py
+++ b/compass/ocean/mesh/remap_topography.py
@@ -170,6 +170,4 @@ def run(self):
 
         ds_out['landIceDraftObserved'] = draft
 
-        ds_out['ssh'] = draft
-
         write_netcdf(ds_out, 'topography_remapped.nc')