From b811b206ee659271479ee0749699bea9e25cd70e Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 11 Apr 2024 19:08:04 -0700 Subject: [PATCH] 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. --- .../mpas-albany-landice/src/Registry.xml | 8 +++++ .../src/mode_forward/mpas_li_advection.F | 34 ++++++++++++++++--- .../src/mode_forward/mpas_li_calving.F | 5 ++- .../mpas_li_time_integration_fe_rk.F | 5 ++- 4 files changed, 42 insertions(+), 10 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 3aa73768f04f..ebf8eefbef16 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 f06c8d1a67e4..81a5b09868c9 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 @@ -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 index de2a523e9339..81f25549da9b 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 @@ -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 index 52e3dd228446..98298f90814d 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 @@ -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