diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index e217bbba2abf..06b625353423 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -1433,6 +1433,14 @@ is the value of that variable from the *previous* time level!
description="temporary copy of cellMask"
persistence="scratch"
/>
+
+
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
index ad3f8f517c50..4cb0f3e132db 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
@@ -193,6 +193,8 @@ subroutine li_advection_thickness_tracers(&
type (field1DInteger), pointer :: &
cellMaskTemporaryField, & ! scratch field containing old values of cellMask
+ edgeMaskTemporaryField, &
+ vertexMaskTemporaryField, &
thermalCellMaskField
! Allocatable arrays need for flux-corrected transport advection
@@ -201,6 +203,7 @@ subroutine li_advection_thickness_tracers(&
integer, dimension(:), pointer :: &
cellMask, & ! integer bitmask for cells
edgeMask, & ! integer bitmask for edges
+ vertexMask, &
thermalCellMask
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -290,6 +293,7 @@ subroutine li_advection_thickness_tracers(&
call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
+ call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1)
call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening)
call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)
@@ -356,6 +360,12 @@ subroutine li_advection_thickness_tracers(&
call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField)
call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.)
+ call mpas_pool_get_field(geometryPool, 'edgeMaskTemporary', edgeMaskTemporaryField)
+ call mpas_allocate_scratch_field(edgeMaskTemporaryField, .true.)
+
+ call mpas_pool_get_field(geometryPool, 'vertexMaskTemporary', vertexMaskTemporaryField)
+ call mpas_allocate_scratch_field(vertexMaskTemporaryField, .true.)
+
call mpas_pool_get_field(scratchPool, 'iceCellMask', thermalCellMaskField)
call mpas_allocate_scratch_field(thermalCellMaskField, .true.)
thermalCellMask => thermalCellMaskField % array
@@ -366,12 +376,14 @@ subroutine li_advection_thickness_tracers(&
! given the old thickness, compute the thickness in each layer
call li_calculate_layerThickness(meshPool, thickness, layerThickness)
- ! update masks
- call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
- err = ior(err, err_tmp)
-
- ! save old copycellMask for determining cells changing from grounded to floating and vice versa
+ ! Save copies of masks because we need to preserve mask
+ ! states prior to advection for accurate time integration.
+ ! A mask update is necessary to calculate grounding line flux,
+ ! after which we will reset the masks to their previous states.
cellMaskTemporaryField % array(:) = cellMask(:)
+ edgeMaskTemporaryField % array(:) = edgeMask(:)
+ vertexMaskTemporaryField % array(:) = vertexMask(:)
+
layerThicknessEdgeFlux(:,:) = 0.0_RKIND
!-----------------------------------------------------------------
@@ -540,10 +552,10 @@ subroutine li_advection_thickness_tracers(&
dynamicThickening = (sum(layerThickness, 1) - thickness) / dt * scyr ! units of m/yr
- ! Update the thickness and cellMask before applying the mass balance.
- ! The update is needed because the SMB and BMB depend on whether ice is present.
+ ! Update the thickness before applying the mass balance, but
+ ! do not update masks because mass balance acts on geometry
+ ! before advection took place.
thickness = sum(layerThickness, 1)
- call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
@@ -627,6 +639,9 @@ subroutine li_advection_thickness_tracers(&
enddo
endif
+ ! We need an updated set of masks to calculate fluxAcrossGroundingLine,
+ ! but we will reset this to the previous state below for accuracy of the
+ ! time integration scheme.
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
err = ior(err, err_tmp)
@@ -667,12 +682,17 @@ subroutine li_advection_thickness_tracers(&
endif
enddo ! edges
+ ! Reset masks to state before advection and mass balance for
+ ! accuracy of time integration scheme.
+ cellMask(:) = cellMaskTemporaryField % array(:)
+ edgeMask(:) = edgeMaskTemporaryField % array(:)
+ vertexMask(:) = vertexMaskTemporaryField % array(:)
! Remap tracers to the standard vertical sigma coordinate
! Note: If tracers are not being advected, then this subroutine simply restores the
! layer thickness to sigma coordinate values.
- call vertical_remap(thickness, cellMask, meshPool, layerThickness, advectedTracers, err_tmp)
+ call vertical_remap(thickness, meshPool, layerThickness, advectedTracers, err_tmp)
err = ior(err, err_tmp)
if (config_print_thickness_advection_info) then
@@ -727,6 +747,8 @@ subroutine li_advection_thickness_tracers(&
call mpas_deallocate_scratch_field(basalTracersField, .true.)
call mpas_deallocate_scratch_field(surfaceTracersField, .true.)
call mpas_deallocate_scratch_field(cellMaskTemporaryField, .true.)
+ call mpas_deallocate_scratch_field(edgeMaskTemporaryField, .true.)
+ call mpas_deallocate_scratch_field(vertexMaskTemporaryField, .true.)
call mpas_deallocate_scratch_field(thermalCellMaskField, .true.)
! === error check
@@ -1939,7 +1961,7 @@ end subroutine li_layer_normal_velocity
!> OpenMP over either blocks or cells.
!
!-----------------------------------------------------------------------
- subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers, err)
+ subroutine vertical_remap(thickness, meshPool, layerThickness, tracers, err)
!-----------------------------------------------------------------
!
@@ -1953,9 +1975,6 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers
real(kind=RKIND), dimension(:), intent(in) :: &
thickness !< Input: ice thickness
- integer, dimension(:), intent(in) :: &
- cellMask !< Input: mask for cells (needed for determining presence/absence of ice)
-
!-----------------------------------------------------------------
!
! input/output variables
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 a460c3224a0e..1637c7d697df 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
@@ -3828,9 +3828,8 @@ subroutine li_finalize_damage_after_advection(domain, err)
call mpas_pool_get_array(geometryPool, 'damageNye', damageNye)
call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
- ! make sure masks are up to date. May not be necessary, but safer to do anyway.
- call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
- err = ior(err, err_tmp)
+ ! Note: In order to preserve accuracy of time integration,
+ ! we do not update masks before finalizing damage.
if (config_preserve_damage) then
do iCell = 1, nCells
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
index b7af7d342100..3f695a0b0f32 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
@@ -250,8 +250,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
err = ior(err, err_tmp)
call mpas_timer_stop("face melting for grounded ice")
-! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? ***
-
! === Basal melting for floating ice ===========
call mpas_timer_start("basal melting for floating ice")
call li_basal_melt_floating_ice(domain, err_tmp)
@@ -364,7 +362,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
// ' is not supported with config_rk_order = $i', &
intArgs=(/config_rk_order/), messageType=MPAS_LOG_ERR)
return
- endif
+ endif
+
+ ! Calculate masks prior to RK loop, but do not update masks within the loop
+ ! to preserve the accuracy of time integration.
+ call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
+ err = ior(err, err_tmp)
+
! *** Start RK loop ***
do rkStage = 1, nRKstages
call mpas_log_write('beginning rk stage $i of $i', &
@@ -401,9 +405,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + &
(1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:)
thickness = sum(layerThickness, 1)
- ! Calculate masks after updating thickness
- call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
- err = ior(err, err_tmp)
+ ! Do not calculate masks after updating thickness! We need to keep masks
+ ! constant for now to preserve accuracy of time integration
if (trim(config_thermal_solver) .ne. 'none') then
do iCell = 1, nCells