diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index dd218f789e4f..b5ca6388df9b 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -1339,6 +1339,10 @@
 					description="If true, use take into account the ssh gradient across previously flux-limited edges in the wetting and drying algorithm"
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_land_ice_overburden_height" type="real" default_value="0.0" units="m"
+					description="The amount of pressure overburden allowed in grounded regions. If ice must float at level of critical thickness, set equal to -1 * critical thickness"
+					possible_values="Any real"
+		/>
 	</nml_record>
 	<nml_record name="ocean_constants">
 		<nml_option name="config_density0" type="real" default_value="1026.0" units="kg m^-3"
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
index a86c24386c8f..1963fc885d1a 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
@@ -76,6 +76,7 @@ module ocn_diagnostics
    !
    !--------------------------------------------------------------------
 
+   logical :: landIcePressureInitialized
    integer :: ke_cell_flag, ke_vertex_flag
    real (kind=RKIND) ::  fCoef
    real (kind=RKIND) ::  landIceTopDragCoeff
@@ -806,6 +807,33 @@ subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, deltaSSH, &
          !$acc exit data copyout(wctEdge) delete(baroclinicThickness, deltaSSH, btrVelocity)
 #endif
 
+      case (thickEdgeFluxConstant)
+
+         ! Use centered (mean) thickness as flux value
+
+         ! Compute mean layer thickness at edge
+#ifdef MPAS_OPENACC
+         !$acc enter data copyin(baroclinicThickness, deltaSSH, wctEdge)
+         !$acc parallel loop &
+         !$acc    present(baroclinicThickness, deltaSSH, wctEdge, cellsOnEdge) &
+         !$acc    private(iEdge, cell1, cell2)
+#else
+         !$omp parallel
+         !$omp do schedule(runtime)
+#endif
+         do iEdge = 1, nEdgesAll
+            if ( maxLevelEdgeTop(iEdge) == 0 ) cycle
+            wctEdge(iEdge) = baroclinicThickness(iEdge)
+         end do
+#ifndef MPAS_OPENACC
+         !$omp end do
+         !$omp end parallel
+#endif
+
+#ifdef MPAS_OPENACC
+         !$acc exit data copyout(wctEdge) delete(baroclinicThickness, deltaSSH)
+#endif
+
       case default
          ! Should have been caught on init
          call mpas_log_write('Thickness flux option unknown', &
@@ -2594,6 +2622,8 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
 
       integer :: iCell ! cell loop index
 
+      real (kind=RKIND) :: h_crit, ssh_crit, landIcePressureMax
+      real (kind=RKIND), dimension(:), allocatable :: landIcePressureApplied
       real (kind=RKIND), dimension(:), pointer :: &
          frazilSurfacePressure, &! frazil ice pressure from forcing
          landIcePressure         ! land ice pressure
@@ -2609,6 +2639,13 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
       call mpas_pool_get_array(forcingPool, 'landIcePressure', &
                                              landIcePressure)
 
+      allocate(landIcePressureApplied(nCellsAll))
+      ! First initialize landIcePressure unchanged, then modify in grounded regions
+      if (landIcePressureOn .and. .not. landIcePressureInitialized) then
+         landIcePressureApplied(:) = landIcePressure(:)
+         landIcePressureInitialized = .true.
+      endif
+
       ! Compute pressure for generalized coordinates.
       ! Pressure at top surface may be due to atmospheric pressure,
       ! sea ice, frazil ice or the weight of an ice shelf
@@ -2652,6 +2689,31 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
       ! Add land ice pressure if needed
 
       if (landIcePressureOn) then
+#ifdef MPAS_OPENACC
+         !$acc parallel loop &
+         !$acc    present(surfacePressure)
+         !! need this eventually, but there is currently
+         !! an issue, probably with the re-retrieval of pointer
+         !!acc    present(surfacePressure, landIcePressure)
+#else
+         !$omp parallel
+         !$omp do schedule(runtime) private(h_crit, ssh_crit, landIcePressureMax)
+#endif
+         do iCell = 1, nCellsAll
+            ! the maximum pressure is the one that results in ssh at the critical
+            ! thickness (gravity is positive), but we also allow some amount of overburden
+            ! because rho may deviate significantly from rho_sw
+            h_crit = config_drying_min_cell_height * (maxLevelCell(iCell) - minLevelCell(iCell) + 1)
+            ssh_crit = bottomDepth(iCell) - h_crit + config_land_ice_overburden_height
+            landIcePressureMax = gravity * rho_sw * ssh_crit
+            landIcePressureApplied(iCell) = &
+               min(landIcePressureMax, landIcePressure(iCell))
+         end do
+#ifndef M_OPENACC
+         !$omp end do
+         !$omp end parallel
+#endif
+
 #ifdef MPAS_OPENACC
          !$acc parallel loop &
          !$acc    present(surfacePressure)
@@ -2663,7 +2725,7 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
 #endif
          do iCell = 1, nCellsAll
             surfacePressure(iCell) = surfacePressure(iCell) &
-                                   + landIcePressure(iCell)
+                                   + landIcePressureApplied(iCell)
          end do
 #ifndef MPAS_OPENACC
          !$omp end do
@@ -4714,6 +4776,8 @@ subroutine ocn_diagnostics_init(domain, err)!{{{
                  & 'Set config_disable_vel_hadv to true for a linear test case. ', &
                  MPAS_LOG_CRIT)
          end if
+         call mpas_log_write('Thickness flux option set to ' //&
+              & trim(config_thickness_flux_type))
       elseif (trim(config_thickness_flux_type) == 'upwind') then
          thickEdgeFluxChoice = thickEdgeFluxUpwind
          call mpas_log_write('Thickness flux option set to ' //&
@@ -4750,6 +4814,12 @@ subroutine ocn_diagnostics_init(domain, err)!{{{
         call mpas_log_write("config_sgr_flux_vertical_location not one of 'bottom', 'uniform', 'top'.", MPAS_LOG_CRIT)
       end if
 
+      if (config_do_restart) then
+         landIcePressureInitialized = .true.
+      else
+         landIcePressureInitialized = .false.
+      endif
+
     end subroutine ocn_diagnostics_init!}}}
 
 !***********************************************************************
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
index abdc7c1ff9b4..f6873b0df2b8 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
@@ -113,6 +113,7 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:), pointer :: surfaceBuoyancyForcing
    real (kind=RKIND), dimension(:), pointer :: surfaceFrictionVelocity
    real (kind=RKIND), dimension(:), pointer :: surfacePressure
+   real (kind=RKIND), dimension(:), pointer :: landIcePressureApplied
 
    real (kind=RKIND), dimension(:,:), pointer :: &
       transportVelocityX, transportVelocityY, transportVelocityZ, transportVelocityZonal, &
@@ -714,6 +715,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       call mpas_pool_get_array(diagnosticsPool, 'dThreshMLD', dThreshMLD)
       call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', &
                 surfacePressure)
+      call mpas_pool_get_array(diagnosticsPool, 'landIcePressureApplied', &
+                landIcePressureApplied)
 
       ! Semi-implicit Array Pointer retrievals
       if (trim(config_time_integrator) == 'split_implicit') then
@@ -966,6 +969,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       !$acc                   thermExpCoeff,                           &
       !$acc                   sfcFlxAttCoeff,                          &
       !$acc                   surfacePressure,                         &
+      !$acc                   landIcePressureApplied,                  &
       !$acc                   potentialDensity,                        &
       !$acc                   displacedDensity,                        &
       !$acc                   boundaryLayerDepth,                      &
@@ -1228,6 +1232,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
       !$acc                   thermExpCoeff,                           &
       !$acc                   sfcFlxAttCoeff,                          &
       !$acc                   surfacePressure,                         &
+      !$acc                   landIcePressureApplied,                  &
       !$acc                   potentialDensity,                        &
       !$acc                   displacedDensity,                        &
       !$acc                   boundaryLayerDepth,                      &
@@ -1433,6 +1438,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
               thermExpCoeff,                           &
               sfcFlxAttCoeff,                          &
               surfacePressure,                         &
+              landIcePressureApplied,                  &
               potentialDensity,                        &
               displacedDensity,                        &
               boundaryLayerDepth,                      &
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
index 39ce178e0ab6..75ceb4b57536 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
@@ -34,6 +34,7 @@ module ocn_wetting_drying
    use ocn_mesh
    use ocn_subgrid
    use ocn_vel_pressure_grad
+   use ocn_surface_land_ice_fluxes
 
    implicit none
    private
@@ -67,6 +68,121 @@ module ocn_wetting_drying
 
 contains
 
+   subroutine ocn_compute_land_ice_ssh_floatation(landIcePressure, ssh_floatation, &
+                                                  err)
+      !-----------------------------------------------------------------
+      !
+      ! input variables
+      !
+      !-----------------------------------------------------------------
+      real (kind=RKIND), dimension(:), intent(in) :: landIcePressure
+
+      !-----------------------------------------------------------------
+      !
+      ! input/output variables
+      !
+      !-----------------------------------------------------------------
+      real (kind=RKIND), dimension(:), intent(inout) :: ssh_floatation
+
+      !-----------------------------------------------------------------
+      !
+      ! output variables
+      !
+      !-----------------------------------------------------------------
+
+      integer, intent(out) :: err !< Output: error flag
+
+      !-----------------------------------------------------------------
+      !
+      ! local variables
+      !
+      !-----------------------------------------------------------------
+      integer :: iCell
+
+      err = 0
+
+      if ( .not. landIcePressureOn ) return
+
+      !$omp parallel
+      !$omp do schedule(runtime)
+      do iCell = 1, nCellsAll
+
+         ! the depth at which floatation occurs with full land ice pressure
+         ! (may be below the bathymetry) (gravity is positive)
+         ssh_floatation(iCell) = - landIcePressure(iCell) / (gravity * rho_sw)
+
+      end do
+      !$omp end do
+      !$omp end parallel
+
+   end subroutine ocn_compute_land_ice_ssh_floatation
+
+   subroutine ocn_compute_land_ice_ssh_adjustment(landIcePressure, landIceMask, &
+                                                  ssh_land_ice_adjustment, err)
+      !-----------------------------------------------------------------
+      !
+      ! input variables
+      !
+      !-----------------------------------------------------------------
+      integer, dimension(:), intent(in) :: landIceMask
+      real (kind=RKIND), dimension(:), intent(in) :: landIcePressure
+
+      !-----------------------------------------------------------------
+      !
+      ! input/output variables
+      !
+      !-----------------------------------------------------------------
+      real (kind=RKIND), dimension(:), intent(inout) :: ssh_land_ice_adjustment
+
+      !-----------------------------------------------------------------
+      !
+      ! output variables
+      !
+      !-----------------------------------------------------------------
+
+      integer, intent(out) :: err !< Output: error flag
+
+      !-----------------------------------------------------------------
+      !
+      ! local variables
+      !
+      !-----------------------------------------------------------------
+      integer :: iCell
+      real (kind=RKIND), dimension(:), allocatable :: landIcePressureMax, rho_mean, ssh_floatation
+      real (kind=RKIND) :: dPressure, ssh_crit, hCrit
+
+      err = 0
+
+      if ( .not. config_zero_drying_velocity ) return
+      if ( .not. landIcePressureOn ) return
+
+      allocate(rho_mean(nCellsAll))
+      allocate(landIcePressureMax(nCellsAll))
+      allocate(ssh_floatation(nCellsAll))
+      rho_mean(:) = rho_sw
+      landIcePressureMax(:) = 0.0_RKIND
+      ssh_floatation(:) = 0.0_RKIND
+
+      call ocn_compute_land_ice_ssh_floatation(landIcePressure, ssh_floatation, err)
+
+      !$omp parallel
+      !$omp do schedule(runtime) private(dPressure, ssh_crit, hCrit)
+      do iCell = 1, nCellsAll
+         hCrit = config_drying_min_cell_height * (maxLevelCell(iCell) - minLevelCell(iCell) + 1)
+         ssh_crit = hCrit - bottomDepth(iCell)
+
+         ! the ssh equivalent by which we adjusted landIcePressure to ensure
+         ! floatation at roughly the point at which the wetting-and-drying
+         ! is active (ssh_crit)
+         ! we need to multiply by landIceMask here for the cases where ssh_crit > 0
+         ssh_land_ice_adjustment(iCell) = landIceMask(iCell) * &
+                                          min(ssh_floatation(iCell) - ssh_crit, 0.0_RKIND)
+      end do
+      !$omp end do
+      !$omp end parallel
+
+   end subroutine ocn_compute_land_ice_ssh_adjustment
+
 !***********************************************************************
 !
 !  routine ocn_wetting_drying_verify
@@ -237,26 +353,32 @@ subroutine ocn_prevent_drying_rk4(domain, block, dt, rkSubstepWeight, config_zer
      !-----------------------------------------------------------------
 
 
-      type (mpas_pool_type), pointer :: tendPool
-      type (mpas_pool_type), pointer :: statePool
-      type (mpas_pool_type), pointer :: provisStatePool
-      real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur
-      real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis
-      real (kind=RKIND), dimension(:, :), pointer :: normalVelocity
+     type (mpas_pool_type), pointer :: tendPool
+     type (mpas_pool_type), pointer :: statePool
+     type (mpas_pool_type), pointer :: provisStatePool
+     type (mpas_pool_type), pointer :: forcingPool
+     real (kind=RKIND), dimension(:), pointer :: landIcePressure, landIceDraft
+     integer, dimension(:), pointer :: landIceMask
+     real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur
+     real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis
+     real (kind=RKIND), dimension(:, :), pointer :: normalVelocity
 
-      integer :: iEdge, k
+     integer :: iEdge, k
 
-      err = 0
+     err = 0
 
      call mpas_pool_get_subpool(block % structs, 'tend', tendPool)
      call mpas_pool_get_subpool(block % structs, 'state', statePool)
      call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool)
+     call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)
 
      call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 1)
      ! use thickness at n because constraint is h_n + dt*T_h > h_min
      call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1)
      call mpas_pool_get_array(provisStatePool, 'layerThickness', layerThicknessProvis, 1)
