From 51802e1bca96f70da0a72ab6786116b88dfa9c67 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 6 Nov 2023 11:47:15 -0800 Subject: [PATCH] Average fluxAcrossGroundingLineOnCells in time integration module Calculate fluxAcrossGroundingLineOnCells in advection module for each RK stage and then take weighted average after RK integration instead of calcuating fluxAcrossGroundingLineOnCells on final state. This is more consistent with how fluxAcrossGroundingLine (on edges) is treated. --- .../src/mode_forward/mpas_li_advection.F | 4 -- .../mpas_li_time_integration_fe_rk.F | 37 ++++++++----------- 2 files changed, 15 insertions(+), 26 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F index 5e0b90846735..56b018cdd399 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F @@ -630,10 +630,6 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) - ! Calculate volume converted from grounded to floating - ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state - call li_grounded_to_floating(cellMaskTemporaryField % array, cellMask, thickness, groundedToFloatingThickness, nCells) - ! Calculate flux across grounding line ! Do this after new thickness & mask have been calculated, including SMB/BMB fluxAcrossGroundingLine(:) = 0.0_RKIND 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 5dda3f792766..6c9874f26486 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 @@ -135,7 +135,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum, & groundedBasalMassBalAccum, & floatingBasalMassBalAccum, & - fluxAcrossGroundingLineAccum + fluxAcrossGroundingLineAccum, & + fluxAcrossGroundingLineOnCellsAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & @@ -202,6 +203,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(groundedBasalMassBalAccum(nCells+1)) allocate(floatingBasalMassBalAccum(nCells+1)) allocate(fluxAcrossGroundingLineAccum(nEdges+1)) + allocate(fluxAcrossGroundingLineOnCellsAccum(nCells+1)) temperaturePrev(:,:) = 0.0_RKIND waterFracPrev(:,:) = 0.0_RKIND @@ -219,6 +221,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedBasalMassBalAccum(:) = 0.0_RKIND floatingBasalMassBalAccum(:) = 0.0_RKIND fluxAcrossGroundingLineAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineOnCellsAccum(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! @@ -454,6 +457,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) 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 @@ -463,6 +467,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) 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") @@ -490,8 +495,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo ! Finalize budget updates - ! Calculate volume converted from grounded to floating - ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state + ! Update masks after RK integration + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + err = ior(err, err_tmp) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) @@ -500,32 +507,18 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'edgeMask') call mpas_timer_stop("halo updates") + ! Calculate volume converted from grounded to floating + ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state call li_grounded_to_floating(cellMaskPrev, cellMask, thickness, groundedToFloatingThickness, nCells) + sfcMassBalApplied(:) = sfcMassBalAccum(:) groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) basalMassBalApplied(:) = basalMassBalAccum(:) groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) fluxAcrossGroundingLine(:) = fluxAcrossGroundingLineAccum(:) - - fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND - do iEdge = 1, nEdges - if (li_mask_is_grounding_line(edgeMask(iEdge))) then - iCell1 = cellsOnEdge(1,iEdge) - iCell2 = cellsOnEdge(2,iEdge) - if (li_mask_is_grounded_ice(cellMask(iCell1))) then - theGroundedCell = iCell1 - else - theGroundedCell = iCell2 - endif - if ( thickness(theGroundedCell) > 0.0_RKIND ) then - fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + & - fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density ! adjust to correct units - else - fluxAcrossGroundingLineOnCells(theGroundedCell) = 0.0_RKIND - endif - endif - enddo ! edges + fluxAcrossGroundingLineOnCells(:) = fluxAcrossGroundingLineOnCellsAccum(:) + ! Reset time step to full length after RK loop deltat = deltatFull