diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py
index a82724d44c..36207f8384 100644
--- a/compass/ocean/iceshelf.py
+++ b/compass/ocean/iceshelf.py
@@ -40,40 +40,6 @@ def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask,
     return land_ice_pressure
 
 
-def compute_land_ice_density_from_draft(land_ice_draft, land_ice_thickness,
-                                        floating_mask, ref_density=None):
-    """
-    Compute the spatially-averaged ice density needed to match the ice draft
-
-    Parameters
-    ----------
-    land_ice_draft : xarray.DataArray
-        The ice draft (sea surface height)
-
-    land_ice_thickness: xarray.DataArray
-        The ice thickness
-
-    floating_mask : xarray.DataArray
-        A mask that is 1 where the ice is assumed in hydrostatic equilibrium
-
-    ref_density : float, optional
-        A reference density for seawater displaced by the ice shelf
-
-    Returns
-    -------
-    land_ice_density: float
-        The ice density
-    """
-    if ref_density is None:
-        ref_density = constants['SHR_CONST_RHOSW']
-    land_ice_draft = numpy.where(floating_mask, land_ice_draft, numpy.nan)
-    land_ice_thickness = numpy.where(floating_mask, land_ice_thickness,
-                                     numpy.nan)
-    land_ice_density = \
-        numpy.nanmean(-ref_density * land_ice_draft / land_ice_thickness)
-    return land_ice_density
-
-
 def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask,
                                          ref_density=None):
     """
diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py
index a81da8832b..1103579d9d 100644
--- a/compass/ocean/tests/isomip_plus/initial_state.py
+++ b/compass/ocean/tests/isomip_plus/initial_state.py
@@ -8,7 +8,7 @@
 from mpas_tools.io import write_netcdf
 
 from compass.ocean.iceshelf import (
-    compute_land_ice_density_from_draft,
+    compute_land_ice_draft_from_pressure,
     compute_land_ice_pressure_from_draft,
     compute_land_ice_pressure_from_thickness,
 )
@@ -101,6 +101,7 @@ def _compute_initial_condition(self):
         section = config['isomip_plus']
         min_land_ice_fraction = section.getfloat('min_land_ice_fraction')
         min_ocean_fraction = section.getfloat('min_ocean_fraction')
+        ice_density = section.getfloat('ice_density')
 
         ds_geom = xr.open_dataset('input_geometry_processed.nc')
         ds_mesh = xr.open_dataset('culled_mesh.nc')
@@ -108,6 +109,8 @@ def _compute_initial_condition(self):
         ds = interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction,
                               thin_film_present)
 
+        ds['bottomDepth'] = -ds.bottomDepthObserved
+
         ds['landIceFraction'] = \
             ds['landIceFraction'].expand_dims(dim='Time', axis=0)
         ds['landIceFloatingFraction'] = \
@@ -129,16 +132,17 @@ def _compute_initial_condition(self):
 
         ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.)
 
-        land_ice_draft = ds.ssh
         if thin_film_present:
             land_ice_thickness = ds.landIceThickness
-            land_ice_density = compute_land_ice_density_from_draft(
-                land_ice_draft, land_ice_thickness,
-                floating_mask=ds.landIceFloatingMask)
             land_ice_pressure = compute_land_ice_pressure_from_thickness(
                 land_ice_thickness=land_ice_thickness, modify_mask=ds.ssh < 0.,
-                land_ice_density=land_ice_density)
+                land_ice_density=ice_density)
+            land_ice_draft = compute_land_ice_draft_from_pressure(
+                land_ice_pressure=land_ice_pressure,
+                modify_mask=ds.bottomDepth > 0.)
+            ds['ssh'] = np.maximum(land_ice_draft, -ds.bottomDepth)
         else:
+            land_ice_draft = ds.ssh
             land_ice_pressure = compute_land_ice_pressure_from_draft(
                 land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.)
 
@@ -147,8 +151,6 @@ def _compute_initial_condition(self):
         if self.time_varying_forcing:
             self._write_time_varying_forcing(ds_init=ds)
 
-        ds['bottomDepth'] = -ds.bottomDepthObserved
-
         section = config['isomip_plus']
 
         # Deepen the bottom depth to maintain the minimum water-column
diff --git a/compass/ocean/tests/isomip_plus/isomip_plus.cfg b/compass/ocean/tests/isomip_plus/isomip_plus.cfg
index 355daa71f8..65caf293aa 100644
--- a/compass/ocean/tests/isomip_plus/isomip_plus.cfg
+++ b/compass/ocean/tests/isomip_plus/isomip_plus.cfg
@@ -78,6 +78,9 @@ min_smoothed_draft_mask = 0.01
 # considered a land-ice cell by MPAS-Ocean (landIceMask == 1).
 min_land_ice_fraction = 0.5
 
+# the density of ice prescribed in ISOMIP+
+ice_density = 918
+
 # the initial temperature at the sea surface
 init_top_temp = -1.9
 # the initial temperature at the sea floor