Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Redefine grounding line mask #37

Merged
merged 6 commits into from
Oct 24, 2022
Merged

Redefine grounding line mask #37

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
bcca081
Update grounding line mask definition for cells, edges, and vertices
trhille May 12, 2022
237a3a0
Update logic for ablation velocity after redefining grounding line
trhille May 12, 2022
b01ceed
Do not melt floating non-dynamic fringe
trhille Oct 20, 2022
1bb4c76
Changes to hydro model to comply with new GL mask change
matthewhoffman Oct 24, 2022
5b8d720
Change two instances of cycle to exit
matthewhoffman Oct 24, 2022
4f8f918
Update comments with month of grounding line mask update
trhille Oct 24, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,8 @@
description="total water flux in subglacial hydrology system" />
<var name="waterFluxMask" type="integer" dimensions="nEdges Time" units="none"
description="mask indicating how to handle fluxes on each edge: 0=calculate based on hydropotential gradient; 1=allow outflow based on hydropotential gradient, but no inflow (NOT YET IMPLEMENTED); 2=zero flux" />
<var name="hydroMarineMarginMask" type="integer" dimensions="nEdges Time" units="none"
description="mask indicating the marine boundary of the active subglacial hydrology domain" />
<var name="waterFluxAdvec" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
description="advective water flux in subglacial hydrology system" />
<var name="waterFluxDiffu" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
Expand Down
6 changes: 4 additions & 2 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2246,8 +2246,10 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
! margin, (2) have at least one neighboring cell without ice, (3) contain
! grounded ice, and (4) have bed topography below sea level.
! OR is adjacent to an inactive floating margin cell
if (li_mask_is_grounding_line(cellMask(iCell)) &
!< GL here means cell is grounded but has a neighbor that is floating or ocean
if (li_mask_is_grounded_ice(cellMask(iCell)) &
! GL here means cell is grounded but has a neighbor that is floating or ocean.
! Note that as of Oct 2022, this is no longer the general
! definition of the grounding line throughout the code.
.and. bedTopography(iCell) < config_sea_level &
.and. li_mask_is_dynamic_ice(cellMask(iCell)) ) then

Expand Down
17 changes: 10 additions & 7 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
nEdgesOnCell, &
edgeMask

integer, dimension(:,:), pointer :: edgesOnCell
integer, dimension(:,:), pointer :: edgesOnCell, cellsOnCell

type (field1dInteger), pointer :: thermalCellMaskField

Expand All @@ -275,7 +275,7 @@ subroutine li_basal_melt_floating_ice(domain, err)

real(kind=RKIND), pointer :: daysSinceStart

integer :: iCell, iEdge, err_tmp
integer :: iCell, jCell, iEdge, iNeighbor, err_tmp

! Local variables for some melt methods

Expand Down Expand Up @@ -509,23 +509,26 @@ subroutine li_basal_melt_floating_ice(domain, err)
call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)

! If config_front_mass_bal_grounded is not none, only apply ice shelf melt to active cells
! and stranded non-dynamic cells.
! and stranded non-dynamic cells. i.e., if a floating non-dynamic
! cell has a grounded neighbor, do not apply ice shelf melt to it.
do iCell = 1, nCellsSolve
if ( li_mask_is_floating_ice(cellMask(iCell)) .and. (.not. li_mask_is_dynamic_ice(cellMask(iCell))) ) then
do iEdge = 1, nEdgesOnCell(iCell)
if ( li_mask_is_grounding_line(edgeMask(edgesOnCell(iEdge,iCell))) ) then

do iNeighbor = 1, nEdgesOnCell(iCell)
jCell = cellsOnCell(iNeighbor, iCell)
if ( li_mask_is_grounded_ice(cellMask(jCell)) ) then
floatingBasalMassBal(iCell) = 0.0_RKIND

exit ! No need to look over other neighbors
endif
enddo
endif
enddo

block => block % next
enddo ! associated(block)
endif
Expand Down
120 changes: 112 additions & 8 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,8 @@ subroutine li_SGH_init(domain, err)
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
err = ior(err, err_tmp)

call calc_hydro_mask(domain)

! remove invalid values - not necessary on restart, but shouldn't hurt
call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
waterThickness = max(0.0_RKIND, waterThickness)
Expand Down Expand Up @@ -294,6 +296,8 @@ subroutine li_SGH_solve(domain, err)
return
endif

call calc_hydro_mask(domain)

