diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py
index 8a69efa617..137361f1b1 100644
--- a/compass/ocean/iceshelf.py
+++ b/compass/ocean/iceshelf.py
@@ -11,34 +11,96 @@
 from compass.model import partition, run_model
 
 
-def compute_land_ice_pressure_and_draft(ssh, modify_mask, ref_density):
+def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask,
+                                             land_ice_density=None):
     """
-    Compute the pressure from and overlying ice shelf and the ice-shelf draft
+    Compute the pressure from an overlying ice shelf from ice thickness
 
     Parameters
     ----------
-    ssh : xarray.DataArray
-        The sea surface height (the ice draft)
+    land_ice_thickness: xarray.DataArray
+        The ice thickness
 
     modify_mask : xarray.DataArray
         A mask that is 1 where ``landIcePressure`` can be deviate from 0
 
-    ref_density : float
+    land_ice_density : float, optional
+        A reference density for land ice
+
+    Returns
+    -------
+    land_ice_pressure : xarray.DataArray
+        The pressure from the overlying land ice on the ocean
+    """
+    gravity = constants['SHR_CONST_G']
+    if land_ice_density is None:
+        land_ice_density = constants['SHR_CONST_RHOICE']
+    land_ice_pressure = modify_mask * \
+        numpy.maximum(land_ice_density * gravity * land_ice_thickness, 0.)
+    return land_ice_pressure
+
+
+def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask,
+                                         ref_density=None):
+    """
+    Compute the pressure from an overlying ice shelf from ice draft
+
+    Parameters
+    ----------
+    land_ice_draft : xarray.DataArray
+        The ice draft (sea surface height)
+
+    modify_mask : xarray.DataArray
+        A mask that is 1 where ``landIcePressure`` can be deviate from 0
+
+    ref_density : float, optional
         A reference density for seawater displaced by the ice shelf
 
     Returns
     -------
-    landIcePressure : xarray.DataArray
+    land_ice_pressure : xarray.DataArray
         The pressure from the overlying land ice on the ocean
+    """
+    gravity = constants['SHR_CONST_G']
+    if ref_density is None:
+        ref_density = constants['SHR_CONST_RHOSW']
+    land_ice_pressure = \
+        modify_mask * numpy.maximum(-ref_density * gravity * land_ice_draft,
+                                    0.)
+    return land_ice_pressure
+
+
+def compute_land_ice_draft_from_pressure(land_ice_pressure, modify_mask,
+                                         ref_density=None):
+    """
+    Compute the ice-shelf draft associated with the pressure from an overlying
+    ice shelf
+
+    Parameters
+    ----------
+    land_ice_pressure : xarray.DataArray
+        The pressure from the overlying land ice on the ocean
+
+    modify_mask : xarray.DataArray
+        A mask that is 1 where ``landIcePressure`` can be deviate from 0
 
-    landIceDraft : xarray.DataArray
-        The ice draft, equal to the initial ``ssh``
+    ref_density : float, optional
+        A reference density for seawater displaced by the ice shelf
+
+    Returns
+    -------
+    land_ice_draft : xarray.DataArray
+        The ice draft
     """
     gravity = constants['SHR_CONST_G']
