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 41614ce60f4d..5571caa865d0 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 @@ -381,47 +381,6 @@ subroutine li_advection_thickness_tracers(& if ( (trim(config_thickness_advection) .eq. 'fo') .or. & (trim(config_thickness_advection) .eq. 'fct') ) then - ! 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, nEdgesSolve - 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 - GLfluxSign = 1.0_RKIND ! edge sign convention is positive from iCell1 to iCell2 on an edge - theGroundedCell = iCell1 - else - GLfluxSign = -1.0_RKIND - theGroundedCell = iCell2 - endif - 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 - call mpas_log_write("thickness at a grounding line is unexepectedly <=0", MPAS_LOG_ERR) - err = ior(err, 1) - return - endif - fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + & - fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density ! adjust to correct units - endif - enddo ! edges - - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLine') - call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLineOnCells') - call mpas_timer_stop("halo updates") if (config_print_thickness_advection_info) then call mpas_log_write('Using first-order upwind for thickness advection') @@ -582,6 +541,48 @@ subroutine li_advection_thickness_tracers(& call mpas_dmpar_field_halo_exch(domain, 'layerThickness') call mpas_timer_stop("halo updates") + ! Calculate flux across grounding line + ! Do this before new masks and thickness have been calculated, and before SMB/BMB + fluxAcrossGroundingLine(:) = 0.0_RKIND + fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND + do iEdge = 1, nEdgesSolve + 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 + GLfluxSign = 1.0_RKIND ! edge sign convention is positive from iCell1 to iCell2 on an edge + theGroundedCell = iCell1 + else + GLfluxSign = -1.0_RKIND + theGroundedCell = iCell2 + endif + 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 + call mpas_log_write("thickness at a grounding line is unexepectedly <=0", MPAS_LOG_ERR) + err = ior(err, 1) + return + endif + fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + & + fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density ! adjust to correct units + endif + enddo ! edges + + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLine') + call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLineOnCells') + call mpas_timer_stop("halo updates") + ! Calculate dynamicThickening (layerThickness is updated by advection at this point, while thickness is still old) dynamicThickening = (sum(layerThickness, 1) - thickness) / dt * scyr ! units of m/yr