-
+     call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask)
+     call mpas_pool_get_array(forcingPool, 'landIceDraft', landIceDraft)
+     call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure)
 
      !$omp parallel
      !$omp do schedule(runtime)
@@ -269,8 +391,10 @@ subroutine ocn_prevent_drying_rk4(domain, block, dt, rkSubstepWeight, config_zer
      ! ensure cells stay wet by selectively damping cells with a damping tendency to make
      ! sure tendency doesn't dry cells
 
-     call ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
-                                             normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err)
+     call ocn_wetting_drying_wettingVelocity( &
+             domain, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
+             normalTransportVelocity, rkSubstepWeight, &
+             landIceMask, landIcePressure, landIceDraft, wettingVelocityFactor, err)
 
      ! prevent drying from happening with selective wettingVelocityFactor
      if (config_zero_drying_velocity) then
@@ -309,7 +433,8 @@ end subroutine ocn_prevent_drying_rk4 !}}}
 !-----------------------------------------------------------------------
 
    subroutine ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
-                                                 normalVelocity, dt, wettingVelocityFactor, err)!{{{
+                                                 normalVelocity, dt, landIceMask, landIcePressure, landIceDraft, &
+                                                 wettingVelocityFactor, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -332,6 +457,11 @@ subroutine ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerT
       real (kind=RKIND), intent(in) :: &
          dt     !< Input: time step
 
+      integer, dimension(:), intent(in) :: landIceMask
+
+      real (kind=RKIND), dimension(:), intent(in) :: &
+         landIcePressure, landIceDraft 
+
       !-----------------------------------------------------------------
       !
       ! input/output variables
@@ -367,10 +497,14 @@ subroutine ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerT
 
       character (len=100) :: log_string
 
-      integer:: cellDummy(2), cellCur, CellNei
+      integer:: cellDummy(2), cellCur, CellNei, nVertLevelsCell
       real (kind=RKIND) :: sshCur, sshNei
       real (kind=RKIND), dimension(:), pointer :: sshCell
       real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCell
+      real (kind=RKIND), dimension(:), allocatable :: sshLandIceAdjustment
+                              !< amount by which to adjust ssh related to
+                              !< the initial adjustment of land ice pressure
+                              !< effects on ssh
 
 
       err = 0
@@ -409,6 +543,22 @@ subroutine ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerT
           end do
 
           layerThicknessCell(k, iCell) = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell)) + dt * divOutFlux
+        end do
+      end do
+      !$omp end do
+      !$omp end parallel
+
+      allocate(sshLandIceAdjustment(nCellsAll))
+      sshLandIceAdjustment(:) = 0.0_RKIND
+      call ocn_compute_land_ice_ssh_adjustment(landIcePressure, landIceMask, sshLandIceAdjustment, err)
+
+      !$omp parallel
+      !$omp do schedule(runtime) private(k, nVertLevelsCell)
+      do iCell = 1, nCellsAll
+        nVertLevelsCell = maxLevelCell(iCell) - minLevelCell(iCell) + 1
+        do k = minLevelCell(iCell), maxLevelCell(iCell)
+          layerThicknessCell(k, iCell) = layerThicknessCell(k, iCell) + &
+                                         sshLandIceAdjustment(iCell) / nVertLevelsCell
           if ( .not. config_use_ssh_gradient_wetting_drying ) then
              call ocn_wetting_velocity_factor_on_cell_edges( &
                       wettingVelocityFactor, layerThicknessCell(k, iCell), normalVelocity, iCell, k)
@@ -458,7 +608,7 @@ subroutine ocn_wetting_drying_wettingVelocity(domain, layerThickEdgeFlux, layerT
             do k = minLevelCell(iCell), maxLevelCell(iCell)
               columnThickness = columnThickness + layerThicknessCell(k, iCell)
             enddo
-            sshCell(iCell) = columnThickness - bottomDepth(iCell)
+            sshCell(iCell) = columnThickness - bottomDepth(iCell) - sshLandIceAdjustment(iCell)
           enddo
           !$omp end do
           !$omp end parallel