-    landIcePressure = \
-        modify_mask * numpy.maximum(-ref_density * gravity * ssh, 0.)
-    landIceDraft = ssh
-    return landIcePressure, landIceDraft
+    if ref_density is None:
+        ref_density = constants['SHR_CONST_RHOSW']
+    land_ice_draft_array = \
+        - (modify_mask.values *
+           land_ice_pressure.values / (ref_density * gravity))
+    land_ice_draft = xarray.DataArray(data=land_ice_draft_array,
+                                      dims=(land_ice_pressure.dims))
+    return land_ice_draft
 
 
 def adjust_ssh(variable, iteration_count, step, update_pio=True,
diff --git a/compass/ocean/suites/isomip_rk4.txt b/compass/ocean/suites/isomip_rk4.txt
new file mode 100644
index 0000000000..ca08984bb5
--- /dev/null
+++ b/compass/ocean/suites/isomip_rk4.txt
@@ -0,0 +1,11 @@
+# ocean/isomip_plus/planar/5km/z-star/Ocean0
+# ocean/isomip_plus/planar/5km/z-star/time_varying_Ocean0
+# ocean/isomip_plus/planar/5km/z-star/thin_film_Ocean0
+# ocean/isomip_plus/planar/5km/z-star/thin_film_time_varying_Ocean0
+
+#ocean/isomip_plus/planar/5km/sigma/thin_film_Ocean0
+# ocean/isomip_plus/planar/5km/sigma/Ocean0
+# ocean/isomip_plus/planar/5km/sigma/thin_film_time_varying_Ocean0
+ocean/isomip_plus/planar/5km/sigma/thin_film_drying_Ocean0
+ocean/isomip_plus/planar/5km/sigma/thin_film_wetting_Ocean0
+ocean/isomip_plus/planar/5km/sigma/thin_film_tidal_forcing_Ocean0
diff --git a/compass/ocean/tests/ice_shelf_2d/initial_state.py b/compass/ocean/tests/ice_shelf_2d/initial_state.py
index d245f3c9c8..cb76da4fbf 100644
--- a/compass/ocean/tests/ice_shelf_2d/initial_state.py
+++ b/compass/ocean/tests/ice_shelf_2d/initial_state.py
@@ -5,7 +5,7 @@
 from mpas_tools.mesh.conversion import convert, cull
 from mpas_tools.planar_hex import make_planar_hex_mesh
 
-from compass.ocean.iceshelf import compute_land_ice_pressure_and_draft
+from compass.ocean.iceshelf import compute_land_ice_pressure_from_draft
 from compass.ocean.vertical import init_vertical_coord
 from compass.step import Step
 
@@ -103,8 +103,10 @@ def run(self):
         landIceFloatingMask = landIceMask.copy()
 
         ref_density = constants['SHR_CONST_RHOSW']
-        landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft(
-            ssh=ds.ssh, modify_mask=modify_mask, ref_density=ref_density)
+        landIceDraft = ds.ssh
+        landIcePressure = compute_land_ice_pressure_from_draft(
+            land_ice_draft=landIceDraft, modify_mask=modify_mask,
+            ref_density=ref_density)
 
         salinity = surface_salinity + ((bottom_salinity - surface_salinity) *
                                        (ds.zMid / (-bottom_depth)))
diff --git a/compass/ocean/tests/isomip_plus/geom.py b/compass/ocean/tests/isomip_plus/geom.py
index ce761a1582..3c52630511 100755
--- a/compass/ocean/tests/isomip_plus/geom.py
+++ b/compass/ocean/tests/isomip_plus/geom.py
@@ -111,11 +111,12 @@ def interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction, thin_film_present):
                                   ds_geom.landIceGroundedFraction)
 
     # mask the topography to the ocean region before interpolation
