Skip to content

Commit

Permalink
Update groundedToFloatingThickness when masks are updated
Browse files Browse the repository at this point in the history 
Update groundedToFloatingThickness whenever masks are updated to
account for cells changing between masks when face-melting and calving
are applied, and not just when sfcMassBal, basalMassBal, and advection
are applied, as was the case previously.
  • Loading branch information
@trhille
trhille committed May 19, 2022
1 parent 616feb3 commit 3c4ba20
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 83 deletions.
82 changes: 0 additions & 82 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,11 +184,6 @@ subroutine li_advection_thickness_tracers(&
surfaceTracersField, & ! scratch field containing values of surface tracers
basalTracersField ! scratch field containing values of basal tracers

type (field1DInteger), pointer :: &
groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget
floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget


integer, dimension(:), pointer :: &
cellMask, & ! integer bitmask for cells
edgeMask, & ! integer bitmask for edges
Expand Down Expand Up @@ -328,23 +323,13 @@ subroutine li_advection_thickness_tracers(&
call mpas_allocate_scratch_field(basalTracersField, .true.)
basalTracers => basalTracersField % array

call mpas_pool_get_field(geometryPool, 'groundedMaskForMassBudgetTemporary', groundedMaskForMassBudgetTemporaryField)
call mpas_pool_get_field(geometryPool, 'floatingMaskForMassBudgetTemporary', floatingMaskForMassBudgetTemporaryField)
call mpas_allocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.)
call mpas_allocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.)


! given the old thickness, compute the thickness in each layer
call li_calculate_layerThickness(meshPool, thickness, layerThickness)

! update masks
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)

! save old mass budget masks for determining cells changing from grounded to floating and vice versa
groundedMaskForMassBudgetTemporaryField % array(:) = groundedMaskForMassBudget(:)
floatingMaskForMassBudgetTemporaryField % array(:) = floatingMaskForMassBudget(:)

!-----------------------------------------------------------------
! Horizontal transport of thickness and tracers
!-----------------------------------------------------------------
Expand Down Expand Up @@ -466,8 +451,6 @@ subroutine li_advection_thickness_tracers(&
thickness = sum(layerThickness, 1)
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)



!-----------------------------------------------------------------
! Add the surface and basal mass balance to the layer thickness
!-----------------------------------------------------------------
Expand Down Expand Up @@ -540,16 +523,6 @@ subroutine li_advection_thickness_tracers(&
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, 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 grounded_to_floating(groundedMaskForMassBudgetTemporaryField % array, &
floatingMaskForMassBudgetTemporaryField % array, &
groundedMaskForMassBudget, &
floatingMaskForMassBudget, &
thickness, &
groundedToFloatingThickness, &
nCells)
! Remap tracers to the standard vertical sigma coordinate
! Note: If tracers are not being advected, then this subroutine simply restores the
! layer thickness to sigma coordinate values.
Expand Down Expand Up @@ -603,8 +576,6 @@ subroutine li_advection_thickness_tracers(&
call mpas_deallocate_scratch_field(layerThicknessOldField, .true.)
call mpas_deallocate_scratch_field(basalTracersField, .true.)
call mpas_deallocate_scratch_field(surfaceTracersField, .true.)
call mpas_deallocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.)
call mpas_deallocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.)
! === error check
if (err > 0) then
Expand Down Expand Up @@ -1934,59 +1905,6 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers
end subroutine vertical_remap
subroutine grounded_to_floating(groundedMaskForMassBudgetOrig, &
floatingMaskForMassBudgetOrig, &
groundedMaskForMassBudgetNew, &
floatingMaskForMassBudgetNew, &
thicknessNew, &
groundedToFloatingThickness, &
nCells)
!-----------------------------------------------------------------
! input variables
!-----------------------------------------------------------------
integer, dimension(:), intent(in) :: &
groundedMaskForMassBudgetOrig, & !< Input: mask for cells before advection
floatingMaskForMassBudgetOrig
integer, dimension(:), intent(in) :: &
groundedMaskForMassBudgetNew, & !< Input: mask for cells after advection
floatingMaskForMassBudgetNew
real(kind=RKIND), dimension(:), intent(in) :: &
thicknessNew !< Input: ice thickness after advection
integer, pointer, intent(in) :: &
nCells
!-----------------------------------------------------------------
! input/output variables
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! output variables
!-----------------------------------------------------------------
real(kind=RKIND), dimension(:), intent(out) :: &
groundedToFloatingThickness !< Input: ice thickness
!-----------------------------------------------------------------
! local variables
!-----------------------------------------------------------------
integer :: iCell
groundedToFloatingThickness(:) = 0.0_RKIND
do iCell = 1, nCells
! find locations that had grounded ice before but now have floating ice
if ( (groundedMaskForMassBudgetOrig(iCell) .eq. 1) .and. &
(floatingMaskForMassBudgetNew(iCell) .eq. 1) ) then
groundedToFloatingThickness(iCell) = thicknessNew(iCell)
! find locations that had floating ice before but now have grounded ice
else if ( (floatingMaskForMassBudgetOrig(iCell) .eq. 1) .and. &
(groundedMaskForMassBudgetNew(iCell) .eq. 1) ) then
groundedToFloatingThickness(iCell) = -1.0_RKIND * thicknessNew(iCell)
endif
enddo
end subroutine grounded_to_floating
!***********************************************************************
end module li_advection
Expand Down
6 changes: 6 additions & 0 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,7 @@ subroutine prepare_advection(domain, err)