call mpas_log_write('Beginning subglacial hydro solve.')
call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_active', config_SGH_chnl_active)
call mpas_pool_get_config(liConfigs, 'config_SGH_till_drainage', Cd)
Expand Down Expand Up @@ -728,6 +732,7 @@ subroutine calc_edge_quantities(block, err)
real (kind=RKIND), dimension(:), pointer :: waterFluxAdvec
real (kind=RKIND), dimension(:), pointer :: waterFluxDiffu
integer, dimension(:), pointer :: waterFluxMask
integer, dimension(:), pointer :: hydroMarineMarginMask
integer, dimension(:,:), pointer :: edgeSignOnCell
integer, dimension(:), pointer :: cellMask
integer, dimension(:), pointer :: edgeMask
Expand Down Expand Up @@ -797,6 +802,7 @@ subroutine calc_edge_quantities(block, err)
call mpas_pool_get_array(hydroPool, 'waterFluxAdvec', waterFluxAdvec)
call mpas_pool_get_array(hydroPool, 'waterFluxDiffu', waterFluxDiffu)
call mpas_pool_get_array(hydroPool, 'waterFluxMask', waterFluxMask)
call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)


Expand Down Expand Up @@ -827,7 +833,7 @@ subroutine calc_edge_quantities(block, err)
! at the edge of the cell in a 1-sided sense
do iEdge = 1, nEdges
if ( (li_mask_is_margin(edgeMask(iEdge)) .and. li_mask_is_grounded_ice(edgeMask(iEdge))) .or. &
(li_mask_is_grounding_line(edgeMask(iEdge)))) then
(hydroMarineMarginMask(iEdge)==1)) then
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
if (li_mask_is_grounded_ice(cellMask(cell1))) then ! cell2 is the icefree cell - replace phi there with cell1 Phig
Expand Down Expand Up @@ -856,7 +862,7 @@ subroutine calc_edge_quantities(block, err)
! (Do this step only after the other hydropotential special cases are treated above.)
do iEdge = 1, nEdges
! Find edges along GL or margin to check for 'backwards' flow
if ((li_mask_is_grounding_line(edgeMask(iEdge))) .or. &
if ((hydroMarineMarginMask(iEdge)==1) .or. &
li_mask_is_margin(edgeMask(iEdge)) ) then
! Now check if flow is backwards
cell1 = cellsOnEdge(1, iEdge)
Expand Down Expand Up @@ -1364,7 +1370,10 @@ subroutine calc_pressure(block, err)
real (kind=RKIND), dimension(:), pointer :: divergenceChannel
real (kind=RKIND), dimension(:), pointer :: channelAreaChangeCell
real (kind=RKIND), dimension(:), pointer :: bedTopography
integer, dimension(:), pointer :: hydroMarineMarginMask
integer, dimension(:), pointer :: cellMask
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnCell
real (kind=RKIND), pointer :: deltatSGH
real (kind=RKIND), pointer :: bedRough, bedRoughMax
real (kind=RKIND), pointer :: rhoi
Expand All @@ -1377,6 +1386,8 @@ subroutine calc_pressure(block, err)
real (kind=RKIND), pointer :: rhoo
integer :: err_tmp
integer :: iCell
integer :: iEdge
logical :: onMarineMargin

err = 0
err_tmp = 0
Expand All @@ -1400,6 +1411,8 @@ subroutine calc_pressure(block, err)

call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)