-    for var in ['Z_bed', 'Z_ice_draft', 'landIceFraction',
+    for var in ['Z_bed', 'Z_ice_surface', 'Z_ice_draft', 'landIceFraction',
                 'landIceFloatingFraction', 'smoothedDraftMask']:
         ds_geom[var] = ds_geom[var] * ds_geom['oceanFraction']
 
     fields = {'bottomDepthObserved': 'Z_bed',
+              'landIceThickness': 'iceThickness',
               'ssh': 'Z_ice_draft',
               'oceanFracObserved': 'oceanFraction',
               'landIceFraction': 'landIceFraction',
@@ -156,7 +157,7 @@ def _get_geom_fields(ds_geom, ds_mesh, thin_film_present):
     y_cell = ds_mesh.yIsomipCell.values
 
     if thin_film_present:
-        ocean_fraction = - ds_geom['landFraction'] + 1.0
+        ocean_fraction = xarray.where(ds_geom['Z_bed'] > 0., 0., 1.)
     else:
         ocean_fraction = (ds_geom['landIceFloatingFraction'] +
                           ds_geom['openOceanFraction'])
diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py
index 955651040b..2c58ac4e7e 100644
--- a/compass/ocean/tests/isomip_plus/initial_state.py
+++ b/compass/ocean/tests/isomip_plus/initial_state.py
@@ -7,7 +7,11 @@
 from mpas_tools.cime.constants import constants
 from mpas_tools.io import write_netcdf
 
-from compass.ocean.iceshelf import compute_land_ice_pressure_and_draft
+from compass.ocean.iceshelf import (
+    compute_land_ice_draft_from_pressure,
+    compute_land_ice_pressure_from_draft,
+    compute_land_ice_pressure_from_thickness,
+)
 from compass.ocean.tests.isomip_plus.geom import interpolate_geom
 from compass.ocean.tests.isomip_plus.viz.plot import MoviePlotter
 from compass.ocean.vertical import init_vertical_coord
@@ -36,7 +40,7 @@ class InitialState(Step):
         Whether the run includes a thin film below grounded ice
     """
     def __init__(self, test_case, resolution, experiment, vertical_coordinate,
-                 time_varying_forcing, thin_film_present):
+                 time_varying_forcing, thin_film_present, include_viz=False):
         """
         Create the step
 
@@ -61,6 +65,7 @@ def __init__(self, test_case, resolution, experiment, vertical_coordinate,
             Whether the run includes a thin film below grounded ice
         """
         super().__init__(test_case=test_case, name='initial_state')
+        self.include_viz = include_viz
         self.resolution = resolution
         self.experiment = experiment
         self.vertical_coordinate = vertical_coordinate
@@ -84,7 +89,8 @@ def run(self):
         """
         ds, frac = self._compute_initial_condition()
         self._compute_restoring(ds, frac)
-        self._plot(ds)
+        if self.include_viz:
+            self._plot(ds)
 
     def _compute_initial_condition(self):
         config = self.config
@@ -97,6 +103,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')
@@ -104,6 +111,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'] = \
@@ -125,22 +134,45 @@ def _compute_initial_condition(self):
 
         ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.)
 
-        ref_density = constants['SHR_CONST_RHOSW']
-        landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft(
-            ssh=ds.ssh, modify_mask=ds.ssh < 0., ref_density=ref_density)
-
-        ds['landIcePressure'] = landIcePressure
-        ds['landIceDraft'] = landIceDraft
+        if thin_film_present:  # compute full pressure in grounded regions
+            land_ice_thickness = ds.landIceThickness
+            land_ice_pressure_unscaled = \
+                compute_land_ice_pressure_from_thickness(
+                    land_ice_thickness=land_ice_thickness,
+                    modify_mask=ds.ssh < 0.,
+                    land_ice_density=ice_density)
+        else:  # assume floatation pressure everywhere
+            land_ice_draft = ds.ssh
+            land_ice_pressure_unscaled = compute_land_ice_pressure_from_draft(
+                land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.)
 
         if self.time_varying_forcing:
-            self._write_time_varying_forcing(ds_init=ds)
+            scales = config.get('isomip_plus_forcing', 'scales')
+            scales = [float(scale)
+                      for scale in scales.replace(',', ' ').split()]
+            land_ice_pressure = land_ice_pressure_unscaled * scales[0]
+
+        if thin_film_present:
+            modify_mask = ds.bottomDepth > 0.
+            land_ice_draft = compute_land_ice_draft_from_pressure(
+                land_ice_pressure=land_ice_pressure,
+                modify_mask=modify_mask)
+            ds['ssh'] = np.maximum(land_ice_draft, -ds.bottomDepth)
+
+        ds['landIceDraft'] = land_ice_draft
+        # We need to add the time dimension for plotting purposes
+        ds['landIcePressure'] = land_ice_pressure.expand_dims(
+            dim='Time', axis=0)
 
-        ds['bottomDepth'] = -ds.bottomDepthObserved
-
-        section = config['isomip_plus']
+        if self.time_varying_forcing:
+            self._write_time_varying_forcing(
+                ds_init=ds,
+                ice_density=ice_density,
+                land_ice_pressure_unscaled=land_ice_pressure_unscaled)
 
         # Deepen the bottom depth to maintain the minimum water-column
         # thickness
+        section = config['isomip_plus']
         min_column_thickness = section.getfloat('min_column_thickness')
         min_layer_thickness = section.getfloat('min_layer_thickness')
         min_levels = section.getint('minimum_levels')
@@ -149,44 +181,52 @@ def _compute_initial_condition(self):
         min_depth = -ds.ssh + min_column_thickness
         ds['bottomDepth'] = np.maximum(ds.bottomDepth, min_depth)
         print(f'Adjusted bottomDepth for '
-              f'{np.sum(ds.bottomDepth.values<min_depth.values)} cells '
+              f'{np.sum(ds.bottomDepth.values < min_depth.values)} cells '
               f'to achieve minimum column thickness of {min_column_thickness}')
 
+        # Initialize vertical coordinate without time dimension
+        ds = ds.isel(Time=0)
         init_vertical_coord(config, ds)
 
-        max_bottom_depth = -config.getfloat('vertical_grid', 'bottom_depth')
-        frac = (0. - ds.zMid) / (0. - max_bottom_depth)
-
         # compute T, S
         init_top_temp = section.getfloat('init_top_temp')
         init_bot_temp = section.getfloat('init_bot_temp')
         init_top_sal = section.getfloat('init_top_sal')
         init_bot_sal = section.getfloat('init_bot_sal')
-        thin_film_mask = np.logical_and(mask.values,
-                                        np.logical_not(floating_mask.values))
-        # These coefficients are hard-coded as the defaults in the namelist
-        # Note that using the land ice pressure rather than the pressure at
-        # floatation will mean that there is a small amount of cooling from
-        # grounding line retreat. However, the thin film should be thin enough
-        # that this effect isn't signicant.
-        freezing_temp = (6.22e-2 +
-                         -5.63e-2 * init_bot_sal +
-                         -7.43e-8 * landIcePressure +
-                         -1.74e-10 * landIcePressure * init_bot_sal)
-        _, thin_film_temp = np.meshgrid(ds.refZMid, freezing_temp)
-        _, thin_film_mask = np.meshgrid(ds.refZMid, thin_film_mask)
-        thin_film_temp = np.expand_dims(thin_film_temp, axis=0)
+
         if self.vertical_coordinate == 'single_layer':
-            ds['temperature'] = init_bot_temp * xr.ones_like(frac)
-            ds['salinity'] = init_bot_sal * xr.ones_like(frac)
+            # Initialize constant T,S
+            ds['temperature'] = init_bot_temp * xr.ones_like(ds.zmid)
+            ds['salinity'] = init_bot_sal * xr.ones_like(ds.zmid)
         else:
+            # Initialize T,S as linear functions with max depth
+            max_bottom_depth = -config.getfloat('vertical_grid',
+                                                'bottom_depth')
+            frac = (0. - ds.zMid) / (0. - max_bottom_depth)
             ds['temperature'] = \
                 (1.0 - frac) * init_top_temp + frac * init_bot_temp
             ds['salinity'] = \
                 (1.0 - frac) * init_top_sal + frac * init_bot_sal
-        # for thin film cells, set temperature to freezing point
-        ds['temperature'] = xr.where(thin_film_mask, thin_film_temp,
-                                     ds.temperature)
+
+        if thin_film_present:
+            # for thin film cells, set temperature to freezing point
+            # TODO consider setting salinity to 0
+            thin_film_mask = np.logical_and(
+                mask.values, np.logical_not(floating_mask.values))
+            # These coefficients are hard-coded as the defaults in the namelist
+            # Note that using the land ice pressure rather than the pressure at
+            # floatation will mean that there is a small amount of cooling from
+            # grounding line retreat. However, the thin film should be thin
+            # enough that this effect isn't signicant.
+            freezing_temp = (6.22e-2 +
+                             -5.63e-2 * init_bot_sal +
+                             -7.43e-8 * land_ice_pressure +
+                             -1.74e-10 * land_ice_pressure * init_bot_sal)
+            _, thin_film_temp = np.meshgrid(ds.refZMid, freezing_temp)
+            _, thin_film_mask = np.meshgrid(ds.refZMid, thin_film_mask)
+            thin_film_temp = np.expand_dims(thin_film_temp, axis=0)
+            ds['temperature'] = xr.where(thin_film_mask, thin_film_temp,
+                                         ds.temperature)
 
         # compute coriolis
         coriolis_parameter = section.getfloat('coriolis_parameter')
@@ -201,6 +241,7 @@ def _compute_initial_condition(self):
         ds['normalVelocity'] = normalVelocity.expand_dims(dim='Time', axis=0)
 
         write_netcdf(ds, 'initial_state.nc')
+
         return ds, frac
 
     def _plot(self, ds):
@@ -227,58 +268,60 @@ def _plot(self, ds):
                                sectionY=section_y, dsMesh=ds, ds=ds,
                                showProgress=show_progress)
 
-        ds['oceanFracObserved'] = \
-            ds['oceanFracObserved'].expand_dims(dim='Time', axis=0)
-        ds['landIcePressure'] = \
-            ds['landIcePressure'].expand_dims(dim='Time', axis=0)
-        ds['landIceGroundedFraction'] = \
-            ds['landIceGroundedFraction'].expand_dims(dim='Time', axis=0)
-        ds['bottomDepth'] = ds['bottomDepth'].expand_dims(dim='Time', axis=0)
-        ds['totalColThickness'] = ds['ssh']
-        ds['totalColThickness'].values = \
-            ds['layerThickness'].sum(dim='nVertLevels')
+        totalColThickness = ds.layerThickness.sum(dim='nVertLevels')
         tol = 1e-10
-        plotter.plot_horiz_series(ds.landIceMask,
+        plotter.plot_horiz_series(ds.landIceMask.expand_dims(
+                                  dim='Time', axis=0),
                                   'landIceMask', 'landIceMask',
                                   True)
-        plotter.plot_horiz_series(ds.landIceFloatingMask,
+        plotter.plot_horiz_series(ds.landIceFloatingMask.expand_dims(
+                                  dim='Time', axis=0),
                                   'landIceFloatingMask', 'landIceFloatingMask',
                                   True)
-        plotter.plot_horiz_series(ds.landIcePressure,
+        plotter.plot_horiz_series(ds.landIcePressure.expand_dims(
+                                  dim='Time', axis=0),
                                   'landIcePressure', 'landIcePressure',
-                                  True, vmin=1, vmax=1e6, cmap_scale='log')
-        plotter.plot_horiz_series(ds.ssh,
-                                  'ssh', 'ssh',
+                                  True, vmin=1e5, vmax=1e7, cmap_scale='log')
+        plotter.plot_horiz_series(ds.landIceThickness.expand_dims(
+                                  dim='Time', axis=0),
+                                  'landIceThickness', 'landIceThickness',
+                                  True, vmin=0, vmax=1e3)
+        plotter.plot_horiz_series(ds.ssh, 'ssh', 'ssh',
                                   True, vmin=-700, vmax=0)
-        plotter.plot_horiz_series(ds.bottomDepth,
+        bottomDepth = ds.bottomDepth.expand_dims(dim='Time', axis=0)
+        plotter.plot_horiz_series(bottomDepth,
                                   'bottomDepth', 'bottomDepth',
                                   True, vmin=0, vmax=700)
-        plotter.plot_horiz_series(ds.ssh + ds.bottomDepth,
+        plotter.plot_horiz_series(ds.ssh + bottomDepth,
                                   'H', 'H', True,
                                   vmin=min_column_thickness + tol, vmax=700,
                                   cmap_set_under='r', cmap_scale='log')
-        plotter.plot_horiz_series(ds.totalColThickness,
+        plotter.plot_horiz_series(totalColThickness,
                                   'totalColThickness', 'totalColThickness',
                                   True, vmin=min_column_thickness + 1e-10,
                                   vmax=700, cmap_set_under='r')
-        plotter.plot_horiz_series(ds.landIceFraction,
+        plotter.plot_horiz_series(ds.landIceFraction.expand_dims(
+                                  dim='Time', axis=0),
                                   'landIceFraction', 'landIceFraction',
                                   True, vmin=0 + tol, vmax=1 - tol,
                                   cmap='cmo.balance',
                                   cmap_set_under='k', cmap_set_over='r')
-        plotter.plot_horiz_series(ds.landIceFloatingFraction,
+        plotter.plot_horiz_series(ds.landIceFloatingFraction.expand_dims(
+                                  dim='Time', axis=0),
                                   'landIceFloatingFraction',
                                   'landIceFloatingFraction',
                                   True, vmin=0 + tol, vmax=1 - tol,
                                   cmap='cmo.balance',
                                   cmap_set_under='k', cmap_set_over='r')
-        plotter.plot_horiz_series(ds.landIceGroundedFraction,
+        plotter.plot_horiz_series(ds.landIceGroundedFraction.expand_dims(
+                                  dim='Time', axis=0),
                                   'landIceGroundedFraction',
                                   'landIceGroundedFraction',
                                   True, vmin=0 + tol, vmax=1 - tol,
                                   cmap='cmo.balance',
                                   cmap_set_under='k', cmap_set_over='r')
-        plotter.plot_horiz_series(ds.oceanFracObserved,
+        plotter.plot_horiz_series(ds.oceanFracObserved.expand_dims(
+                                  dim='Time', axis=0),
                                   'oceanFracObserved', 'oceanFracObserved',
                                   True, vmin=0 + tol, vmax=1 - tol,
                                   cmap='cmo.balance',
@@ -363,51 +406,79 @@ def _compute_restoring(self, ds, frac):
 
         write_netcdf(ds_forcing, 'init_mode_forcing_data.nc')
 
-    def _write_time_varying_forcing(self, ds_init):
+    def _write_time_varying_forcing(self, ds_init, ice_density,
+                                    land_ice_pressure_unscaled):
         """
-        Write time-varying land-ice forcing and update the initial condition
+        Write time-varying land-ice forcing
         """
 
         config = self.config
+        ds_forcing = xr.Dataset()
         dates = config.get('isomip_plus_forcing', 'dates')
-        dates = [date.ljust(64) for date in dates.replace(',', ' ').split()]
+        dates = [date.ljust(64)
+                 for date in dates.replace(',', ' ').split()]
         scales = config.get('isomip_plus_forcing', 'scales')
-        scales = [float(scale) for scale in scales.replace(',', ' ').split()]
-
-        ds_out = xr.Dataset()
-        ds_out['xtime'] = ('Time', dates)
-        ds_out['xtime'] = ds_out.xtime.astype('S')
-
-        landIceDraft = list()
-        landIcePressure = list()
-        landIceFraction = list()
-        landIceFloatingFraction = list()
-
-        for scale in scales:
-            landIceDraft.append(scale * ds_init.landIceDraft)
-            landIcePressure.append(scale * ds_init.landIcePressure)
-            landIceFraction.append(ds_init.landIceFraction)
-            landIceFloatingFraction.append(ds_init.landIceFloatingFraction)
-
-        ds_out['landIceDraftForcing'] = xr.concat(landIceDraft, 'Time')
-        ds_out.landIceDraftForcing.attrs['units'] = 'm'
-        ds_out.landIceDraftForcing.attrs['long_name'] = \
+        scales = [float(scale)
+                  for scale in scales.replace(',', ' ').split()]
+
+        # The initial condition already has the first scale value applied
+        land_ice_pressure_forcing = ds_init.landIcePressure.copy()
+        land_ice_fraction_forcing = ds_init.landIceFraction.copy()
+        land_ice_floating_fraction_forcing = \
+            ds_init.landIceFloatingFraction.copy()
+
+        # We add additional time slices for the remaining scale values
+        for scale in scales[1:]:
+            land_ice_pressure_forcing = xr.concat(
+                [land_ice_pressure_forcing,
+                 scale * land_ice_pressure_unscaled],
+                'Time')
+            land_ice_fraction_forcing = xr.concat(
+                [land_ice_fraction_forcing, ds_init.landIceFraction],
+                'Time')
+            # Since floating fraction does not change, none of the thin film
+            # cases allow for the area undergoing melting to change
+            land_ice_floating_fraction_forcing = xr.concat(
+                [land_ice_floating_fraction_forcing,
+                 ds_init.landIceFloatingFraction],
+                'Time')
+
+        if self.thin_film_present:
+            # Set the maximum land ice draft in grounded regions
+            land_ice_draft_forcing = compute_land_ice_draft_from_pressure(
+                land_ice_pressure=land_ice_pressure_forcing,
+                modify_mask=ds_init.bottomDepth > 0.)
+            land_ice_draft_forcing = np.maximum(land_ice_draft_forcing,
+                                                -ds_init.bottomDepth)
+        else:
+            # Just scale draft in the same manner as pressure
+            land_ice_draft_forcing = ds_init.landIceDraft.copy() * scales[0]
+            for scale in scales[1:]:
+                land_ice_draft_forcing = xr.concat(
+                    [land_ice_draft_forcing,
+                     scale * land_ice_pressure_unscaled],
+                    'Time')
+
+        ds_forcing['xtime'] = xr.DataArray(data=dates,
+                                           dims=('Time')).astype('S')
+        ds_forcing['landIceDraftForcing'] = land_ice_draft_forcing
+        ds_forcing.landIceDraftForcing.attrs['units'] = 'm'
+        ds_forcing.landIceDraftForcing.attrs['long_name'] = \
             'The approximate elevation of the land ice-ocean interface'
-        ds_out['landIcePressureForcing'] = \
-            xr.concat(landIcePressure, 'Time')
-        ds_out.landIcePressureForcing.attrs['units'] = 'm'
-        ds_out.landIcePressureForcing.attrs['long_name'] = \
-            'Pressure from the weight of land ice at the ice-ocean interface'
-        ds_out['landIceFractionForcing'] = \
-            xr.concat(landIceFraction, 'Time')
-        ds_out.landIceFractionForcing.attrs['long_name'] = \
+        ds_forcing['landIcePressureForcing'] = land_ice_pressure_forcing
+        ds_forcing.landIcePressureForcing.attrs['units'] = 'm'
+        ds_forcing.landIcePressureForcing.attrs['long_name'] = \
+            'Pressure from the weight of land ice at the ice-ocean ' \
+            'interface'
+        ds_forcing['landIceFractionForcing'] = land_ice_fraction_forcing
+        ds_forcing.landIceFractionForcing.attrs['long_name'] = \
             'The fraction of each cell covered by land ice'
-        ds_out['landIceFloatingFractionForcing'] = \
-            xr.concat(landIceFloatingFraction, 'Time')
-        ds_out.landIceFloatingFractionForcing.attrs['long_name'] = \
+        ds_forcing['landIceFloatingFractionForcing'] = \
+            land_ice_floating_fraction_forcing
+        ds_forcing.landIceFloatingFractionForcing.attrs['long_name'] = \
             'The fraction of each cell covered by floating land ice'
-        write_netcdf(ds_out, 'land_ice_forcing.nc')
 
-        ds_init['landIceDraft'] = scales[0] * ds_init.landIceDraft
-        ds_init['ssh'] = ds_init.landIceDraft
-        ds_init['landIcePressure'] = scales[0] * ds_init.landIcePressure
+        ds_forcing.encoding['unlimited_dims'] = {'Time'}
+        # write_netcdf is not used here because it does not yet support
+        # multiple time levels
+        ds_forcing.to_netcdf('land_ice_forcing.nc')
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
diff --git a/compass/ocean/tests/isomip_plus/isomip_plus_test.py b/compass/ocean/tests/isomip_plus/isomip_plus_test.py
index cf17731a01..b701816661 100644
--- a/compass/ocean/tests/isomip_plus/isomip_plus_test.py
+++ b/compass/ocean/tests/isomip_plus/isomip_plus_test.py
@@ -140,7 +140,8 @@ def __init__(self, test_group, resolution, experiment,
                          experiment=experiment,
                          vertical_coordinate=vertical_coordinate,
                          time_varying_forcing=time_varying_forcing,
-                         thin_film_present=thin_film_present))
+                         thin_film_present=thin_film_present,
+                         include_viz=thin_film_present))
         self.add_step(
             SshAdjustment(test_case=self, resolution=resolution,
                           vertical_coordinate=vertical_coordinate,
diff --git a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust
index 76c2443f7c..047f5db62e 100644
--- a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust
+++ b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust
@@ -1,7 +1,17 @@
 config_time_integrator='RK4'
-config_dt='00:00:10'
+config_thickness_flux_type='upwind'
+config_land_ice_flux_mode='pressure_only'
 config_use_wetting_drying=.true.
+config_zero_drying_velocity_ramp=.true.
+config_zero_drying_velocity_ramp_hmax=1e1
+config_use_ssh_gradient_wetting_drying=.false.
 config_prevent_drying=.true.
 config_drying_min_cell_height=1.0e-3
 config_zero_drying_velocity = .true.
-config_thickness_flux_type='upwind'
+config_use_land_ice_ssh_adjustment = .false.
+config_use_land_ice_forcing=.false.
+config_land_ice_dynamic_height=3.1e-3
+config_use_land_ice_pressure_relaxation = .true.
+config_land_ice_pressure_increment = 10.
+config_use_land_ice_overburden_height=.true.
+config_land_ice_overburden_height=-3e-3
diff --git a/compass/ocean/tests/isomip_plus/process_geom.py b/compass/ocean/tests/isomip_plus/process_geom.py
index 1293cc78cd..3f07250b11 100644
--- a/compass/ocean/tests/isomip_plus/process_geom.py
+++ b/compass/ocean/tests/isomip_plus/process_geom.py
@@ -171,10 +171,18 @@ def _smooth_geometry(land_fraction, ds, filter_sigma, threshold=0.01):
                                       mode='constant', cval=0.)
         bed[mask] /= smoothed_mask[mask]
 
+        # We only use the ice surface where the ocean is present to compute
+        # land ice pressure
+        ice_thickness = filters.gaussian_filter(
+            ds.iceThickness * ocean_fraction, filter_sigma, mode='constant',
+            cval=0.)
+        ice_thickness[mask] /= smoothed_mask[mask]
+
         smoothed_draft_mask = filters.gaussian_filter(
             ds.landIceFloatingFraction, filter_sigma, mode='constant', cval=0.)
         smoothed_draft_mask[mask] /= smoothed_mask[mask]
 
         ds['Z_ice_draft'] = (('y', 'x'), draft)
         ds['Z_bed'] = (('y', 'x'), bed)
+        ds['iceThickness'] = (('y', 'x'), ice_thickness)
         ds['smoothedDraftMask'] = (('y', 'x'), smoothed_draft_mask)
diff --git a/compass/ocean/tests/isomip_plus/streams.forward.template b/compass/ocean/tests/isomip_plus/streams.forward.template
index 4eff72b570..e6d26d6f52 100644
--- a/compass/ocean/tests/isomip_plus/streams.forward.template
+++ b/compass/ocean/tests/isomip_plus/streams.forward.template
@@ -42,6 +42,7 @@
     <var name="density"/>
     <var name="potentialDensity"/>
     <var name="ssh"/>
+    <var name="surfacePressure"/>
     <var name="pressure"/>
     <var name="normalVelocity"/>
     <var name="velocityX"/>
@@ -64,6 +65,7 @@
     <var name="topDragMagnitude"/>
     <var name="accumulatedLandIceMass"/>
     <var name="accumulatedLandIceHeat"/>
+    <var name="wettingVelocityFactor"/>
 </stream>
 
 <stream name="output"
@@ -85,14 +87,15 @@
     <var name="layerThickness"/>
     <var name="zMid"/>
     <var name="ssh"/>
+    <var name="surfacePressure"/>
     <var name="landIcePressure"/>
     <var name="seaIcePressure"/>
-    <var name="landIcePressure"/>
     <var name="atmosphericPressure"/>
     <var name="surfaceStressMagnitude"/>
     <var name="surfaceThicknessFlux"/>
     <var name="velocityX"/>
     <var name="velocityY"/>
+    <var name="wettingVelocityFactor"/>
 </stream>
 
 
diff --git a/compass/ocean/tests/isomip_plus/viz/__init__.py b/compass/ocean/tests/isomip_plus/viz/__init__.py
index 5b9d745e3a..5bca3c48f4 100644
--- a/compass/ocean/tests/isomip_plus/viz/__init__.py
+++ b/compass/ocean/tests/isomip_plus/viz/__init__.py
@@ -99,6 +99,27 @@ def run(self):
                                   True, vmin=min_column_thickness + 1e-10,
                                   vmax=700, cmap_set_under='r',
                                   cmap_scale='log')
+        plotter.plot_horiz_series(dsOut.surfacePressure,
+                                  'surfacePressure', 'surfacePressure',
+                                  True, vmin=1e5, vmax=1e7, cmap_scale='log')
+        delSurfacePressure = dsOut.landIcePressure - dsOut.surfacePressure
+        plotter.plot_horiz_series(delSurfacePressure,
+                                  'delSurfacePressure', 'delSurfacePressure',
+                                  True, vmin=1e5, vmax=1e7, cmap_scale='log')
+        # Cannot currently plot because it is on edges
+        for t in range(dsOut.sizes['Time']):
+            fig = plt.figure()
+            sc = plt.scatter(
+                dsOut.xEdge, dsOut.yEdge, 10,
+                dsOut.wettingVelocityFactor.isel(nVertLevels=0, Time=t),
+                vmin=0, vmax=1, cmap='cmo.thermal')
+            current_cmap = sc.get_cmap()
+            current_cmap.set_under('k')
+            current_cmap.set_over('r')
+            plt.colorbar()
+            fig.gca().set_aspect('equal')
+            plt.savefig(f'wettingVelocityFactor_{t:02g}.png')
+            plt.close(fig)
 
         if 'tidal' in expt:
             delssh = dsOut.ssh - dsOut.ssh[0, :]
@@ -111,10 +132,10 @@ def run(self):
         wct_thin = wct[:, idx_thin]
         wct_mean = wct_thin.mean(dim='nCells').values
         time = dsOut.daysSinceStartOfSim.values
-        fig = plt.figure()
+        plt.figure()
         plt.plot(time, wct_mean, '.')
-        fig.set_xlabel('Time (days)')
-        fig.set_ylabel('Mean thickness of thin film (m)')
+        plt.xlabel('Time (days)')
+        plt.ylabel('Mean thickness of thin film (m)')
         plt.savefig('wct_thin_t.png')
         plt.close()
 
diff --git a/compass/ocean/tests/isomip_plus/viz/plot.py b/compass/ocean/tests/isomip_plus/viz/plot.py
index e231443a5c..02e3471261 100644
--- a/compass/ocean/tests/isomip_plus/viz/plot.py
+++ b/compass/ocean/tests/isomip_plus/viz/plot.py
@@ -84,6 +84,8 @@ def plot_melt_time_series(self, sshMax=None):
         mean thermal driving, mean friction velocity
         """
 
+        if 'timeMonthly_avg_landIceFreshwaterFlux' not in self.ds.keys():
+            return
         rho_fw = 1000.
         secPerYear = 365 * 24 * 60 * 60
 
@@ -239,10 +241,12 @@ def __init__(self, inFolder, streamfunctionFolder, outFolder, expt,
         self.sectionY = sectionY
         self.showProgress = showProgress
 
+        if 'Time' in dsMesh.dims:
+            dsMesh = dsMesh.isel(Time=0)
         self.dsMesh = dsMesh
         self.ds = ds
 
-        landIceMask = self.dsMesh.landIceMask.isel(Time=0) > 0
+        landIceMask = self.dsMesh.landIceMask > 0
         self.oceanMask = self.dsMesh.maxLevelCell - 1 >= 0
         self.cavityMask = numpy.logical_and(self.oceanMask, landIceMask)
 
@@ -349,6 +353,8 @@ def plot_melt_rates(self, vmin=-100., vmax=100.):
         vmin, vmax : float, optional
             The minimum and maximum values for the colorbar
         """
+        if 'timeMonthly_avg_landIceFreshwaterFlux' not in self.ds.keys():
+            return
         rho_fw = 1000.
         secPerYear = 365 * 24 * 60 * 60
 
@@ -367,6 +373,8 @@ def plot_ice_shelf_boundary_variables(self):
         ice
         """
 
+        if 'timeMonthly_avg_landIceFreshwaterFlux' not in self.ds.keys():
+            return
         self.plot_horiz_series(self.ds.timeMonthly_avg_landIceHeatFlux,
                                'heat flux from ocean to ice-ocean interface',
                                prefix='oceanHeatFlux',
@@ -508,7 +516,9 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain,
             The time indices at which to plot. If not provided, set to all.
         """
 
-        nTime = self.ds.sizes['Time']
+        if 'Time' not in da.dims:
+            da = da.expand_dims(dim='Time', axis=0)
+        nTime = da.sizes['Time']
         if self.showProgress:
             widgets = ['plotting {}: '.format(nameInTitle),
                        progressbar.Percentage(), ' ',
diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst
index d7f0d3e05d..c45836a79d 100644
--- a/docs/developers_guide/ocean/api.rst
+++ b/docs/developers_guide/ocean/api.rst
@@ -1037,7 +1037,9 @@ ocean framework
 
    haney.compute_haney_number
 
-   iceshelf.compute_land_ice_pressure_and_draft
+   iceshelf.compute_land_ice_draft_from_pressure
+   iceshelf.compute_land_ice_pressure_from_draft
+   iceshelf.compute_land_ice_pressure_from_thickness
    iceshelf.adjust_ssh
 
    particles.write
diff --git a/docs/users_guide/ocean/test_groups/isomip_plus.rst b/docs/users_guide/ocean/test_groups/isomip_plus.rst
index 567fa600b5..e53a4940b5 100644
--- a/docs/users_guide/ocean/test_groups/isomip_plus.rst
+++ b/docs/users_guide/ocean/test_groups/isomip_plus.rst
@@ -125,6 +125,9 @@ The ``isomip_plus`` test cases share the following config options:
     # 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