From a975404e388db810ee7a2576a32279b93347ac69 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 11 Apr 2024 19:08:04 -0700 Subject: [PATCH 1/3] 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 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..ccb5a1a43537 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 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..30142138d0fa 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 @@ -401,9 +401,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 From f135bdf0d0c1d5ebd0cfa1d18e975c6ee82febef Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 11 Apr 2024 19:28:52 -0700 Subject: [PATCH 2/3] Remove cellMask from vertical_remap Remove cellMask from vertical_remap, as it is not used in that routine. --- .../src/mode_forward/mpas_li_advection.F | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) 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 ccb5a1a43537..9e983c4b45aa 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 @@ -696,7 +696,7 @@ subroutine li_advection_thickness_tracers(& ! 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 @@ -1965,7 +1965,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) !----------------------------------------------------------------- ! @@ -1979,9 +1979,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 From f6a6653f31570f777f2552c4da5552f15580fb60 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 19 Aug 2024 15:43:06 -0700 Subject: [PATCH 3/3] Update masks before RK loop, but not at start of each advection stage. Update masks before RK loop, but not at start of each advection stage. --- .../src/mode_forward/mpas_li_advection.F | 4 ---- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 10 +++++++--- 2 files changed, 7 insertions(+), 7 deletions(-) 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 9e983c4b45aa..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 @@ -376,10 +376,6 @@ 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 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, 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 30142138d0fa..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', &