From b92225df59e6f56e3616143d95760eaef840cffe Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 11 Apr 2022 11:20:27 -0600 Subject: [PATCH] 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. --- .../src/mode_forward/mpas_li_calving.F | 67 +++++++++++++++++-- 1 file changed, 62 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 00088eaa3316..f9f3a12206c5 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/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)