Do not update cell, edge, or vertex masks during RK loop

Remove calls to li_calculate_mask where possible within the RK loop. Where we need masks for budget calculations, reset masks to their pre-advection states once the necessary calculation is complete.
MALI-Dev · Apr 12, 2024 · b811b20 · b811b20
1 parent 3aa29f5
commit b811b20
Show file tree

Hide file tree

Showing 4 changed files with 42 additions and 10 deletions.
diff --git 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"
 	       />
+	       <var name="edgeMaskTemporary" type="integer" dimensions="nEdges" units="none"
+               description="temporary copy of edgeMask"
+               persistence="scratch"
+	       />
+	       <var name="vertexMaskTemporary" type="integer" dimensions="nVertices" units="none"
+               description="temporary copy of vertexMask"
+               persistence="scratch"
+	       />
         </var_struct>
 
 <!-- ================ -->

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
@@ -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
@@ -370,8 +380,14 @@ subroutine li_advection_thickness_tracers(&
       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 +556,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 +643,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,6 +686,11 @@ 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
@@ -727,6 +751,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

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
@@ -3827,9 +3827,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
@@ -399,9 +399,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