Skip to content

Commit

Permalink
Make grounding line flux calculation consistent with FCT
Browse files Browse the repository at this point in the history
Make grounding line flux calculation consistent with FCT advection.
Also update edgeMask halos before calculating grounding line flux.
  • Loading branch information
trhille committed Oct 20, 2023
1 parent 8fe44f0 commit b43726e
Showing 1 changed file with 22 additions and 15 deletions.
37 changes: 22 additions & 15 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
Original file line number Diff line number Diff line change
Expand Up @@ -369,38 +369,45 @@ subroutine li_advection_thickness_tracers(&
! update masks
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
! Update halo on edgeMask for calculating flux across grounding line.
call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'edgeMask')
call mpas_timer_stop("halo updates")

layerThicknessEdgeFlux(:,:) = 0.0_RKIND
!-----------------------------------------------------------------
! Horizontal transport of thickness and tracers
!-----------------------------------------------------------------

if ( (trim(config_thickness_advection) .eq. 'fo') .or. &
(trim(config_thickness_advection) .eq. 'fct') ) then

! Calculate flux across grounding line before applying
! sfc/basal mass balance and advection
! Calculate flux across grounding line
! Do this after new thickness & mask have been calculated, including SMB/BMB
fluxAcrossGroundingLine(:) = 0.0_RKIND
fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND
do iEdge = 1, nEdges
if (li_mask_is_grounding_line(edgeMask(iEdge))) then
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
if (li_mask_is_grounded_ice(cellMask(iCell1))) then
GLfluxSign = 1.0_RKIND ! edge sign convention is positive from iCell1 to iCell2 on an edge
theGroundedCell = iCell1
else
else
GLfluxSign = -1.0_RKIND
theGroundedCell = iCell2
endif
do k = 1, nVertLevels
thicknessFluxEdge = layerNormalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge)
fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * thicknessFluxEdge / dvEdge(iEdge)
enddo
if (trim(config_thickness_advection) == 'fct') then
do k = 1, nVertLevels
fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge)
enddo
else
do k = 1, nVertLevels
layerThicknessEdgeFlux(k,iEdge) = layerNormalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge)
fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge)
enddo
endif
! assign to grounded cell in fluxAcrossGroundingLineOnCells
if (thickness(theGroundedCell) <= 0.0_RKIND) then
! This should never be the case, but checking to avoid
! possible divide by zero
if (thickness(theGroundedCell) <= 0.0_RKIND) then
! This should never be the case, but checking to avoid possible divide by zero
call mpas_log_write("thickness at a grounding line is unexepectedly <=0", MPAS_LOG_ERR)
err = ior(err, 1)
return
Expand Down Expand Up @@ -474,7 +481,7 @@ subroutine li_advection_thickness_tracers(&
! copy old (input) values of layer thickness and tracers to scratch arrays
layerThicknessOld(:,:) = layerThickness(:,:)
advectedTracersOld(:,:,:) = advectedTracers(:,:,:)

layerThicknessEdgeFlux(:,:) = 0.0_RKIND
! compute new values of layer thickness and tracers
if ((trim(config_thickness_advection) .eq. 'fo') .and. &
( (trim(config_tracer_advection) .eq. 'fo') .or. &
Expand Down

0 comments on commit b43726e

Please sign in to comment.