Skip to content

Commit

Permalink
Account for multiple updates of calvingThickness within a timestep
Browse files Browse the repository at this point in the history 
Add a very short subroutine update_calving_budget that is called to
recalculate groundedCalvingThickness and floatingCalvingThickness
after calvingThickness is applied and before masks are updated. Also
treat grounded and floating calving specially within the remove_icebergs
and remove_small_islands routines to deal with the face that
calvingThickness is updated multiple times per timestep. Thus, it
is okay for a cell to have both nonzero groundedCalvingThickness
and floatingCalvingThickness.
  • Loading branch information
@trhille
trhille committed May 17, 2022
1 parent a16814e commit a6c5c9d
Showing 1 changed file with 108 additions and 40 deletions.
148 changes: 108 additions & 40 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,11 @@ subroutine li_calve_ice(domain, err)
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
calvingThickness(:) = 0.0_RKIND
groundedCalvingThickness(:) = 0.0_RKIND
floatingCalvingThickness(:) = 0.0_RKIND

block => block % next
end do
Expand Down Expand Up @@ -288,37 +292,13 @@ subroutine li_calve_ice(domain, err)
! now also remove any icebergs
call remove_icebergs(domain)

block => domain % blocklist
do while (associated(block))
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)

groundedCalvingThickness(:) = 0.0_RKIND
floatingCalvingThickness(:) = 0.0_RKIND
where (groundedMaskForMassBudget .eq. 1)
groundedCalvingThickness = calvingThickness
elsewhere (floatingMaskForMassBudget .eq. 1)
floatingCalvingThickness = calvingThickness
elsewhere
groundedCalvingThickness = 0.0_RKIND
floatingCalvingThickness = 0.0_RKIND
end where

block => block % next
end do

! Final operations after calving has been applied.
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)

! In data calving mode we just calculate what should be calved but don't actually calve it.
Expand Down Expand Up @@ -596,6 +576,7 @@ subroutine li_restore_calving_front(domain, err)
block => block % next
enddo
call update_calving_budget(domain)
! Update mask and geometry
block => domain % blocklist
do while (associated(block))
Expand Down Expand Up @@ -895,6 +876,7 @@ subroutine thickness_calving(domain, calvingFraction, err)
! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)
call update_calving_budget(domain)
call remove_small_islands(meshPool, geometryPool)
block => block % next
Expand Down Expand Up @@ -968,6 +950,7 @@ subroutine floating_calving(domain, calvingFraction, err)
! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)
call update_calving_budget(domain)
call remove_small_islands(meshPool, geometryPool)
block => block % next
Expand All @@ -994,10 +977,14 @@ subroutine remove_small_islands(meshPool, geometryPool)
logical, pointer :: config_remove_small_islands
real(kind=RKIND), pointer :: config_sea_level
real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine)
real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine)
groundedCalvingThickness, &
floatingCalvingThickness
real (kind=RKIND), dimension(:), pointer :: thickness
real (kind=RKIND), dimension(:), pointer :: bedTopography
integer, dimension(:), pointer :: cellMask
integer, dimension(:), pointer :: cellMask, &
groundedMaskForMassBudget, &
floatingMaskForMassBudget
integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell
integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell
integer, pointer :: nCellsSolve
Expand All @@ -1014,6 +1001,10 @@ subroutine remove_small_islands(meshPool, geometryPool)
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
Expand All @@ -1037,6 +1028,12 @@ subroutine remove_small_islands(meshPool, geometryPool)
! If this is a single cell of ice surrounded by open ocean, kill this location
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell)
thickness(iCell) = 0.0_RKIND
! Update calving budgets
if (groundedMaskForMassBudget(iCell) .eq. 1) then
groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell)
elseif (floatingMaskForMassBudget(iCell) .eq. 1) then
floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell)
endif
elseif (nIceNeighbors == 1) then
! check if this neighbor has any additional neighbors with ice
nIceNeighbors2 = 0
Expand All @@ -1055,8 +1052,21 @@ subroutine remove_small_islands(meshPool, geometryPool)
! kill both cells
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell)
thickness(iCell) = 0.0_RKIND
! Update calving budgets
if (groundedMaskForMassBudget(iCell) .eq. 1) then
groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell)
elseif (floatingMaskForMassBudget(iCell) .eq. 1) then
floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell)
endif
calvingThickness(neighborWithIce) = calvingThickness(neighborWithIce) + thickness(neighborWithIce)
thickness(neighborWithIce) = 0.0_RKIND
! Update calving budgets
if (groundedMaskForMassBudget(neighborWithIce) .eq. 1) then
groundedCalvingThickness(neighborWithIce) = groundedCalvingThickness(neighborWithIce) + thickness(neighborWithIce)
elseif (floatingMaskForMassBudget(neighborWithIce) .eq. 1) then
floatingCalvingThickness(neighborWithIce) = floatingCalvingThickness(neighborWithIce) + thickness(neighborWithIce)
endif
endif
endif ! check on nIceNeighbors
Expand Down Expand Up @@ -1137,6 +1147,7 @@ subroutine topographic_calving(domain, calvingFraction, err)
! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)
call update_calving_budget(domain)
call remove_small_islands(meshPool, geometryPool)
block => block % next
Expand Down Expand Up @@ -1316,6 +1327,7 @@ subroutine eigencalving(domain, err)
call mpas_timer_stop("halo updates")
! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)
call update_calving_budget(domain)

! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
Expand All @@ -1340,7 +1352,7 @@ subroutine eigencalving(domain, err)
endif
enddo
! TODO: global reduce & reporting on amount of calving generated in this step

call update_calving_budget(domain)
! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
Expand All @@ -1365,7 +1377,7 @@ subroutine eigencalving(domain, err)
endif
enddo
! TODO: global reduce & reporting on amount of calving generated in this step
call update_calving_budget(domain)
call remove_small_islands(meshPool, geometryPool)
block => block % next
Expand Down Expand Up @@ -1494,7 +1506,7 @@ subroutine specified_calving_velocity(domain, err)
! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)
call update_calving_budget(domain)
! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
Expand All @@ -1518,7 +1530,7 @@ subroutine specified_calving_velocity(domain, err)
endif
enddo
! TODO: global reduce & reporting on amount of calving generated in this step
call update_calving_budget(domain)
! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
Expand All @@ -1543,7 +1555,7 @@ subroutine specified_calving_velocity(domain, err)
endif
enddo
! TODO: global reduce & reporting on amount of calving generated in this step

call update_calving_budget(domain)
call remove_small_islands(meshPool, geometryPool)

block => block % next
Expand Down Expand Up @@ -1821,13 +1833,13 @@ subroutine von_Mises_calving(domain, err)

! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)

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

call remove_small_islands(meshPool, geometryPool)

block => block % next

enddo ! associated(block)
Expand Down Expand Up @@ -2086,7 +2098,7 @@ subroutine ismip6_retreat(domain, err)

! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)

call update_calving_budget(domain)
! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
Expand Down Expand Up @@ -2999,7 +3011,7 @@ subroutine damage_calving(domain, err)

! === apply calving ===
thickness(:) = thickness(:) - calvingThickness(:)

call update_calving_budget(domain)
! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
Expand All @@ -3015,7 +3027,7 @@ subroutine damage_calving(domain, err)
thickness(iCell) = 0.0_RKIND
endif
enddo

call update_calving_budget(domain)
! update mask
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
err = ior(err, err_tmp)
Expand All @@ -3040,9 +3052,8 @@ subroutine damage_calving(domain, err)
endif
enddo
! TODO: global reduce & reporting on amount of calving generated in this step
call update_calving_budget(domain)
call remove_small_islands(meshPool, geometryPool)
block => block % next
enddo
Expand Down Expand Up @@ -3790,9 +3801,13 @@ subroutine remove_icebergs(domain)
integer, dimension(:), pointer :: seedMask, growMask !masks to pass to flood-fill routine
!integer, dimension(:), allocatable :: seedMaskOld !mask of where icebergs are removed, for debugging
real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine)
real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine)
groundedCalvingThickness, &
floatingCalvingThickness
real (kind=RKIND), dimension(:), pointer :: thickness
integer, dimension(:), pointer :: cellMask
integer, dimension(:), pointer :: cellMask, &
groundedMaskForMassBudget, &
floatingMaskForMassBudget
integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell
integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell
Expand Down Expand Up @@ -3877,6 +3892,10 @@ subroutine remove_icebergs(domain)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
call mpas_pool_get_dimension(geometryPool, 'nCells', nCells)
call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve)
!allocate(seedMaskOld(nCells+1)) ! debug: make this a mask of where icebergs were removed
Expand All @@ -3887,6 +3906,11 @@ subroutine remove_icebergs(domain)
if (seedMask(iCell) == 0 .and. li_mask_is_floating_ice(cellMask(iCell))) then
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) ! remove any remaining ice here
thickness(iCell) = 0.0_RKIND
if (groundedMaskForMassBudget(iCell) .eq. 1) then
groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell)
elseif (floatingMaskForMassBudget(iCell) .eq. 1) then
floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell)
endif
localIcebergCellCount = localIcebergCellCount + 1
!seedMaskOld(iCell) = 1 ! debug: make this a mask of where icebergs were removed
endif
Expand Down Expand Up @@ -4063,6 +4087,50 @@ subroutine li_flood_fill(seedMask, growMask, domain)
end subroutine li_flood_fill
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! routine update_calving_budget
!
!> \brief Keep a running total of calving applied to grounded and floating
!> ice, respectively.
!> \author Trevor Hillebrand
!> \date May 2022
!> \details This routine should be called after each time time calvingThickness
!> is applied, but before masks are updated which often happens multiple times
!> in a timestep.
!-----------------------------------------------------------------------
subroutine update_calving_budget(domain)
!-----------------------------------------------------------------
! input/output variables
!-----------------------------------------------------------------
type (domain_type), intent(inout) :: domain
!-----------------------------------------------------------------
! local variables
!-----------------------------------------------------------------
type (mpas_pool_type), pointer :: meshPool, geometryPool
integer, dimension(:), pointer :: groundedMaskForMassBudget, & ! binary masks for mass budget
floatingMaskForMassBudget
real (kind=RKIND), dimension(:), pointer :: calvingThickness, &
groundedCalvingThickness, & ! Grounded and floating components for mass budget
floatingCalvingThickness
call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
where (groundedMaskForMassBudget .eq. 1)
groundedCalvingThickness = calvingThickness
elsewhere (floatingMaskForMassBudget .eq. 1)
floatingCalvingThickness = calvingThickness
end where
end subroutine update_calving_budget
end module li_calving

0 comments on commit a6c5c9d

Please sign in to comment.