real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity
real (kind=RKIND), dimension(:), pointer :: groundedToFloatingThickness
real (kind=RKIND), pointer :: calvingCFLdt
integer, pointer :: processLimitingTimestep
integer, dimension(:), pointer :: edgeMask
Expand Down Expand Up @@ -374,8 +375,13 @@ subroutine prepare_advection(domain, err)

call mpas_pool_get_array(velocityPool, 'normalVelocity', normalVelocity)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)
call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity)

! groundedToFloatingThickness is updated throughout
! the timestep, so necessary to zero it at
! the beginning of every timestep.
groundedToFloatingThickness(:) = 0.0_RKIND
! compute normal velocities and advective CFL limit for this block

call li_layer_normal_velocity( &
Expand Down
28 changes: 27 additions & 1 deletion components/mpas-albany-landice/src/shared/mpas_li_mask.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err)
!-----------------------------------------------------------------
integer, pointer :: nCells, nVertices, nEdges, vertexDegree
integer, pointer :: nVertInterfaces
real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography
real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography, groundedToFloatingThickness
integer, dimension(:), pointer :: nEdgesOnCell, cellMask, vertexMask, edgeMask, &
groundedMaskForMassBudget, floatingMaskForMassBudget
integer, dimension(:,:), pointer :: cellsOnCell, cellsOnVertex, cellsOnEdge, dirichletVelocityMask
Expand All @@ -310,6 +310,9 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err)
integer :: iCellNeighbor
logical :: openOceanNeighbor
logical :: updateMassBudgetMasks=.true.
type (field1DInteger), pointer :: & ! used to calculate groundedToFloatingThickness
groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget
floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget

call mpas_timer_start('calculate mask')

Expand All @@ -332,9 +335,15 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)

call mpas_pool_get_field(geometryPool, 'groundedMaskForMassBudgetTemporary', groundedMaskForMassBudgetTemporaryField)
call mpas_pool_get_field(geometryPool, 'floatingMaskForMassBudgetTemporary', floatingMaskForMassBudgetTemporaryField)
call mpas_allocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.)
call mpas_allocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.)

call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density)
call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density)
call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
Expand All @@ -360,6 +369,11 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err)
cellMask = ior(cellMask, li_mask_ValueIce)
end where


! save old mass budget masks for determining cells changing from grounded
! to floating and vice versa
groundedMaskForMassBudgetTemporaryField % array(:) = groundedMaskForMassBudget(:)
floatingMaskForMassBudgetTemporaryField % array(:) = floatingMaskForMassBudget(:)
! Is it floating? (ice thickness equal to floatation is considered floating)
! For now floating ice and grounded ice are mutually exclusive.
! This may change if a grounding line parameterization is added.
Expand Down Expand Up @@ -518,6 +532,18 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err)
endif
enddo

! Update groundedToFloatingThickness based on new masks
do iCell = 1, nCells
! find locations that had grounded ice before but now have floating ice
if ( (groundedMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. &
(floatingMaskForMassBudget(iCell) .eq. 1) ) then
groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) + thickness(iCell)
! find locations that had floating ice before but now have grounded ice
else if ( (floatingMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. &
(groundedMaskForMassBudget(iCell) .eq. 1) ) then
groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) - thickness(iCell)
endif
enddo
! Now use flood fill to identify all floating ice connected to the open
! ocean. Subglacial lakes are any floating ice that are not in the
! resulting flood fill mask. Commented code below does not yet work
Expand Down

0 comments on commit 3c4ba20

Please sign in to comment.