call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
Expand All @@ -1417,6 +1430,7 @@ subroutine calc_pressure(block, err)
call mpas_pool_get_array(hydroPool, 'divergence', divergence)
call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel)
call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell)
call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed)
call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
Expand Down Expand Up @@ -1472,11 +1486,20 @@ subroutine calc_pressure(block, err)
((.not. li_mask_is_ice(cellMask(iCell))) .and. (bedTopography(iCell) < config_sea_level) ) ) then
! set pressure correctly under floating ice and open ocean
waterPressure(iCell) = rhoo * gravity * (config_sea_level - bedTopography(iCell))
elseif (li_mask_is_grounding_line(cellMask(iCell))) then
! At GL, don't let water pressure fall below ocean pressure
! TODO: Not sure if this should include the water layer thickness term. Leaving it off.
if (waterPressure(iCell) < rhoo * gravity * (config_sea_level - bedTopography(iCell))) then
waterPressure(iCell) = rhoo * gravity * (config_sea_level - bedTopography(iCell))
else
onMarineMargin = .false.
do iEdge = 1, nEdgesOnCell(iCell)
if (hydroMarineMarginMask(edgesOnCell(iEdge, iCell)) == 1) then
onMarineMargin = .true.
exit
endif
enddo
if (onMarineMargin) then
! At marine margin, don't let water pressure fall below ocean pressure
! TODO: Not sure if this should include the water layer thickness term. Leaving it off.
if (waterPressure(iCell) < rhoo * gravity * (config_sea_level - bedTopography(iCell))) then
waterPressure(iCell) = rhoo * gravity * (config_sea_level - bedTopography(iCell))
endif
endif
endif
enddo
Expand Down Expand Up @@ -1633,6 +1656,7 @@ subroutine update_channel(block, err)
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: channelDiffusivity
integer, dimension(:), pointer :: waterFluxMask
integer, dimension(:), pointer :: hydroMarineMarginMask
integer, dimension(:), pointer :: edgeMask
real (kind=RKIND), dimension(:,:), pointer :: flowParamA
integer, dimension(:,:), pointer :: cellsOnEdge
Expand Down Expand Up @@ -1680,6 +1704,7 @@ subroutine update_channel(block, err)
call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'waterFluxMask', waterFluxMask)
call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
call mpas_pool_get_array(hydroPool, 'channelDiffusivity', channelDiffusivity)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)

Expand All @@ -1699,7 +1724,7 @@ subroutine update_channel(block, err)

! Note: an edge with only one grounded cell neighbor is called floating, so this logic retains channel vars
! on those edges to allow channel discharge across GL
where (.not. ( (li_mask_is_grounded_ice(edgeMask)) .or. (li_mask_is_grounding_line(edgeMask)) ) )
where (.not. ( (li_mask_is_grounded_ice(edgeMask)) .or. (hydroMarineMarginMask==1) ) )
channelArea = 0.0_RKIND
channelDischarge = 0.0_RKIND
end where
Expand Down Expand Up @@ -2049,4 +2074,83 @@ subroutine ocean_connection_N(domain)
end subroutine ocean_connection_N


!***********************************************************************
!
! routine calc_hydro_mask
!
!> \brief Calculate the boundaries of the active hydrology domain
!> \author Matt Hoffman
!> \date 24 October 2022
!> \details
!> This routine calculates a mask of the boundaries of the active hydrology domain
!-----------------------------------------------------------------------
subroutine calc_hydro_mask(domain)

!-----------------------------------------------------------------
! input variables
!-----------------------------------------------------------------

!-----------------------------------------------------------------
! input/output variables
!-----------------------------------------------------------------
type (domain_type), intent(inout) :: domain !< Input/Output: domain object

!-----------------------------------------------------------------
! output variables
!-----------------------------------------------------------------

!-----------------------------------------------------------------
! local variables
!-----------------------------------------------------------------
type (block_type), pointer :: block
type (mpas_pool_type), pointer :: hydroPool
type (mpas_pool_type), pointer :: geometryPool
type (mpas_pool_type), pointer :: meshPool
real (kind=RKIND), dimension(:), pointer :: bedTopography
integer, dimension(:), pointer :: hydroMarineMarginMask
integer, dimension(:,:), pointer :: cellsOnEdge
integer, dimension(:), pointer :: cellMask
integer, pointer :: nEdgesSolve
integer :: cell1, cell2, iEdge
real (kind=RKIND), pointer :: config_sea_level

call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)

block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)

hydroMarineMarginMask(:) = 0
do iEdge = 1, nEdgesSolve
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
! We are looking for edges with 1 edge grounded ice and the other edge floating ice or open ocean
if ( (li_mask_is_grounded_ice(cellMask(cell1))) .and. &
(li_mask_is_floating_ice(cellMask(cell2)) .or. &
((bedTopography(cell2) < config_sea_level) .and. (.not. li_mask_is_ice(cellMask(cell2)))) ) ) then
hydroMarineMarginMask(iEdge) = 1
elseif ( (li_mask_is_grounded_ice(cellMask(cell2))) .and. &
(li_mask_is_floating_ice(cellMask(cell1)) .or. &
((bedTopography(cell1) < config_sea_level) .and. (.not. li_mask_is_ice(cellMask(cell1)))) ) ) then
hydroMarineMarginMask(iEdge) = 1
endif
enddo

block => block % next
end do

call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'hydroMarineMarginMask')
call mpas_timer_stop("halo updates")

!--------------------------------------------------------------------
end subroutine calc_hydro_mask

end module li_subglacial_hydro
Loading