From 9dc07048b7b035f970c70482aeb27977f73cdab4 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 29 Aug 2025 12:36:26 -0700 Subject: [PATCH 01/19] Reset masks in RK loop before velocity solve Also remove mask update from velocity solver and write to logs when updating and resetting masks. --- .../mpas_li_time_integration_fe_rk.F | 21 +++++++++++++++---- .../src/mode_forward/mpas_li_velocity.F | 1 - .../src/shared/mpas_li_mask.F | 2 +- 3 files changed, 18 insertions(+), 6 deletions(-) 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 9581e32adc2d..b8c9c3c39bd5 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 @@ -146,7 +146,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) fluxAcrossGroundingLineOnCellsAccum, & calvingThicknessAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions - integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection + integer, dimension(:), allocatable :: cellMaskPrev, edgeMaskPrev, vertexMaskPrev ! cell mask before advection real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & temperaturePrev, & @@ -156,10 +156,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damage3d, & damage3dPrev integer, pointer :: nVertLevels - integer, pointer :: nCells, nEdges + integer, pointer :: nCells, nEdges, nVertices integer, pointer :: albanyVelocityError integer :: iCell1, iCell2, iEdge, theGroundedCell - integer, dimension(:), pointer :: edgeMask, cellMask + integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask real (kind=RKIND), pointer :: deltat, config_ice_density real (kind=RKIND) :: deltatFull real (kind=RKIND), dimension(4) :: rkSubstepWeights @@ -196,6 +196,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) + call mpas_pool_get_dimension(meshPool, 'nVertices', nVertices) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) call mpas_pool_get_array(velocityPool, 'albanyVelocityError', albanyVelocityError) @@ -220,6 +221,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) allocate(passiveTracer3d(nVertLevels, nCells+1)) allocate(cellMaskPrev(nCells+1)) + allocate(edgeMaskPrev(nEdges+1)) + allocate(vertexMaskPrev(nVertices+1)) allocate(sfcMassBalAccum(nCells+1)) allocate(basalMassBalAccum(nCells+1)) @@ -238,6 +241,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damage3d(:,:) = 0.0_RKIND passiveTracer3dPrev(:,:) = 0.0_RKIND cellMaskPrev(:) = 0 + edgeMaskPrev(:) = 0 + vertexMaskPrev(:) = 0 passiveTracer3d(:,:) = 0.0_RKIND sfcMassBalAccum(:) = 0.0_RKIND @@ -293,6 +298,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) call mpas_pool_get_array(geometryPool, 'damage', damage) 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(thermalPool, 'temperature', temperature) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) @@ -301,6 +309,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) waterFracPrev(:,:) = waterFrac(:,:) layerThicknessPrev(:,:) = layerThickness(:,:) cellMaskPrev(:) = cellMask(:) + edgeMaskPrev(:) = edgeMask(:) + vertexMaskPrev(:) = vertexMask(:) do k = 1, nVertLevels damage3dPrev(k,:) = damage(:) passiveTracer3dPrev(k,:) = passiveTracer2d(:) @@ -532,7 +542,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! This frequently results in icebergs that causes intermediate ! RK stage velocity solves to fail. call li_remove_icebergs(domain) + call mpas_log_write('Resetting masks') + cellMask(:) = cellMaskPrev(:) + edgeMask(:) = edgeMaskPrev(:) + vertexMask(:) = vertexMaskPrev(:) call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) err = ior(err, err_tmp) endif @@ -1281,7 +1295,6 @@ subroutine advection_solver(domain, err) call li_calculate_layerThickness(meshPool, thickness, layerThickness) endif - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) call li_update_geometry(geometryPool) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F index 4d07dbbfdb84..f75779c25def 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F @@ -333,7 +333,6 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) call li_update_geometry(geometryPool) block => block % next diff --git a/components/mpas-albany-landice/src/shared/mpas_li_mask.F b/components/mpas-albany-landice/src/shared/mpas_li_mask.F index 64050cb83adf..6f97e5ed874f 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -354,7 +354,7 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) ! ==== ! Calculate cellMask values=========================== ! ==== - + call mpas_log_write("Updating masks") !call mpas_timer_start('calculate mask cell') ! Set mask to 0 everywhere, but need to preserve bits the initial ice extent bit do i=1, nCells From dd75c61a5ff7c843e0d9927931a4ad9eeb7d3171 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 4 Mar 2026 12:52:44 -0800 Subject: [PATCH 02/19] Add missing deallocate calls Deallocate edgeMaskPrev and vertexMaskPrev where cellMaskPrev is deallocated. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 2 ++ 1 file changed, 2 insertions(+) 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 b8c9c3c39bd5..5b94d0b95679 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 @@ -663,6 +663,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) deallocate(cellMaskPrev) + deallocate(edgeMaskPrev) + deallocate(vertexMaskPrev) deallocate(sfcMassBalAccum) deallocate(basalMassBalAccum) deallocate(groundedSfcMassBalAccum) From b8739a77ba0248bdd1eb74115f43380e7d814eb6 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 4 Mar 2026 13:00:36 -0800 Subject: [PATCH 03/19] Remove outdated comment Remove outdated comment about updating masks before velocity solve. --- .../mpas-albany-landice/src/mode_forward/mpas_li_velocity.F | 1 - 1 file changed, 1 deletion(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F index f75779c25def..20a049cefc43 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F @@ -326,7 +326,6 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_dmpar_field_halo_exch(domain, 'thickness') call mpas_timer_stop("halo updates") - ! Update mask just to be safe (might be redundant) block => domain % blocklist do while (associated(block)) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) From 2bdeb340015efd5f0d560c4d5ada880e3bc3771f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 4 Mar 2026 13:11:16 -0800 Subject: [PATCH 04/19] Add another missing deallocation call Add call to deallocate fluxAcrossGroundingLineOnCellsAccum --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 1 + 1 file changed, 1 insertion(+) 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 5b94d0b95679..14525173a22c 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 @@ -671,6 +671,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(groundedBasalMassBalAccum) deallocate(floatingBasalMassBalAccum) deallocate(fluxAcrossGroundingLineAccum) + deallocate(fluxAcrossGroundingLineOnCellsAccum) deallocate(calvingThicknessAccum) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta From aff72f22e055a767227eb0f6bfadf735a7b44aed Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 4 Mar 2026 19:34:32 -0800 Subject: [PATCH 05/19] Control mask updates before velocity solve Control mask updates before velocity solve with a logical flag in the call to li_velocity_solve --- .../mpas-albany-landice/src/mode_forward/mpas_li_core.F | 2 +- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 4 ++-- .../src/mode_forward/mpas_li_velocity.F | 9 +++++++-- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index 1d11ad87bfda..32d4ab083ea7 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -819,7 +819,7 @@ function li_core_initial_solve(domain) result(err) call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool) call mpas_pool_get_array(velocityPool, 'albanyVelocityError', albanyVelocityError) albanyVelocityError = 0 - call li_velocity_solve(domain, solveVelo=solveVelo, err=err_tmp) + call li_velocity_solve(domain, solveVelo=solveVelo, updateMask=.true., err=err_tmp) err = ior(err, err_tmp) ! Initial calving solution 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 14525173a22c..6510164eb317 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 @@ -547,7 +547,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) cellMask(:) = cellMaskPrev(:) edgeMask(:) = edgeMaskPrev(:) vertexMask(:) = vertexMaskPrev(:) - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + call li_velocity_solve(domain, solveVelo=.true., updateMask=.false., err=err_tmp) err = ior(err, err_tmp) endif ! *** end RK loop *** @@ -637,7 +637,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! === Solve velocity for final state ===================== if (solveVeloAfterCalving) then - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + call li_velocity_solve(domain, solveVelo=.true., updateMask=.true., err=err_tmp) err = ior(err, err_tmp) endif diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F index 20a049cefc43..6b01e5124ba4 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F @@ -222,7 +222,7 @@ end subroutine li_velocity_block_init !> This routine calls velocity solvers. ! !----------------------------------------------------------------------- - subroutine li_velocity_solve(domain, solveVelo, err) + subroutine li_velocity_solve(domain, solveVelo, updateMask, err) use mpas_vector_reconstruction use li_mask @@ -236,6 +236,8 @@ subroutine li_velocity_solve(domain, solveVelo, err) !< velocity should be solved, or if an existing solution should be !< used to calculate the diagnostic fields related to velocty !< (e.g. on a restart) + logical, intent(in) :: updateMask !< Input: update mask before + !< velocity solve only if outside of RK loop. !----------------------------------------------------------------- ! input/output variables @@ -331,7 +333,10 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - + ! Only update mask outside of RK loop + if ( updateMask ) then + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + endif call li_update_geometry(geometryPool) block => block % next From 1da70688a366a292981102f2e8d165bda7ae94d3 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 12:58:23 -0700 Subject: [PATCH 06/19] Add RK pairs for embedded SSP RK time integration Add RK pairs for embedded SSP RK time integration. These pairs were suggested by Co-pilot and are not necesarily aligned with pairs suggested in the literature. --- .../mpas-albany-landice/src/Registry.xml | 32 + .../mpas_li_time_integration_fe_rk.F | 573 ++++++++++++++---- 2 files changed, 483 insertions(+), 122 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index a14205e0d783..0295c2ac9279 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -558,6 +558,38 @@ description="Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2." possible_values="3 or 4" /> + + + + + + + + 1) .and. (rkStage < 4) ) then - call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) - call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) - call mpas_pool_get_array(geometryPool, 'damage', damage) + includeThermalForError = trim(config_thermal_solver) .ne. 'none' + maxStepAttempts = 1 + if (useEmbeddedSSPRK) then + maxStepAttempts = max(1, config_embedded_rk_max_retries + 1) + endif - call mpas_pool_get_array(thermalPool, 'temperature', temperature) - call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + stepAccepted = .false. + do stepAttempt = 1, maxStepAttempts + ! Reset prognostic state to pre-RK values before each attempt. + layerThickness(:,:) = layerThicknessPrev(:,:) + thickness(:) = sum(layerThickness, 1) + temperature(:,:) = temperaturePrev(:,:) + waterFrac(:,:) = waterFracPrev(:,:) + damage(:) = damagePrev(:) + passiveTracer2d(:) = passiveTracer2dPrev(:) + cellMask(:) = cellMaskPrev(:) + edgeMask(:) = edgeMaskPrev(:) + vertexMask(:) = vertexMaskPrev(:) + + sfcMassBalAccum(:) = 0.0_RKIND + basalMassBalAccum(:) = 0.0_RKIND + groundedSfcMassBalAccum(:) = 0.0_RKIND + groundedBasalMassBalAccum(:) = 0.0_RKIND + floatingBasalMassBalAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineOnCellsAccum(:) = 0.0_RKIND + calvingThicknessAccum(:) = 0.0_RKIND + calvingThickness(:) = 0.0_RKIND + calvingThicknessFromThreshold(:) = 0.0_RKIND + calvingVelocity(:) = 0.0_RKIND + + if (useEmbeddedSSPRK) then + layerThicknessEmbedded(:,:) = layerThicknessPrev(:,:) + temperatureEmbedded(:,:) = temperaturePrev(:,:) + waterFracEmbedded(:,:) = waterFracPrev(:,:) + damageEmbedded(:) = damagePrev(:) + passiveTracer2dEmbedded(:) = passiveTracer2dPrev(:) + endif - layerThicknessTmp(:,:) = layerThickness(:,:) - layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & - (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) - thickness = sum(layerThickness, 1) - ! Do not calculate masks after updating thickness! We need to keep masks - ! constant for now to preserve accuracy of time integration +! *** Start RK loop *** + do rkStage = 1, nRKstages + call mpas_log_write('beginning rk stage $i of $i', & + intArgs=(/rkStage, nRKstages/)) + deltat = deltatFull * rkSubstepWeights(rkStage) - if (trim(config_thermal_solver) .ne. 'none') then - do iCell = 1, nCells - do k = 1, nVertLevels - if (layerThickness(k,iCell) > 0.0_RKIND) then - temperature(k,iCell) = ( rkSSPweights(rkStage) * temperaturePrev(k,iCell) * layerThicknessPrev(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * temperature(k,iCell) * & - layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) - waterFrac(k,iCell) = ( rkSSPweights(rkStage) * waterFracPrev(k,iCell) * layerThicknessPrev(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * waterFrac(k,iCell) * & - layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) - endif - enddo - enddo + ! === calculate damage =========== + if (config_calculate_damage) then + call mpas_timer_start("damage") + call li_calculate_damage(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("damage") endif - if (config_calculate_damage) then + ! === Compute new state for prognostic variables === + call mpas_timer_start("advect thickness and tracers") + call advection_solver(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("advect thickness and tracers") + + ! If using SSP RK, then update thickness and tracers incrementally. + ! For first RK stage, thickness and tracer updates above are sufficient. + ! Likewise, for the 4-stage SSP RK3 the last stage is just a forward euler update. + if ( (rkStage > 1) .and. (rkStage < 4) ) then + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) + call mpas_pool_get_array(geometryPool, 'damage', damage) + + call mpas_pool_get_array(thermalPool, 'temperature', temperature) + call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + + layerThicknessTmp(:,:) = layerThickness(:,:) + layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & + (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) + thickness = sum(layerThickness, 1) + ! 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 + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + temperature(k,iCell) = ( rkSSPweights(rkStage) * temperaturePrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * temperature(k,iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + waterFrac(k,iCell) = ( rkSSPweights(rkStage) * waterFracPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * waterFrac(k,iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + endif + enddo + enddo + endif + + if (config_calculate_damage) then + do iCell = 1, nCells + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + else + damage3d(k,iCell) = 0.0_RKIND + endif + enddo + damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) + enddo + endif + do iCell = 1, nCells do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then - damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * & - layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) else - damage3d(k,iCell) = 0.0_RKIND + passiveTracer3d(k,iCell) = 0.0_RKIND endif enddo - damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) + passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo + + endif + + ! === Ensure damage is within bounds before velocity solve === + if ( config_finalize_damage_after_advection ) then + call mpas_timer_start("finalize damage") + call li_finalize_damage_after_advection(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("finalize damage") endif - do iCell = 1, nCells - do k = 1, nVertLevels - if (layerThickness(k,iCell) > 0.0_RKIND) then - passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * & - layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) - else - passiveTracer3d(k,iCell) = 0.0_RKIND + call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + + ! update budgets + sfcMassBalAccum = sfcMassBalAccum + rkTendWeights(rkStage) * sfcMassBalApplied + groundedSfcMassBalAccum = groundedSfcMassBalAccum + rkTendWeights(rkStage) * groundedSfcMassBalApplied + basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied + groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied + floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied + fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine + fluxAcrossGroundingLineOnCellsAccum = fluxAcrossGroundingLineOnCellsAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLineOnCells + + ! Halo updates + call mpas_timer_start("halo updates") + + call mpas_dmpar_field_halo_exch(domain, 'thickness') + call mpas_dmpar_field_halo_exch(domain, 'temperature') + call mpas_dmpar_field_halo_exch(domain, 'waterFrac') + call mpas_dmpar_field_halo_exch(domain, 'damage') + call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') + + call mpas_timer_stop("halo updates") + + if (useEmbeddedSSPRK) then + if (rkStage == 1) then + layerThicknessStage1(:,:) = layerThickness(:,:) + if (includeThermalForError) then + temperatureStage1(:,:) = temperature(:,:) + waterFracStage1(:,:) = waterFrac(:,:) endif - enddo - passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) - enddo + damageStage1(:) = damage(:) + passiveTracer2dStage1(:) = passiveTracer2d(:) + elseif (rkStage == 2) then + layerThicknessStage2(:,:) = layerThickness(:,:) + if (includeThermalForError) then + temperatureStage2(:,:) = temperature(:,:) + waterFracStage2(:,:) = waterFrac(:,:) + endif + damageStage2(:) = damage(:) + passiveTracer2dStage2(:) = passiveTracer2d(:) + endif + endif - endif - - ! === Ensure damage is within bounds before velocity solve === - if ( config_finalize_damage_after_advection ) then - call mpas_timer_start("finalize damage") - call li_finalize_damage_after_advection(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("finalize damage") + ! Update velocity for each RK step + ! === Solve Velocity ===================== + if ( config_update_velocity_before_calving .or. ( (.not. config_update_velocity_before_calving) & + .and. (rkStage < nRKstages) ) ) then + + if (config_restore_calving_front) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + calvingThicknessAccum = calvingThicknessAccum + rkTendWeights(rkStage) * calvingThickness + endif + ! We need to remove icebergs between RK stages because the + ! main calving routine is not called until after the RK loop. + ! This frequently results in icebergs that causes intermediate + ! RK stage velocity solves to fail. + call li_remove_icebergs(domain) + call mpas_log_write('Resetting masks') + + cellMask(:) = cellMaskPrev(:) + edgeMask(:) = edgeMaskPrev(:) + vertexMask(:) = vertexMaskPrev(:) + call li_velocity_solve(domain, solveVelo=.true., updateMask=.false., err=err_tmp) + err = ior(err, err_tmp) + endif +! *** end RK loop *** + enddo + + if (.not. useEmbeddedSSPRK) then + stepAccepted = .true. + exit endif - call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) - call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) - call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) - call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) - call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) - call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) - call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) - call mpas_pool_get_array(geometryPool, 'thickness', thickness) - - ! update budgets - sfcMassBalAccum = sfcMassBalAccum + rkTendWeights(rkStage) * sfcMassBalApplied - groundedSfcMassBalAccum = groundedSfcMassBalAccum + rkTendWeights(rkStage) * groundedSfcMassBalApplied - basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied - groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied - floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied - fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine - fluxAcrossGroundingLineOnCellsAccum = fluxAcrossGroundingLineOnCellsAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLineOnCells - - ! Halo updates - call mpas_timer_start("halo updates") + if (config_rk_order == 2) then + ! SSPRK2(2,2) with embedded RK1 uses stage-1 Forward Euler endpoint. + layerThicknessEmbedded(:,:) = layerThicknessStage1(:,:) + if (includeThermalForError) then + temperatureEmbedded(:,:) = temperatureStage1(:,:) + waterFracEmbedded(:,:) = waterFracStage1(:,:) + endif + damageEmbedded(:) = damageStage1(:) + passiveTracer2dEmbedded(:) = passiveTracer2dStage1(:) + elseif (config_rk3_stages == 3) then + ! SSPRK3(3,3) embedded RK2: y* = y0 + h/2*(k1+k2) = 2*y_stage2 - y0. + layerThicknessEmbedded(:,:) = 2.0_RKIND * layerThicknessStage2(:,:) - layerThicknessPrev(:,:) + if (includeThermalForError) then + temperatureEmbedded(:,:) = 2.0_RKIND * temperatureStage2(:,:) - temperaturePrev(:,:) + waterFracEmbedded(:,:) = 2.0_RKIND * waterFracStage2(:,:) - waterFracPrev(:,:) + endif + damageEmbedded(:) = 2.0_RKIND * damageStage2(:) - damagePrev(:) + passiveTracer2dEmbedded(:) = 2.0_RKIND * passiveTracer2dStage2(:) - passiveTracer2dPrev(:) + else + ! SSPRK3(4,3) embedded RK2 midpoint: y* = y0 + h*k2 = y0 + 2*(y_stage2 - y_stage1). + layerThicknessEmbedded(:,:) = layerThicknessPrev(:,:) + & + 2.0_RKIND * (layerThicknessStage2(:,:) - layerThicknessStage1(:,:)) + if (includeThermalForError) then + temperatureEmbedded(:,:) = temperaturePrev(:,:) + & + 2.0_RKIND * (temperatureStage2(:,:) - temperatureStage1(:,:)) + waterFracEmbedded(:,:) = waterFracPrev(:,:) + & + 2.0_RKIND * (waterFracStage2(:,:) - waterFracStage1(:,:)) + endif + damageEmbedded(:) = damagePrev(:) + 2.0_RKIND * (damageStage2(:) - damageStage1(:)) + passiveTracer2dEmbedded(:) = passiveTracer2dPrev(:) + & + 2.0_RKIND * (passiveTracer2dStage2(:) - passiveTracer2dStage1(:)) + endif - call mpas_dmpar_field_halo_exch(domain, 'thickness') - call mpas_dmpar_field_halo_exch(domain, 'temperature') - call mpas_dmpar_field_halo_exch(domain, 'waterFrac') - call mpas_dmpar_field_halo_exch(domain, 'damage') - call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') + call li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & + waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & + damage, damageEmbedded, nVertLevels, nCells, includeThermalForError, & + config_calculate_damage, config_embedded_rk_rtol, & + config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, rkErrorNorm) + + if (rkErrorNorm <= 1.0_RKIND) then + stepAccepted = .true. + call mpas_log_write('Embedded SSPRK accepted attempt $i with normalized error $r', & + intArgs=(/stepAttempt/), realArgs=(/rkErrorNorm/)) + exit + endif - call mpas_timer_stop("halo updates") + if (stepAttempt == maxStepAttempts) then + err = 1 + call mpas_log_write('Embedded SSPRK failed to satisfy tolerance after $i attempts (error=$r)', & + intArgs=(/maxStepAttempts/), realArgs=(/rkErrorNorm/), & + messageType=MPAS_LOG_ERR) + return + endif - ! Update velocity for each RK step - ! === Solve Velocity ===================== - if ( config_update_velocity_before_calving .or. ( (.not. config_update_velocity_before_calving) & - .and. (rkStage < nRKstages) ) ) then + dtScale = config_embedded_rk_safety_factor * rkErrorNorm ** (-1.0_RKIND / real(embeddedHighOrder, RKIND)) + dtScale = min(config_embedded_rk_max_growth, dtScale) + dtScale = max(config_embedded_rk_min_shrink, dtScale) + deltatFull = deltatFull * dtScale - if (config_restore_calving_front) then - ! restore the calving front to its initial position before velocity solve. - call li_restore_calving_front(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) - calvingThicknessAccum = calvingThicknessAccum + rkTendWeights(rkStage) * calvingThickness - endif - ! We need to remove icebergs between RK stages because the - ! main calving routine is not called until after the RK loop. - ! This frequently results in icebergs that causes intermediate - ! RK stage velocity solves to fail. - call li_remove_icebergs(domain) - call mpas_log_write('Resetting masks') - - cellMask(:) = cellMaskPrev(:) - edgeMask(:) = edgeMaskPrev(:) - vertexMask(:) = vertexMaskPrev(:) - call li_velocity_solve(domain, solveVelo=.true., updateMask=.false., err=err_tmp) - err = ior(err, err_tmp) + if (deltatFull < config_min_adaptive_timestep) then + err = 1 + call mpas_log_write('Embedded SSPRK requested deltat below config_min_adaptive_timestep', & + messageType=MPAS_LOG_ERR) + return endif -! *** end RK loop *** + + call mpas_log_write('Embedded SSPRK rejected attempt $i; reducing deltat by factor $r to $r seconds', & + intArgs=(/stepAttempt/), realArgs=(/dtScale, deltatFull/)) enddo + if (.not. stepAccepted) then + err = 1 + call mpas_log_write('Embedded SSPRK did not accept a step', messageType=MPAS_LOG_ERR) + return + endif + ! Finalize budget updates ! Update masks after RK integration call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) @@ -658,13 +898,30 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(waterFracPrev) deallocate(layerThicknessPrev) deallocate(layerThicknessTmp) + deallocate(layerThicknessEmbedded) + deallocate(layerThicknessStage1) + deallocate(layerThicknessStage2) deallocate(damage3dPrev) deallocate(damage3d) deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) + deallocate(temperatureEmbedded) + deallocate(waterFracEmbedded) + deallocate(temperatureStage1) + deallocate(temperatureStage2) + deallocate(waterFracStage1) + deallocate(waterFracStage2) deallocate(cellMaskPrev) deallocate(edgeMaskPrev) deallocate(vertexMaskPrev) + deallocate(damagePrev) + deallocate(passiveTracer2dPrev) + deallocate(damageEmbedded) + deallocate(passiveTracer2dEmbedded) + deallocate(damageStage1) + deallocate(damageStage2) + deallocate(passiveTracer2dStage1) + deallocate(passiveTracer2dStage2) deallocate(sfcMassBalAccum) deallocate(basalMassBalAccum) deallocate(groundedSfcMassBalAccum) @@ -1573,6 +1830,78 @@ subroutine calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, er !-------------------------------------------------------------------- end subroutine calculate_layerThicknessEdge +!*********************************************************************** +! +! routine li_embedded_ssprk_error_norm +! +!> \brief Compute normalized error estimate for embedded SSPRK control +!> \details +!> Uses the current RK endpoint and a formally reconstructed embedded +!> lower-order endpoint to compute a max-norm over thickness and +!> tracer-like fields. +!----------------------------------------------------------------------- + subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & + waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & + damage, damageEmbedded, nVertLevels, nCells, includeThermal, includeDamage, & + relTol, absTolThickness, absTolTracer, errorNorm) + + real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness, layerThicknessEmbedded + real (kind=RKIND), dimension(:,:), intent(in) :: temperature, temperatureEmbedded + real (kind=RKIND), dimension(:,:), intent(in) :: waterFrac, waterFracEmbedded + real (kind=RKIND), dimension(:), intent(in) :: passiveTracer2d, passiveTracer2dEmbedded + real (kind=RKIND), dimension(:), intent(in) :: damage, damageEmbedded + integer, intent(in) :: nVertLevels, nCells + logical, intent(in) :: includeThermal, includeDamage + real (kind=RKIND), intent(in) :: relTol, absTolThickness, absTolTracer + real (kind=RKIND), intent(out) :: errorNorm + + integer :: iCell, k + real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr + + errorNorm = 0.0_RKIND + + do iCell = 1, nCells + yPrimary = sum(layerThickness(:, iCell)) + yEmbedded = sum(layerThicknessEmbedded(:, iCell)) + scaleVal = absTolThickness + relTol * max(abs(yPrimary), abs(yEmbedded)) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(yPrimary - yEmbedded) / scaleVal + errorNorm = max(errorNorm, thisErr) + endif + + scaleVal = absTolTracer + relTol * max(abs(passiveTracer2d(iCell)), abs(passiveTracer2dEmbedded(iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(passiveTracer2d(iCell) - passiveTracer2dEmbedded(iCell)) / scaleVal + errorNorm = max(errorNorm, thisErr) + endif + + if (includeDamage) then + scaleVal = absTolTracer + relTol * max(abs(damage(iCell)), abs(damageEmbedded(iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(damage(iCell) - damageEmbedded(iCell)) / scaleVal + errorNorm = max(errorNorm, thisErr) + endif + endif + + if (includeThermal) then + do k = 1, nVertLevels + scaleVal = absTolTracer + relTol * max(abs(temperature(k,iCell)), abs(temperatureEmbedded(k,iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(temperature(k,iCell) - temperatureEmbedded(k,iCell)) / scaleVal + errorNorm = max(errorNorm, thisErr) + endif + + scaleVal = absTolTracer + relTol * max(abs(waterFrac(k,iCell)), abs(waterFracEmbedded(k,iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(waterFrac(k,iCell) - waterFracEmbedded(k,iCell)) / scaleVal + errorNorm = max(errorNorm, thisErr) + endif + enddo + endif + enddo + + end subroutine li_embedded_ssprk_error_norm + !*********************************************************************** ! ! routine advance_clock From 11a114134acd91222e804a8b9320339e8b43e814 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 13:04:41 -0700 Subject: [PATCH 07/19] Only allocate arrays when embedded pairs are enabled Only allocate extra arrays used in embedded methods when embedded RK is actually being used. This will save memory overhead. --- .../mpas_li_time_integration_fe_rk.F | 96 ++++++++++--------- 1 file changed, 51 insertions(+), 45 deletions(-) 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 6c0886513cb4..603d9e7772f9 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 @@ -262,30 +262,33 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(waterFracPrev(nVertLevels, nCells+1)) allocate(layerThicknessPrev(nVertLevels, nCells+1)) allocate(layerThicknessTmp(nVertLevels, nCells+1)) - allocate(layerThicknessEmbedded(nVertLevels, nCells+1)) - allocate(layerThicknessStage1(nVertLevels, nCells+1)) - allocate(layerThicknessStage2(nVertLevels, nCells+1)) allocate(damage3dPrev(nVertLevels, nCells+1)) allocate(damage3d(nVertLevels, nCells+1)) allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) allocate(passiveTracer3d(nVertLevels, nCells+1)) - allocate(temperatureEmbedded(nVertLevels, nCells+1)) - allocate(waterFracEmbedded(nVertLevels, nCells+1)) - allocate(temperatureStage1(nVertLevels, nCells+1)) - allocate(temperatureStage2(nVertLevels, nCells+1)) - allocate(waterFracStage1(nVertLevels, nCells+1)) - allocate(waterFracStage2(nVertLevels, nCells+1)) allocate(cellMaskPrev(nCells+1)) allocate(edgeMaskPrev(nEdges+1)) allocate(vertexMaskPrev(nVertices+1)) allocate(damagePrev(nCells+1)) allocate(passiveTracer2dPrev(nCells+1)) - allocate(damageEmbedded(nCells+1)) - allocate(passiveTracer2dEmbedded(nCells+1)) - allocate(damageStage1(nCells+1)) - allocate(damageStage2(nCells+1)) - allocate(passiveTracer2dStage1(nCells+1)) - allocate(passiveTracer2dStage2(nCells+1)) + + if (config_use_embedded_ssprk) then + allocate(layerThicknessEmbedded(nVertLevels, nCells+1)) + allocate(layerThicknessStage1(nVertLevels, nCells+1)) + allocate(layerThicknessStage2(nVertLevels, nCells+1)) + allocate(temperatureEmbedded(nVertLevels, nCells+1)) + allocate(waterFracEmbedded(nVertLevels, nCells+1)) + allocate(temperatureStage1(nVertLevels, nCells+1)) + allocate(temperatureStage2(nVertLevels, nCells+1)) + allocate(waterFracStage1(nVertLevels, nCells+1)) + allocate(waterFracStage2(nVertLevels, nCells+1)) + allocate(damageEmbedded(nCells+1)) + allocate(passiveTracer2dEmbedded(nCells+1)) + allocate(damageStage1(nCells+1)) + allocate(damageStage2(nCells+1)) + allocate(passiveTracer2dStage1(nCells+1)) + allocate(passiveTracer2dStage2(nCells+1)) + endif allocate(sfcMassBalAccum(nCells+1)) allocate(basalMassBalAccum(nCells+1)) @@ -300,9 +303,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) waterFracPrev(:,:) = 0.0_RKIND layerThicknessPrev(:,:) = 0.0_RKIND layerThicknessTmp(:,:) = 0.0_RKIND - layerThicknessEmbedded(:,:) = 0.0_RKIND - layerThicknessStage1(:,:) = 0.0_RKIND - layerThicknessStage2(:,:) = 0.0_RKIND damage3dPrev(:,:) = 0.0_RKIND damage3d(:,:) = 0.0_RKIND passiveTracer3dPrev(:,:) = 0.0_RKIND @@ -310,20 +310,26 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) edgeMaskPrev(:) = 0 vertexMaskPrev(:) = 0 passiveTracer3d(:,:) = 0.0_RKIND - temperatureEmbedded(:,:) = 0.0_RKIND - waterFracEmbedded(:,:) = 0.0_RKIND - temperatureStage1(:,:) = 0.0_RKIND - temperatureStage2(:,:) = 0.0_RKIND - waterFracStage1(:,:) = 0.0_RKIND - waterFracStage2(:,:) = 0.0_RKIND damagePrev(:) = 0.0_RKIND passiveTracer2dPrev(:) = 0.0_RKIND - damageEmbedded(:) = 0.0_RKIND - passiveTracer2dEmbedded(:) = 0.0_RKIND - damageStage1(:) = 0.0_RKIND - damageStage2(:) = 0.0_RKIND - passiveTracer2dStage1(:) = 0.0_RKIND - passiveTracer2dStage2(:) = 0.0_RKIND + + if (config_use_embedded_ssprk) then + layerThicknessEmbedded(:,:) = 0.0_RKIND + layerThicknessStage1(:,:) = 0.0_RKIND + layerThicknessStage2(:,:) = 0.0_RKIND + temperatureEmbedded(:,:) = 0.0_RKIND + waterFracEmbedded(:,:) = 0.0_RKIND + temperatureStage1(:,:) = 0.0_RKIND + temperatureStage2(:,:) = 0.0_RKIND + waterFracStage1(:,:) = 0.0_RKIND + waterFracStage2(:,:) = 0.0_RKIND + damageEmbedded(:) = 0.0_RKIND + passiveTracer2dEmbedded(:) = 0.0_RKIND + damageStage1(:) = 0.0_RKIND + damageStage2(:) = 0.0_RKIND + passiveTracer2dStage1(:) = 0.0_RKIND + passiveTracer2dStage2(:) = 0.0_RKIND + endif sfcMassBalAccum(:) = 0.0_RKIND basalMassBalAccum(:) = 0.0_RKIND @@ -898,30 +904,30 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(waterFracPrev) deallocate(layerThicknessPrev) deallocate(layerThicknessTmp) - deallocate(layerThicknessEmbedded) - deallocate(layerThicknessStage1) - deallocate(layerThicknessStage2) deallocate(damage3dPrev) deallocate(damage3d) deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) - deallocate(temperatureEmbedded) - deallocate(waterFracEmbedded) - deallocate(temperatureStage1) - deallocate(temperatureStage2) - deallocate(waterFracStage1) - deallocate(waterFracStage2) deallocate(cellMaskPrev) deallocate(edgeMaskPrev) deallocate(vertexMaskPrev) deallocate(damagePrev) deallocate(passiveTracer2dPrev) - deallocate(damageEmbedded) - deallocate(passiveTracer2dEmbedded) - deallocate(damageStage1) - deallocate(damageStage2) - deallocate(passiveTracer2dStage1) - deallocate(passiveTracer2dStage2) + if (allocated(layerThicknessEmbedded)) deallocate(layerThicknessEmbedded) + if (allocated(layerThicknessStage1)) deallocate(layerThicknessStage1) + if (allocated(layerThicknessStage2)) deallocate(layerThicknessStage2) + if (allocated(temperatureEmbedded)) deallocate(temperatureEmbedded) + if (allocated(waterFracEmbedded)) deallocate(waterFracEmbedded) + if (allocated(temperatureStage1)) deallocate(temperatureStage1) + if (allocated(temperatureStage2)) deallocate(temperatureStage2) + if (allocated(waterFracStage1)) deallocate(waterFracStage1) + if (allocated(waterFracStage2)) deallocate(waterFracStage2) + if (allocated(damageEmbedded)) deallocate(damageEmbedded) + if (allocated(passiveTracer2dEmbedded)) deallocate(passiveTracer2dEmbedded) + if (allocated(damageStage1)) deallocate(damageStage1) + if (allocated(damageStage2)) deallocate(damageStage2) + if (allocated(passiveTracer2dStage1)) deallocate(passiveTracer2dStage1) + if (allocated(passiveTracer2dStage2)) deallocate(passiveTracer2dStage2) deallocate(sfcMassBalAccum) deallocate(basalMassBalAccum) deallocate(groundedSfcMassBalAccum) From 0faa297feb46db72e5411d46c0535b754f9cfbcf Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 13:30:07 -0700 Subject: [PATCH 08/19] Synchronize embedded norm across ranks Synchronize embedded norm across ranks so we don't get a crash when different ranks have different norms, and thus each set a different dt. --- .../mpas_li_time_integration_fe_rk.F | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) 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 603d9e7772f9..0c1c3f96effd 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 @@ -110,6 +110,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) type (block_type), pointer :: block integer :: err_tmp type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool, velocityPool + type (mpas_pool_type), pointer :: meshPoolTmp logical, pointer :: config_restore_calving_front, & config_restore_calving_front_prevent_retreat @@ -187,7 +188,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) integer :: iCell1, iCell2, iEdge, theGroundedCell integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask real (kind=RKIND), pointer :: deltat, config_ice_density + real (kind=RKIND), pointer :: deltatTmp real (kind=RKIND) :: deltatFull + real (kind=RKIND) :: deltatStage real (kind=RKIND), dimension(4) :: rkSubstepWeights real (kind=RKIND), dimension(4) :: rkSSPweights real (kind=RKIND), dimension(4) :: rkTendWeights ! Weights used for calculating budget terms @@ -195,7 +198,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) logical :: solveVeloAfterCalving logical :: useEmbeddedSSPRK, includeThermalForError, stepAccepted integer :: stepAttempt, maxStepAttempts, embeddedHighOrder - real (kind=RKIND) :: rkErrorNorm, dtScale + real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, dtScale err = 0 err_tmp = 0 @@ -207,6 +210,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) maxStepAttempts = 1 embeddedHighOrder = 3 rkErrorNorm = 0.0_RKIND + rkErrorNormGlobal = 0.0_RKIND dtScale = 1.0_RKIND call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) @@ -553,7 +557,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) do rkStage = 1, nRKstages call mpas_log_write('beginning rk stage $i of $i', & intArgs=(/rkStage, nRKstages/)) - deltat = deltatFull * rkSubstepWeights(rkStage) + deltatStage = deltatFull * rkSubstepWeights(rkStage) + ! Keep deltat consistent on all blocks for retry attempts. + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPoolTmp) + call mpas_pool_get_array(meshPoolTmp, 'deltat', deltatTmp) + deltatTmp = deltatStage + block => block % next + enddo + deltat = deltatStage ! === calculate damage =========== if (config_calculate_damage) then @@ -762,6 +775,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) config_calculate_damage, config_embedded_rk_rtol, & config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, rkErrorNorm) + call mpas_dmpar_max_real(domain % dminfo, rkErrorNorm, rkErrorNormGlobal) + rkErrorNorm = rkErrorNormGlobal + if (rkErrorNorm <= 1.0_RKIND) then stepAccepted = .true. call mpas_log_write('Embedded SSPRK accepted attempt $i with normalized error $r', & From d1fc6449dd6465e0c57d477b670410aefc17ec8f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 13:43:41 -0700 Subject: [PATCH 09/19] Do not use damage to calculate error. Ignor damage and (nearly) ice-free regions for error norm calculation. --- .../mode_forward/mpas_li_time_integration_fe_rk.F | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) 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 0c1c3f96effd..e7eeedb1a364 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 @@ -196,7 +196,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) real (kind=RKIND), dimension(4) :: rkTendWeights ! Weights used for calculating budget terms integer :: nRKstages logical :: solveVeloAfterCalving - logical :: useEmbeddedSSPRK, includeThermalForError, stepAccepted + logical :: useEmbeddedSSPRK, includeThermalForError, includeDamageForError, stepAccepted integer :: stepAttempt, maxStepAttempts, embeddedHighOrder real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, dtScale @@ -205,6 +205,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) solveVeloAfterCalving = .false. useEmbeddedSSPRK = .false. includeThermalForError = .false. + includeDamageForError = .false. stepAccepted = .false. stepAttempt = 1 maxStepAttempts = 1 @@ -515,6 +516,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) embeddedHighOrder = config_rk_order includeThermalForError = trim(config_thermal_solver) .ne. 'none' + includeDamageForError = .false. + if (useEmbeddedSSPRK .and. config_calculate_damage) then + call mpas_log_write('Embedded SSPRK error control excludes damage field to avoid non-smooth limiter dominance.') + endif maxStepAttempts = 1 if (useEmbeddedSSPRK) then maxStepAttempts = max(1, config_embedded_rk_max_retries + 1) @@ -772,7 +777,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & damage, damageEmbedded, nVertLevels, nCells, includeThermalForError, & - config_calculate_damage, config_embedded_rk_rtol, & + includeDamageForError, config_embedded_rk_rtol, & config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, rkErrorNorm) call mpas_dmpar_max_real(domain % dminfo, rkErrorNorm, rkErrorNormGlobal) @@ -1879,6 +1884,7 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, integer :: iCell, k real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr + real (kind=RKIND), parameter :: minIceForError = 1.0e-12_RKIND errorNorm = 0.0_RKIND @@ -1891,6 +1897,10 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, errorNorm = max(errorNorm, thisErr) endif + if (max(yPrimary, yEmbedded) <= minIceForError) then + cycle + endif + scaleVal = absTolTracer + relTol * max(abs(passiveTracer2d(iCell)), abs(passiveTracer2dEmbedded(iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(passiveTracer2d(iCell) - passiveTracer2dEmbedded(iCell)) / scaleVal @@ -1907,6 +1917,7 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, if (includeThermal) then do k = 1, nVertLevels + if (max(layerThickness(k,iCell), layerThicknessEmbedded(k,iCell)) <= minIceForError) cycle scaleVal = absTolTracer + relTol * max(abs(temperature(k,iCell)), abs(temperatureEmbedded(k,iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(temperature(k,iCell) - temperatureEmbedded(k,iCell)) / scaleVal From 6b83c9ef91a998181bd941d4b6e5cbda8ae9d915 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 14:00:24 -0700 Subject: [PATCH 10/19] Make clock advancement consistent with new deltat Make clock advancement consistent with new deltat when deltat is set by the embedded RK error norm. This means moving the clock update to a different point in the time step, which may cause changes relative to previous behavior. --- .../mpas_li_time_integration_fe_rk.F | 57 +++++++++++-------- 1 file changed, 32 insertions(+), 25 deletions(-) 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 e7eeedb1a364..f852a3c932db 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 @@ -199,6 +199,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) logical :: useEmbeddedSSPRK, includeThermalForError, includeDamageForError, stepAccepted integer :: stepAttempt, maxStepAttempts, embeddedHighOrder real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, dtScale + real (kind=RKIND), save :: previousEmbeddedAcceptedDt = -1.0_RKIND err = 0 err_tmp = 0 @@ -352,12 +353,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advection prep") -! === Advance the clock before all other physics happen =========== - call mpas_timer_start("advancing clock") - call advance_clock(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("advancing clock") - !TODO: Determine whether grounded melting should in fact be called first ! === Ocean forcing extrapolation into ice-shelf cavities =========== call mpas_timer_start("ocean forcing extrapolation") @@ -520,6 +515,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) if (useEmbeddedSSPRK .and. config_calculate_damage) then call mpas_log_write('Embedded SSPRK error control excludes damage field to avoid non-smooth limiter dominance.') endif + + if (useEmbeddedSSPRK .and. (previousEmbeddedAcceptedDt > 0.0_RKIND)) then + deltatFull = min(deltatFull, previousEmbeddedAcceptedDt * config_embedded_rk_max_growth) + endif maxStepAttempts = 1 if (useEmbeddedSSPRK) then maxStepAttempts = max(1, config_embedded_rk_max_retries + 1) @@ -820,6 +819,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) return endif + if (useEmbeddedSSPRK) previousEmbeddedAcceptedDt = deltatFull + + ! Advance model time only once, after this step has been accepted. + call mpas_timer_start("advancing clock") + call advance_clock(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("advancing clock") + ! Finalize budget updates ! Update masks after RK integration call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) @@ -858,7 +865,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) endif ! Reset time step to full length after RK loop - deltat = deltatFull + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPoolTmp) + call mpas_pool_get_array(meshPoolTmp, 'deltat', deltatTmp) + deltatTmp = deltatFull + block => block % next + enddo + deltat = deltatFull ! === Update subglacial hydrology =========== ! It's not clear where the best place to call this should be. @@ -1977,8 +1991,6 @@ subroutine advance_clock(domain, err) type (MPAS_TimeInterval_type) :: timeStepInterval !< the current time step as an interval type (MPAS_Time_type) :: simulationStartTime_timeType - logical, pointer :: config_adaptive_timestep - character (len=StrKIND), pointer :: xtime character (len=StrKIND), pointer :: simulationStartTime character(len=StrKIND) :: timeStamp !< current time as a string @@ -1990,24 +2002,19 @@ subroutine advance_clock(domain, err) err = 0 - call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep) - ! === - ! === Adaptive timestep: update clock information + ! === Update clock timestep from mesh deltat ! === - ! Set time step in clock object since the time step could have changed - ! Need to get value out of a block, but all blocks should have the same value, so just use first block - if (config_adaptive_timestep) then - block => domain % blocklist - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_array(meshPool, 'deltat', deltat) - ! convert dt in seconds to timeInterval type - call mpas_set_timeInterval(timeStepInterval, dt=deltat, ierr=err_tmp) - err = ior(err,err_tmp) - ! update the clock with the timeInterval - call mpas_set_clock_timestep(domain % clock, timeStepInterval, err_tmp) - err = ior(err,err_tmp) - endif + ! Need to get value out of a block, but all blocks should have the same value, so just use first block. + block => domain % blocklist + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array(meshPool, 'deltat', deltat) + ! convert dt in seconds to timeInterval type + call mpas_set_timeInterval(timeStepInterval, dt=deltat, ierr=err_tmp) + err = ior(err,err_tmp) + ! update the clock with the timeInterval + call mpas_set_clock_timestep(domain % clock, timeStepInterval, err_tmp) + err = ior(err,err_tmp) ! === ! === Update clock information From ecae4a4e8ceda5b93a78a99dec1c7ebf6bbcab8d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 14:24:03 -0700 Subject: [PATCH 11/19] Fix timekeeping issue in set_timestep Fix issue in set_timestep that was causing crash when using dt determined by embedded RK. --- .../mpas_li_time_integration_fe_rk.F | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) 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 f852a3c932db..3ad4dc648a93 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 @@ -1038,6 +1038,7 @@ subroutine prepare_advection(domain, err) logical, pointer :: config_print_thickness_advection_info logical, pointer :: config_adaptive_timestep logical, pointer :: config_adaptive_timestep_include_DCFL + logical, pointer :: config_use_embedded_ssprk character (len=StrKIND), pointer :: & config_thickness_advection, & ! method for advecting thickness and tracers @@ -1084,6 +1085,7 @@ subroutine prepare_advection(domain, err) real (kind=RKIND), pointer :: deltat ! variable in blocks real (kind=RKIND) :: dtSeconds ! local variable + real (kind=RKIND) :: embeddedRKDt integer :: err_tmp @@ -1092,6 +1094,7 @@ subroutine prepare_advection(domain, err) call mpas_pool_get_config(liConfigs, 'config_print_thickness_advection_info', config_print_thickness_advection_info) call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep) call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep_include_DCFL', config_adaptive_timestep_include_DCFL) + call mpas_pool_get_config(liConfigs, 'config_use_embedded_ssprk', config_use_embedded_ssprk) call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration) call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order) @@ -1262,10 +1265,19 @@ subroutine prepare_advection(domain, err) ! Set adaptive timestep if needed if (config_adaptive_timestep) then + block => domain % blocklist + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array(meshPool, 'deltat', deltat) + embeddedRKDt = deltat + call mpas_pool_get_array(geometryPool, 'calvingCFLdt', calvingCFLdt) call mpas_pool_get_array(geometryPool, 'faceMeltingCFLdt', faceMeltingCFLdt) call mpas_pool_get_array(meshPool, 'processLimitingTimestep', processLimitingTimestep) - call set_timestep(allowableAdvecDtAllProcs, allowableDiffDtAllProcs, calvingCFLdt, faceMeltingCFLdt, domain % clock, dtSeconds, processLimitingTimestep, err_tmp) + if (config_use_embedded_ssprk) then + call set_timestep(allowableAdvecDtAllProcs, allowableDiffDtAllProcs, calvingCFLdt, faceMeltingCFLdt, domain % clock, dtSeconds, processLimitingTimestep, err_tmp, embeddedRKDt) + else + call set_timestep(allowableAdvecDtAllProcs, allowableDiffDtAllProcs, calvingCFLdt, faceMeltingCFLdt, domain % clock, dtSeconds, processLimitingTimestep, err_tmp) + endif err = ior(err,err_tmp) ! Set new value on all blocks block => domain % blocklist @@ -1622,7 +1634,7 @@ end subroutine advection_solver !> This routine sdjusts the time step based on the CFL condition. ! !----------------------------------------------------------------------- - subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMeltingCFLdt, clock, dtSeconds, CFLprocess, err) + subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMeltingCFLdt, clock, dtSeconds, CFLprocess, err, embeddedRKDt) use mpas_timekeeping !----------------------------------------------------------------- @@ -1632,6 +1644,7 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel real (kind=RKIND), intent(in) :: allowableDiffDt real (kind=RKIND), intent(in) :: calvingCFLdt, faceMeltingCFLdt type (MPAS_Clock_type), intent(in) :: clock + real (kind=RKIND), intent(in), optional :: embeddedRKDt !----------------------------------------------------------------- ! output variables @@ -1656,6 +1669,7 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel type (MPAS_Time_type) :: nextForceTime, currTime type (MPAS_TimeInterval_type) :: intervalToNextForceTime real (kind=RKIND) :: secondsToNextForceTime + real (kind=RKIND) :: fracSecondsToNextForceTime real (kind=RKIND) :: allowableDt real (kind=RKIND) :: proposedDt integer :: err_tmp @@ -1730,6 +1744,14 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel ! Take minimum of the max adaptive timestep setting and the allowable dt from the CFL condition proposedDt = min(allowableDt, config_max_adaptive_timestep) + ! If embedded SSPRK is active, do not let set_timestep increase dt beyond + ! the current embedded-RK-controlled timestep. + if (present(embeddedRKDt)) then + if (embeddedRKDt > 0.0_RKIND) then + proposedDt = min(proposedDt, embeddedRKDt) + endif + endif + ! Round down the proposed dt to avoid complications with fractional seconds ! (some timekeeping-related functionality, like restarts, don't support them) ! (need to perform floor with an 8-bit integer to allow up 293 billion years in seconds) From 1177905e15cc970d9e9b2dd9aa8aa75324fc0b64 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 14:46:05 -0700 Subject: [PATCH 12/19] Fix deltat update --- .../mpas_li_time_integration_fe_rk.F | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) 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 3ad4dc648a93..0ffaa6858c11 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 @@ -821,6 +821,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) if (useEmbeddedSSPRK) previousEmbeddedAcceptedDt = deltatFull + ! Restore full accepted timestep on all blocks before advancing the clock. + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPoolTmp) + call mpas_pool_get_array(meshPoolTmp, 'deltat', deltatTmp) + deltatTmp = deltatFull + block => block % next + enddo + deltat = deltatFull + ! Advance model time only once, after this step has been accepted. call mpas_timer_start("advancing clock") call advance_clock(domain, err_tmp) @@ -864,16 +874,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) advMask2ndOrder) endif -! Reset time step to full length after RK loop - block => domain % blocklist - do while (associated(block)) - call mpas_pool_get_subpool(block % structs, 'mesh', meshPoolTmp) - call mpas_pool_get_array(meshPoolTmp, 'deltat', deltatTmp) - deltatTmp = deltatFull - block => block % next - enddo - deltat = deltatFull - ! === Update subglacial hydrology =========== ! It's not clear where the best place to call this should be. ! Seems sensible to put it after thermal evolution is complete to get updated basal melting source term. From bab0fdcf01193bbb665f8615e972d9f3313c6efa Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 15:22:51 -0700 Subject: [PATCH 13/19] Floor dt before advancing clock Floord RK-determined dt before advancing clock, which is necessary to avoid throwing an error with non-integer interval to next forcing time. --- .../mpas_li_time_integration_fe_rk.F | 35 +++++++------------ 1 file changed, 12 insertions(+), 23 deletions(-) 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 0ffaa6858c11..c2a52c239860 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 @@ -819,6 +819,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) return endif + if (useEmbeddedSSPRK) then + deltatFull = real(floor(deltatFull, KIND=8), RKIND) + if (deltatFull < config_min_adaptive_timestep) then + err = 1 + call mpas_log_write('Embedded SSPRK accepted deltat falls below config_min_adaptive_timestep after flooring', & + messageType=MPAS_LOG_ERR) + return + endif + endif + if (useEmbeddedSSPRK) previousEmbeddedAcceptedDt = deltatFull ! Restore full accepted timestep on all blocks before advancing the clock. @@ -1038,7 +1048,6 @@ subroutine prepare_advection(domain, err) logical, pointer :: config_print_thickness_advection_info logical, pointer :: config_adaptive_timestep logical, pointer :: config_adaptive_timestep_include_DCFL - logical, pointer :: config_use_embedded_ssprk character (len=StrKIND), pointer :: & config_thickness_advection, & ! method for advecting thickness and tracers @@ -1085,7 +1094,6 @@ subroutine prepare_advection(domain, err) real (kind=RKIND), pointer :: deltat ! variable in blocks real (kind=RKIND) :: dtSeconds ! local variable - real (kind=RKIND) :: embeddedRKDt integer :: err_tmp @@ -1094,7 +1102,6 @@ subroutine prepare_advection(domain, err) call mpas_pool_get_config(liConfigs, 'config_print_thickness_advection_info', config_print_thickness_advection_info) call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep) call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep_include_DCFL', config_adaptive_timestep_include_DCFL) - call mpas_pool_get_config(liConfigs, 'config_use_embedded_ssprk', config_use_embedded_ssprk) call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration) call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order) @@ -1265,19 +1272,10 @@ subroutine prepare_advection(domain, err) ! Set adaptive timestep if needed if (config_adaptive_timestep) then - block => domain % blocklist - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_array(meshPool, 'deltat', deltat) - embeddedRKDt = deltat - call mpas_pool_get_array(geometryPool, 'calvingCFLdt', calvingCFLdt) call mpas_pool_get_array(geometryPool, 'faceMeltingCFLdt', faceMeltingCFLdt) call mpas_pool_get_array(meshPool, 'processLimitingTimestep', processLimitingTimestep) - if (config_use_embedded_ssprk) then - call set_timestep(allowableAdvecDtAllProcs, allowableDiffDtAllProcs, calvingCFLdt, faceMeltingCFLdt, domain % clock, dtSeconds, processLimitingTimestep, err_tmp, embeddedRKDt) - else - call set_timestep(allowableAdvecDtAllProcs, allowableDiffDtAllProcs, calvingCFLdt, faceMeltingCFLdt, domain % clock, dtSeconds, processLimitingTimestep, err_tmp) - endif + call set_timestep(allowableAdvecDtAllProcs, allowableDiffDtAllProcs, calvingCFLdt, faceMeltingCFLdt, domain % clock, dtSeconds, processLimitingTimestep, err_tmp) err = ior(err,err_tmp) ! Set new value on all blocks block => domain % blocklist @@ -1634,7 +1632,7 @@ end subroutine advection_solver !> This routine sdjusts the time step based on the CFL condition. ! !----------------------------------------------------------------------- - subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMeltingCFLdt, clock, dtSeconds, CFLprocess, err, embeddedRKDt) + subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMeltingCFLdt, clock, dtSeconds, CFLprocess, err) use mpas_timekeeping !----------------------------------------------------------------- @@ -1644,7 +1642,6 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel real (kind=RKIND), intent(in) :: allowableDiffDt real (kind=RKIND), intent(in) :: calvingCFLdt, faceMeltingCFLdt type (MPAS_Clock_type), intent(in) :: clock - real (kind=RKIND), intent(in), optional :: embeddedRKDt !----------------------------------------------------------------- ! output variables @@ -1744,14 +1741,6 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel ! Take minimum of the max adaptive timestep setting and the allowable dt from the CFL condition proposedDt = min(allowableDt, config_max_adaptive_timestep) - ! If embedded SSPRK is active, do not let set_timestep increase dt beyond - ! the current embedded-RK-controlled timestep. - if (present(embeddedRKDt)) then - if (embeddedRKDt > 0.0_RKIND) then - proposedDt = min(proposedDt, embeddedRKDt) - endif - endif - ! Round down the proposed dt to avoid complications with fractional seconds ! (some timekeeping-related functionality, like restarts, don't support them) ! (need to perform floor with an 8-bit integer to allow up 293 billion years in seconds) From b822df809908c9e645c2959f06301ace6ce9f6e3 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 15:50:43 -0700 Subject: [PATCH 14/19] Use global rms error instead of max Use global rms error instead of max, which was causing dt to go to zero. --- .../mpas_li_time_integration_fe_rk.F | 53 ++++++++++++++----- 1 file changed, 40 insertions(+), 13 deletions(-) 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 c2a52c239860..660e6e671a44 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 @@ -198,7 +198,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) logical :: solveVeloAfterCalving logical :: useEmbeddedSSPRK, includeThermalForError, includeDamageForError, stepAccepted integer :: stepAttempt, maxStepAttempts, embeddedHighOrder - real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, dtScale + integer :: rkErrorCountLocal, rkErrorCountGlobal + real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, rkErrorSumSqLocal, rkErrorSumSqGlobal, dtScale real (kind=RKIND), save :: previousEmbeddedAcceptedDt = -1.0_RKIND err = 0 @@ -213,6 +214,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) embeddedHighOrder = 3 rkErrorNorm = 0.0_RKIND rkErrorNormGlobal = 0.0_RKIND + rkErrorCountLocal = 0 + rkErrorCountGlobal = 0 + rkErrorSumSqLocal = 0.0_RKIND + rkErrorSumSqGlobal = 0.0_RKIND dtScale = 1.0_RKIND call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) @@ -777,9 +782,17 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & damage, damageEmbedded, nVertLevels, nCells, includeThermalForError, & includeDamageForError, config_embedded_rk_rtol, & - config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, rkErrorNorm) - - call mpas_dmpar_max_real(domain % dminfo, rkErrorNorm, rkErrorNormGlobal) + config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, & + rkErrorNorm, rkErrorCountLocal) + + rkErrorSumSqLocal = rkErrorNorm * rkErrorNorm * real(rkErrorCountLocal, RKIND) + call mpas_dmpar_sum_real(domain % dminfo, rkErrorSumSqLocal, rkErrorSumSqGlobal) + call mpas_dmpar_sum_int(domain % dminfo, rkErrorCountLocal, rkErrorCountGlobal) + if (rkErrorCountGlobal > 0) then + rkErrorNormGlobal = sqrt(rkErrorSumSqGlobal / real(rkErrorCountGlobal, RKIND)) + else + rkErrorNormGlobal = 0.0_RKIND + endif rkErrorNorm = rkErrorNormGlobal if (rkErrorNorm <= 1.0_RKIND) then @@ -1889,13 +1902,13 @@ end subroutine calculate_layerThicknessEdge !> \brief Compute normalized error estimate for embedded SSPRK control !> \details !> Uses the current RK endpoint and a formally reconstructed embedded -!> lower-order endpoint to compute a max-norm over thickness and -!> tracer-like fields. +!> lower-order endpoint to compute an RMSE-like normalized error over +!> thickness and tracer-like fields. !----------------------------------------------------------------------- subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & damage, damageEmbedded, nVertLevels, nCells, includeThermal, includeDamage, & - relTol, absTolThickness, absTolTracer, errorNorm) + relTol, absTolThickness, absTolTracer, errorNorm, errorCount) real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness, layerThicknessEmbedded real (kind=RKIND), dimension(:,:), intent(in) :: temperature, temperatureEmbedded @@ -1906,12 +1919,15 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, logical, intent(in) :: includeThermal, includeDamage real (kind=RKIND), intent(in) :: relTol, absTolThickness, absTolTracer real (kind=RKIND), intent(out) :: errorNorm + integer, intent(out) :: errorCount integer :: iCell, k - real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr + real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr, errorSumSq real (kind=RKIND), parameter :: minIceForError = 1.0e-12_RKIND errorNorm = 0.0_RKIND + errorCount = 0 + errorSumSq = 0.0_RKIND do iCell = 1, nCells yPrimary = sum(layerThickness(:, iCell)) @@ -1919,7 +1935,8 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, scaleVal = absTolThickness + relTol * max(abs(yPrimary), abs(yEmbedded)) if (scaleVal > 0.0_RKIND) then thisErr = abs(yPrimary - yEmbedded) / scaleVal - errorNorm = max(errorNorm, thisErr) + errorSumSq = errorSumSq + thisErr * thisErr + errorCount = errorCount + 1 endif if (max(yPrimary, yEmbedded) <= minIceForError) then @@ -1929,14 +1946,16 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, scaleVal = absTolTracer + relTol * max(abs(passiveTracer2d(iCell)), abs(passiveTracer2dEmbedded(iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(passiveTracer2d(iCell) - passiveTracer2dEmbedded(iCell)) / scaleVal - errorNorm = max(errorNorm, thisErr) + errorSumSq = errorSumSq + thisErr * thisErr + errorCount = errorCount + 1 endif if (includeDamage) then scaleVal = absTolTracer + relTol * max(abs(damage(iCell)), abs(damageEmbedded(iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(damage(iCell) - damageEmbedded(iCell)) / scaleVal - errorNorm = max(errorNorm, thisErr) + errorSumSq = errorSumSq + thisErr * thisErr + errorCount = errorCount + 1 endif endif @@ -1946,18 +1965,26 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, scaleVal = absTolTracer + relTol * max(abs(temperature(k,iCell)), abs(temperatureEmbedded(k,iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(temperature(k,iCell) - temperatureEmbedded(k,iCell)) / scaleVal - errorNorm = max(errorNorm, thisErr) + errorSumSq = errorSumSq + thisErr * thisErr + errorCount = errorCount + 1 endif scaleVal = absTolTracer + relTol * max(abs(waterFrac(k,iCell)), abs(waterFracEmbedded(k,iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(waterFrac(k,iCell) - waterFracEmbedded(k,iCell)) / scaleVal - errorNorm = max(errorNorm, thisErr) + errorSumSq = errorSumSq + thisErr * thisErr + errorCount = errorCount + 1 endif enddo endif enddo + if (errorCount > 0) then + errorNorm = sqrt(errorSumSq / real(errorCount, RKIND)) + else + errorNorm = 0.0_RKIND + endif + end subroutine li_embedded_ssprk_error_norm !*********************************************************************** From 0f341ec263686a436f18d789094a2c67762cfcc5 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 16:35:12 -0700 Subject: [PATCH 15/19] Apply face-melt speed after advection Apply face-melt speed after advection, so it uses the dt determined by the embedded RK method. --- .../src/mode_forward/mpas_li_iceshelf_melt.F | 76 +++++++++++++------ .../mpas_li_time_integration_fe_rk.F | 10 ++- 2 files changed, 60 insertions(+), 26 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F index c9294f1d2c6f..f6d2890df3af 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F @@ -78,7 +78,7 @@ module li_iceshelf_melt !----------------------------------------------------------------------- - subroutine li_face_melt_grounded_ice(domain, err) + subroutine li_face_melt_grounded_ice(domain, err, applyFaceMeltNow, computeFaceMeltNow) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- @@ -88,6 +88,8 @@ subroutine li_face_melt_grounded_ice(domain, err) !----------------------------------------------------------------- type (domain_type), intent(inout) :: & domain !< Input/Output: domain object + logical, intent(in), optional :: applyFaceMeltNow !< If present and false, compute melt fields/CFL only + logical, intent(in), optional :: computeFaceMeltNow !< If present and false, apply existing faceMeltingThickness field !----------------------------------------------------------------- ! output variables @@ -111,9 +113,15 @@ subroutine li_face_melt_grounded_ice(domain, err) integer, dimension(:), pointer :: cellMask integer :: err_tmp logical :: applyToFloating, applyToGrounded, applyToGroundingLine + logical :: doApplyFaceMelt + logical :: doComputeFaceMelt err = 0 err_tmp = 0 + doApplyFaceMelt = .true. + doComputeFaceMelt = .true. + if (present(applyFaceMeltNow)) doApplyFaceMelt = applyFaceMeltNow + if (present(computeFaceMeltNow)) doComputeFaceMelt = computeFaceMeltNow call mpas_pool_get_config(liConfigs, 'config_front_mass_bal_grounded', config_front_mass_bal_grounded) call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_float', config_basal_mass_bal_float) @@ -133,8 +141,11 @@ subroutine li_face_melt_grounded_ice(domain, err) if ( trim(config_front_mass_bal_grounded) == 'ismip6' & .or. trim(config_front_mass_bal_grounded) == 'uniform' ) then call grounded_face_melt_ismip6(domain, applyToGrounded, & - applyToFloating, applyToGroundingLine, err_tmp) + applyToFloating, applyToGroundingLine, err_tmp, & + computeFaceMeltSpeedNow=doComputeFaceMelt, applyFrontAblationNow=doApplyFaceMelt) err = ior(err, err_tmp) + + if (.not. doApplyFaceMelt) return block => domain % blocklist do while (associated(block)) @@ -1405,7 +1416,7 @@ end subroutine iceshelf_melt_ismip6 subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & - applyToFloating, applyToGroundingLine, err) + applyToFloating, applyToGroundingLine, err, computeFaceMeltSpeedNow, applyFrontAblationNow) use li_calving use li_constants, only: oceanFreezingTempDepthDependence @@ -1416,6 +1427,8 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & type (domain_type), intent(inout) :: & domain !< Input/Output: domain object logical, intent(in) :: applyToFloating, applyToGrounded, applyToGroundingLine + logical, intent(in), optional :: computeFaceMeltSpeedNow + logical, intent(in), optional :: applyFrontAblationNow !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- @@ -1464,11 +1477,16 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & logical, pointer :: config_calculate_thermal_forcing logical, pointer :: config_ocean_data_extrapolation integer, pointer :: config_ocean_data_extrap_ncells_extra + logical :: doComputeFaceMeltSpeed, doApplyFrontAblation real (kind=RKIND) :: dzAccumulated, dz integer :: err_tmp, kk, iLayer err = 0 call mpas_log_write('Starting face melt routine') + doComputeFaceMeltSpeed = .true. + doApplyFrontAblation = .true. + if (present(computeFaceMeltSpeedNow)) doComputeFaceMeltSpeed = computeFaceMeltSpeedNow + if (present(applyFrontAblationNow)) doApplyFrontAblation = applyFrontAblationNow ! Get sea level, bedTopography, ice density call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi) @@ -1511,7 +1529,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(geometryPool, 'faceMeltingCFLdt', faceMeltingCFLdt) call mpas_pool_get_array(geometryPool, 'dtFaceMeltingCFLratio', dtFaceMeltingCFLratio) - if ( config_use_3d_thermal_forcing_for_face_melt ) then + if (doComputeFaceMeltSpeed .and. config_use_3d_thermal_forcing_for_face_melt ) then call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', zOcean) @@ -1624,29 +1642,37 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & TFocean = TFocean + ismip6shelfMelt_deltaT !apply dT correction endif - faceMeltSpeed(:) = 0.0_RKIND - allocate(faceMeltSpeedVertAvg(nCells+1)) - faceMeltSpeedVertAvg(:) = 0.0_RKIND + if (doComputeFaceMeltSpeed) then + faceMeltSpeed(:) = 0.0_RKIND - ! Calculate face melting for each cell - do iCell = 1, nCells + ! Calculate face-melt speed for each cell + do iCell = 1, nCells - if ( bedTopography(iCell) < 0.0_RKIND ) then - waterDepth = seaLevel - bedTopography(iCell) - else - waterDepth = 0.0_RKIND - endif + if ( bedTopography(iCell) < 0.0_RKIND ) then + waterDepth = seaLevel - bedTopography(iCell) + else + waterDepth = 0.0_RKIND + endif - if (trim(config_front_mass_bal_grounded) == 'ismip6') then - ! Calculate ice front melt rate at each cell - faceMeltSpeed(iCell) = (aSubglacial * waterDepth * & ! m s^-1 - (ismip6Runoff(iCell) / rhow * secPerDay)**alphaSubglacial & - + B) * max(0.0_RKIND, TFocean(iCell) + addTFocean)**betaTF / secPerDay - elseif (trim(config_front_mass_bal_grounded) == 'uniform') then - faceMeltSpeed(iCell) = config_uniform_face_melt_rate - endif + if (trim(config_front_mass_bal_grounded) == 'ismip6') then + ! Calculate ice front melt rate at each cell + faceMeltSpeed(iCell) = (aSubglacial * waterDepth * & ! m s^-1 + (ismip6Runoff(iCell) / rhow * secPerDay)**alphaSubglacial & + + B) * max(0.0_RKIND, TFocean(iCell) + addTFocean)**betaTF / secPerDay + elseif (trim(config_front_mass_bal_grounded) == 'uniform') then + faceMeltSpeed(iCell) = config_uniform_face_melt_rate + endif - enddo + enddo + endif + + if (.not. doApplyFrontAblation) then + block => block % next + cycle + endif + + allocate(faceMeltSpeedVertAvg(nCells+1)) + faceMeltSpeedVertAvg(:) = 0.0_RKIND where ( (thickness > 0.0_RKIND) .and. (lowerSurface < 0.0_RKIND) ) faceMeltSpeedVertAvg(:) = faceMeltSpeed(:) * abs(lowerSurface(:)) / thickness(:) @@ -1660,12 +1686,12 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & maxDt=faceMeltingCFLdt, CFLratio=dtFaceMeltingCFLratio, err=err_tmp) err = ior(err, err_tmp) + deallocate(faceMeltSpeedVertAvg) + block => block % next enddo ! associated(block) - deallocate(faceMeltSpeedVertAvg) - end subroutine grounded_face_melt_ismip6 !----------------------------------------------------------------------- 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 660e6e671a44..2c364453fa31 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 @@ -367,7 +367,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! === Face melting for grounded ice =========== call mpas_timer_start("face melting for grounded ice") - call li_face_melt_grounded_ice(domain, err_tmp) + ! Compute face-melt tendencies/CFL diagnostics now, but defer applying + ! the thickness update until after RK accepts the timestep. + call li_face_melt_grounded_ice(domain, err_tmp, applyFaceMeltNow=.false.) err = ior(err, err_tmp) call mpas_timer_stop("face melting for grounded ice") @@ -854,6 +856,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo deltat = deltatFull + ! Apply grounded face melt once, using the accepted timestep. + call mpas_timer_start("face melting for grounded ice") + call li_face_melt_grounded_ice(domain, err_tmp, computeFaceMeltNow=.false.) + err = ior(err, err_tmp) + call mpas_timer_stop("face melting for grounded ice") + ! Advance model time only once, after this step has been accepted. call mpas_timer_start("advancing clock") call advance_clock(domain, err_tmp) From 064c1b7fdbe98559436b48c7c86d095ce588e92f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 12 Mar 2026 17:48:34 -0700 Subject: [PATCH 16/19] Use SSP for embedded pair to RK(4,3) Use a SSP RK(4,2) for the embedded pair to RK(4,3), which prevents the time step from going to 0. --- .../mpas_li_time_integration_fe_rk.F | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) 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 2c364453fa31..7f434d8bcab0 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 @@ -766,18 +766,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damageEmbedded(:) = 2.0_RKIND * damageStage2(:) - damagePrev(:) passiveTracer2dEmbedded(:) = 2.0_RKIND * passiveTracer2dStage2(:) - passiveTracer2dPrev(:) else - ! SSPRK3(4,3) embedded RK2 midpoint: y* = y0 + h*k2 = y0 + 2*(y_stage2 - y_stage1). - layerThicknessEmbedded(:,:) = layerThicknessPrev(:,:) + & - 2.0_RKIND * (layerThicknessStage2(:,:) - layerThicknessStage1(:,:)) + ! SSPRK3(4,3): use stage-2 endpoint as SSP-consistent lower-order companion. + layerThicknessEmbedded(:,:) = layerThicknessStage2(:,:) if (includeThermalForError) then - temperatureEmbedded(:,:) = temperaturePrev(:,:) + & - 2.0_RKIND * (temperatureStage2(:,:) - temperatureStage1(:,:)) - waterFracEmbedded(:,:) = waterFracPrev(:,:) + & - 2.0_RKIND * (waterFracStage2(:,:) - waterFracStage1(:,:)) + temperatureEmbedded(:,:) = temperatureStage2(:,:) + waterFracEmbedded(:,:) = waterFracStage2(:,:) endif - damageEmbedded(:) = damagePrev(:) + 2.0_RKIND * (damageStage2(:) - damageStage1(:)) - passiveTracer2dEmbedded(:) = passiveTracer2dPrev(:) + & - 2.0_RKIND * (passiveTracer2dStage2(:) - passiveTracer2dStage1(:)) + damageEmbedded(:) = damageStage2(:) + passiveTracer2dEmbedded(:) = passiveTracer2dStage2(:) endif call li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & From f581f119aad1af20c18babd424c01c40d7b8524d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 13 Mar 2026 11:35:55 -0700 Subject: [PATCH 17/19] Change embedded pairs for RK2 and RK(3,3) Use the optimized embedded pair from Fekete et al. (2022) for RK2, and the second-stage endpoint as a pair for RK(3,3). These work much better than the previous pairs; for a 4km Thwaites mesh with historical ISMIP6 forcing and 2nd/3rd order advection, they result in a small but stable time step of about 1.5 days. Meanwhile, the pair for RK(4,3) is allowing multi-month time steps, so it is unclear why these are so much more restrictive. --- .../mpas_li_time_integration_fe_rk.F | 33 ++++++++++++------- 1 file changed, 21 insertions(+), 12 deletions(-) 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 7f434d8bcab0..a54d094fc972 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 @@ -200,6 +200,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) integer :: stepAttempt, maxStepAttempts, embeddedHighOrder integer :: rkErrorCountLocal, rkErrorCountGlobal real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, rkErrorSumSqLocal, rkErrorSumSqGlobal, dtScale + real (kind=RKIND) :: rk2EmbeddedW1, rk2EmbeddedW2 real (kind=RKIND), save :: previousEmbeddedAcceptedDt = -1.0_RKIND err = 0 @@ -218,6 +219,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkErrorCountGlobal = 0 rkErrorSumSqLocal = 0.0_RKIND rkErrorSumSqGlobal = 0.0_RKIND + rk2EmbeddedW1 = 0.694021459207626_RKIND + rk2EmbeddedW2 = 0.305978540792374_RKIND dtScale = 1.0_RKIND call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) @@ -748,23 +751,29 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) endif if (config_rk_order == 2) then - ! SSPRK2(2,2) with embedded RK1 uses stage-1 Forward Euler endpoint. - layerThicknessEmbedded(:,:) = layerThicknessStage1(:,:) + ! SSPRK2(2,2) embedded companion from Fekete et al. (2022), + ! using convex stage-state weights. + layerThicknessEmbedded(:,:) = rk2EmbeddedW1 * layerThicknessStage1(:,:) + & + rk2EmbeddedW2 * layerThickness(:,:) if (includeThermalForError) then - temperatureEmbedded(:,:) = temperatureStage1(:,:) - waterFracEmbedded(:,:) = waterFracStage1(:,:) + temperatureEmbedded(:,:) = rk2EmbeddedW1 * temperatureStage1(:,:) + & + rk2EmbeddedW2 * temperature(:,:) + waterFracEmbedded(:,:) = rk2EmbeddedW1 * waterFracStage1(:,:) + & + rk2EmbeddedW2 * waterFrac(:,:) endif - damageEmbedded(:) = damageStage1(:) - passiveTracer2dEmbedded(:) = passiveTracer2dStage1(:) + damageEmbedded(:) = rk2EmbeddedW1 * damageStage1(:) + & + rk2EmbeddedW2 * damage(:) + passiveTracer2dEmbedded(:) = rk2EmbeddedW1 * passiveTracer2dStage1(:) + & + rk2EmbeddedW2 * passiveTracer2d(:) elseif (config_rk3_stages == 3) then - ! SSPRK3(3,3) embedded RK2: y* = y0 + h/2*(k1+k2) = 2*y_stage2 - y0. - layerThicknessEmbedded(:,:) = 2.0_RKIND * layerThicknessStage2(:,:) - layerThicknessPrev(:,:) + ! SSPRK3(3,3): use stage-2 endpoint as lower-order companion. + layerThicknessEmbedded(:,:) = layerThicknessStage2(:,:) if (includeThermalForError) then - temperatureEmbedded(:,:) = 2.0_RKIND * temperatureStage2(:,:) - temperaturePrev(:,:) - waterFracEmbedded(:,:) = 2.0_RKIND * waterFracStage2(:,:) - waterFracPrev(:,:) + temperatureEmbedded(:,:) = temperatureStage2(:,:) + waterFracEmbedded(:,:) = waterFracStage2(:,:) endif - damageEmbedded(:) = 2.0_RKIND * damageStage2(:) - damagePrev(:) - passiveTracer2dEmbedded(:) = 2.0_RKIND * passiveTracer2dStage2(:) - passiveTracer2dPrev(:) + damageEmbedded(:) = damageStage2(:) + passiveTracer2dEmbedded(:) = passiveTracer2dStage2(:) else ! SSPRK3(4,3): use stage-2 endpoint as SSP-consistent lower-order companion. layerThicknessEmbedded(:,:) = layerThicknessStage2(:,:) From dad6e5b4e058c2bb61ee80b5e25add45beae5a25 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 13 Mar 2026 15:34:44 -0600 Subject: [PATCH 18/19] Remove declaration of unused variable Remove declaration of unused variable Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 1 - 1 file changed, 1 deletion(-) 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 a54d094fc972..c7eaaaa08714 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 @@ -1692,7 +1692,6 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel type (MPAS_Time_type) :: nextForceTime, currTime type (MPAS_TimeInterval_type) :: intervalToNextForceTime real (kind=RKIND) :: secondsToNextForceTime - real (kind=RKIND) :: fracSecondsToNextForceTime real (kind=RKIND) :: allowableDt real (kind=RKIND) :: proposedDt integer :: err_tmp From 9be8aa0e8df5c3fe8fdcc13320504cc90cea246c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 13 Mar 2026 14:51:00 -0700 Subject: [PATCH 19/19] Weight error calculation by cellArea Weight error calculation in embedded RK method by cellArea --- .../mpas_li_time_integration_fe_rk.F | 62 +++++++++---------- 1 file changed, 28 insertions(+), 34 deletions(-) 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 c7eaaaa08714..5aca473bb557 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 @@ -155,7 +155,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) fluxAcrossGroundingLineAccum, & fluxAcrossGroundingLineOnCellsAccum, & calvingThicknessAccum - real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions + real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions, areaCell integer, dimension(:), allocatable :: cellMaskPrev, edgeMaskPrev, vertexMaskPrev ! cell mask before advection real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & @@ -198,8 +198,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) logical :: solveVeloAfterCalving logical :: useEmbeddedSSPRK, includeThermalForError, includeDamageForError, stepAccepted integer :: stepAttempt, maxStepAttempts, embeddedHighOrder - integer :: rkErrorCountLocal, rkErrorCountGlobal real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, rkErrorSumSqLocal, rkErrorSumSqGlobal, dtScale + real (kind=RKIND) :: rkErrorWeightLocal, rkErrorWeightGlobal real (kind=RKIND) :: rk2EmbeddedW1, rk2EmbeddedW2 real (kind=RKIND), save :: previousEmbeddedAcceptedDt = -1.0_RKIND @@ -215,10 +215,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) embeddedHighOrder = 3 rkErrorNorm = 0.0_RKIND rkErrorNormGlobal = 0.0_RKIND - rkErrorCountLocal = 0 - rkErrorCountGlobal = 0 rkErrorSumSqLocal = 0.0_RKIND rkErrorSumSqGlobal = 0.0_RKIND + rkErrorWeightLocal = 0.0_RKIND + rkErrorWeightGlobal = 0.0_RKIND rk2EmbeddedW1 = 0.694021459207626_RKIND rk2EmbeddedW2 = 0.305978540792374_RKIND dtScale = 1.0_RKIND @@ -253,6 +253,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) call mpas_pool_get_array(meshPool, 'deltat', deltat) call mpas_pool_get_array(meshPool, 'layerCenterSigma', layerCenterSigma) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) @@ -787,16 +788,15 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & - damage, damageEmbedded, nVertLevels, nCells, includeThermalForError, & + damage, damageEmbedded, areaCell, nVertLevels, nCells, includeThermalForError, & includeDamageForError, config_embedded_rk_rtol, & config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, & - rkErrorNorm, rkErrorCountLocal) + rkErrorSumSqLocal, rkErrorWeightLocal) - rkErrorSumSqLocal = rkErrorNorm * rkErrorNorm * real(rkErrorCountLocal, RKIND) call mpas_dmpar_sum_real(domain % dminfo, rkErrorSumSqLocal, rkErrorSumSqGlobal) - call mpas_dmpar_sum_int(domain % dminfo, rkErrorCountLocal, rkErrorCountGlobal) - if (rkErrorCountGlobal > 0) then - rkErrorNormGlobal = sqrt(rkErrorSumSqGlobal / real(rkErrorCountGlobal, RKIND)) + call mpas_dmpar_sum_real(domain % dminfo, rkErrorWeightLocal, rkErrorWeightGlobal) + if (rkErrorWeightGlobal > 0.0_RKIND) then + rkErrorNormGlobal = sqrt(rkErrorSumSqGlobal / rkErrorWeightGlobal) else rkErrorNormGlobal = 0.0_RKIND endif @@ -1914,32 +1914,32 @@ end subroutine calculate_layerThicknessEdge !> \brief Compute normalized error estimate for embedded SSPRK control !> \details !> Uses the current RK endpoint and a formally reconstructed embedded -!> lower-order endpoint to compute an RMSE-like normalized error over +!> lower-order endpoint to compute an area-weighted normalized error over !> thickness and tracer-like fields. !----------------------------------------------------------------------- subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & - damage, damageEmbedded, nVertLevels, nCells, includeThermal, includeDamage, & - relTol, absTolThickness, absTolTracer, errorNorm, errorCount) + damage, damageEmbedded, areaCell, nVertLevels, nCells, includeThermal, & + includeDamage, relTol, absTolThickness, absTolTracer, errorSumSq, errorWeight) real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness, layerThicknessEmbedded real (kind=RKIND), dimension(:,:), intent(in) :: temperature, temperatureEmbedded real (kind=RKIND), dimension(:,:), intent(in) :: waterFrac, waterFracEmbedded real (kind=RKIND), dimension(:), intent(in) :: passiveTracer2d, passiveTracer2dEmbedded real (kind=RKIND), dimension(:), intent(in) :: damage, damageEmbedded + real (kind=RKIND), dimension(:), intent(in) :: areaCell integer, intent(in) :: nVertLevels, nCells logical, intent(in) :: includeThermal, includeDamage real (kind=RKIND), intent(in) :: relTol, absTolThickness, absTolTracer - real (kind=RKIND), intent(out) :: errorNorm - integer, intent(out) :: errorCount + real (kind=RKIND), intent(out) :: errorSumSq + real (kind=RKIND), intent(out) :: errorWeight integer :: iCell, k - real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr, errorSumSq + real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr real (kind=RKIND), parameter :: minIceForError = 1.0e-12_RKIND - errorNorm = 0.0_RKIND - errorCount = 0 errorSumSq = 0.0_RKIND + errorWeight = 0.0_RKIND do iCell = 1, nCells yPrimary = sum(layerThickness(:, iCell)) @@ -1947,8 +1947,8 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, scaleVal = absTolThickness + relTol * max(abs(yPrimary), abs(yEmbedded)) if (scaleVal > 0.0_RKIND) then thisErr = abs(yPrimary - yEmbedded) / scaleVal - errorSumSq = errorSumSq + thisErr * thisErr - errorCount = errorCount + 1 + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) endif if (max(yPrimary, yEmbedded) <= minIceForError) then @@ -1958,16 +1958,16 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, scaleVal = absTolTracer + relTol * max(abs(passiveTracer2d(iCell)), abs(passiveTracer2dEmbedded(iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(passiveTracer2d(iCell) - passiveTracer2dEmbedded(iCell)) / scaleVal - errorSumSq = errorSumSq + thisErr * thisErr - errorCount = errorCount + 1 + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) endif if (includeDamage) then scaleVal = absTolTracer + relTol * max(abs(damage(iCell)), abs(damageEmbedded(iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(damage(iCell) - damageEmbedded(iCell)) / scaleVal - errorSumSq = errorSumSq + thisErr * thisErr - errorCount = errorCount + 1 + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) endif endif @@ -1977,26 +1977,20 @@ subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, scaleVal = absTolTracer + relTol * max(abs(temperature(k,iCell)), abs(temperatureEmbedded(k,iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(temperature(k,iCell) - temperatureEmbedded(k,iCell)) / scaleVal - errorSumSq = errorSumSq + thisErr * thisErr - errorCount = errorCount + 1 + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) endif scaleVal = absTolTracer + relTol * max(abs(waterFrac(k,iCell)), abs(waterFracEmbedded(k,iCell))) if (scaleVal > 0.0_RKIND) then thisErr = abs(waterFrac(k,iCell) - waterFracEmbedded(k,iCell)) / scaleVal - errorSumSq = errorSumSq + thisErr * thisErr - errorCount = errorCount + 1 + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) endif enddo endif enddo - if (errorCount > 0) then - errorNorm = sqrt(errorSumSq / real(errorCount, RKIND)) - else - errorNorm = 0.0_RKIND - endif - end subroutine li_embedded_ssprk_error_norm !***********************************************************************