Distribute unablatedVolumeDynCell among neighboring cells

If a cell is left with unablated volume after step 2c in li_apply_front_ablation_velocity,
distribute that volume among its dynamic neighbors. Testing with von Mises calving on the MISMIP+
domain indicates that this greatly reduces inaccuracies in the calving routine, but that very large
calving events can still trigger the >10% calving inaccuracy error. This could potentially
be alleviated by adding a do-while loop over the new section 3 in li_apply_front_ablation_velocity,
but we may want to control for extremely large calving events anyway.

components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F

@@ -2059,7 +2059,7 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
       real(kind=RKIND), dimension(6) :: localInfo, globalInfo
       real(kind=RKIND) :: edgeLengthScaling
       real(kind=RKIND), parameter :: ablationSmallThk = 1.0e-8_RKIND ! in meters, a small thickness threshold
-      integer :: err_tmp
+      integer :: err_tmp, iter, nDynNeighbors
 
       err = 0
       err_tmp = 0
@@ -2493,8 +2493,6 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
             if (iCell <= nCellsSolve) ablationSubtotal2 = ablationSubtotal2 + removeVolumeHere
          endif
       enddo
-      !call mpas_log_write("Done calculating ablation for dynamic floating cells. Removed $r m^3", realArgs=(/calvingSubtotal2/))
-
 
       ! Clean up to account for roundoff level errors that can occur
       do iCell = 1, nCells
@@ -2503,7 +2501,66 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
          endif
       enddo
 
-      ! Clean up - zap any stranded ice.  This only needs to be considered for floating ice.
+      ! 3. Distribute unablatedVolumeDynCell among neighboring cells. This is
+      ! necesssary when a cell is fully depleted but still has not ablated as
+      ! much ice as required by the ablation parameterization. This still may
+      ! result in large error for very high ablation rates, which could be
+      ! alleviated by adding a do-while loop.
+
+      ! a. Calculate volume to be distributed evenly among neighboring dynamic
+      ! cells
+      do iCell = 1, nCells
+         if ( (unablatedVolumeDynCell(iCell) > 0.0_RKIND) ) then
+            nDynNeighbors = 0
+            ! Count neighbors between which to evenly distribute
+            ! unablatedVolumeDynCell
+            do iEdge = 1, nEdgesOnCell(iCell)
+               iNeighbor = cellsOnCell(iEdge, iCell)
+               if ((li_mask_is_dynamic_ice(cellMask(iNeighbor))) .and. (bedTopography(iNeighbor) <= config_sea_level) &
+                    .and. (cellVolume(iNeighbor) > 0.0_RKIND) ) then
+                  nDynNeighbors = nDynNeighbors + 1
+               endif
+            enddo
+            if ( nDynNeighbors > 0 ) then
+               ! Now distribute unablatedVolumeDynCell between those neighbors
+               do iEdge = 1, nEdgesOnCell(iCell)
+                  iNeighbor = cellsOnCell(iEdge, iCell)
+                  if ((li_mask_is_dynamic_ice(cellMask(iNeighbor))) .and. (bedTopography(iNeighbor) <= config_sea_level) &
+                    .and. (cellVolume(iNeighbor) > 0.0_RKIND) ) then
+                     unablatedVolumeDynCell(iNeighbor) = unablatedVolumeDynCell(iNeighbor) + &
+                                                             unablatedVolumeDynCell(iCell) / nDynNeighbors
+                  endif
+               enddo
+               ! This cell now has no more unablatedVolume
+               unablatedVolumeDynCell(iCell) = 0.0_RKIND
+            endif
+         endif
+      enddo
+
+      ! b. Apply ablation
+      do iCell = 1, nCells
+          if ((li_mask_is_dynamic_ice(cellMask(iCell))) .and. (bedTopography(iCell) <= config_sea_level) &
+               .and. (cellVolume(iCell) > 0.0_RKIND) )then
+               removeVolumeHere = min(cellVolume(iCell), unablatedVolumeDynCell(iCell))  ! Don't use more than available
+               ablationThickness(iCell) = ablationThickness(iCell) + removeVolumeHere / areaCell(iCell) ! add to ablationThickness calculated in 2c
+               ablatedVolumeDynCell(iCell) = ablatedVolumeDynCell(iCell) + removeVolumeHere
+               unablatedVolumeDynCell(iCell) = unablatedVolumeDynCell(iCell) - removeVolumeHere
+               cellVolume(iCell) = cellVolume(iCell) - removeVolumeHere
+
+               if (iCell <= nCellsSolve) ablationSubtotal2 = ablationSubtotal2 + removeVolumeHere
+           endif
+      enddo
+      !call mpas_log_write("Done calculating ablation for dynamic floating cells. Removed $r m^3", realArgs=(/calvingSubtotal2/))
+
+      ! 4. Clean up after applying ablation velocity from parameterization
+      ! a. Clean up to account for roundoff level errors that can occur
+      do iCell = 1, nCells
+         if (abs(ablationThickness(iCell) - thickness(iCell)) < ablationSmallThk) then
+            ablationThickness(iCell) = thickness(iCell)
+         endif
+      enddo
+
+      ! b. Zap any stranded ice.  This only needs to be considered for floating ice.
       ablationSubtotal3 = 0.0_RKIND
       do iCell = 1, nCells
          if (li_mask_is_floating_ice(cellMask(iCell)) .and. (.not. li_mask_is_dynamic_ice(cellMask(iCell)))) then
@@ -2526,7 +2583,7 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
          endif
       enddo
 
-      ! Clean up to account for roundoff level errors that can occur
+      ! c. Clean up to account for roundoff level errors that can occur
       do iCell = 1, nCells
          if (abs(ablationThickness(iCell) - thickness(iCell)) < ablationSmallThk) then
             ablationThickness(iCell) = thickness(iCell)