Skip to content

Commit a6c5c9d

Browse files
trhilletrhille
committed
Account for multiple updates of calvingThickness within a timestep
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.
1 parent a16814e commit a6c5c9d

File tree

1 file changed

+108
-40
lines changed

1 file changed

+108
-40
lines changed

‎components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F‎

Lines changed: 108 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -217,7 +217,11 @@ subroutine li_calve_ice(domain, err)
217217
do while (associated(block))
218218
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
219219
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
220+
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
221+
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
220222
calvingThickness(:) = 0.0_RKIND
223+
groundedCalvingThickness(:) = 0.0_RKIND
224+
floatingCalvingThickness(:) = 0.0_RKIND
221225

222226
block => block % next
223227
end do
@@ -288,37 +292,13 @@ subroutine li_calve_ice(domain, err)
288292
! now also remove any icebergs
289293
call remove_icebergs(domain)
290294

291-
block => domain % blocklist
292-
do while (associated(block))
293-
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
294-
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
295-
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
296-
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
297-
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
298-
299-
groundedCalvingThickness(:) = 0.0_RKIND
300-
floatingCalvingThickness(:) = 0.0_RKIND
301-
where (groundedMaskForMassBudget .eq. 1)
302-
groundedCalvingThickness = calvingThickness
303-
elsewhere (floatingMaskForMassBudget .eq. 1)
304-
floatingCalvingThickness = calvingThickness
305-
elsewhere
306-
groundedCalvingThickness = 0.0_RKIND
307-
floatingCalvingThickness = 0.0_RKIND
308-
end where
309-
310-
block => block % next
311-
end do
312-
313295
! Final operations after calving has been applied.
314296
block => domain % blocklist
315297
do while (associated(block))
316298
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
317299
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
318300
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
319301
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
320-
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
321-
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
322302
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
323303

324304
! In data calving mode we just calculate what should be calved but don't actually calve it.
@@ -596,6 +576,7 @@ subroutine li_restore_calving_front(domain, err)
596576
block => block % next
597577
enddo
598578
579+
call update_calving_budget(domain)
599580
! Update mask and geometry
600581
block => domain % blocklist
601582
do while (associated(block))
@@ -895,6 +876,7 @@ subroutine thickness_calving(domain, calvingFraction, err)
895876
! === apply calving ===
896877
thickness(:) = thickness(:) - calvingThickness(:)
897878
879+
call update_calving_budget(domain)
898880
call remove_small_islands(meshPool, geometryPool)
899881
900882
block => block % next
@@ -968,6 +950,7 @@ subroutine floating_calving(domain, calvingFraction, err)
968950
! === apply calving ===
969951
thickness(:) = thickness(:) - calvingThickness(:)
970952
953+
call update_calving_budget(domain)
971954
call remove_small_islands(meshPool, geometryPool)
972955
973956
block => block % next
@@ -994,10 +977,14 @@ subroutine remove_small_islands(meshPool, geometryPool)
994977
995978
logical, pointer :: config_remove_small_islands
996979
real(kind=RKIND), pointer :: config_sea_level
997-
real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine)
980+
real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine)
981+
groundedCalvingThickness, &
982+
floatingCalvingThickness
998983
real (kind=RKIND), dimension(:), pointer :: thickness
999984
real (kind=RKIND), dimension(:), pointer :: bedTopography
1000-
integer, dimension(:), pointer :: cellMask
985+
integer, dimension(:), pointer :: cellMask, &
986+
groundedMaskForMassBudget, &
987+
floatingMaskForMassBudget
1001988
integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell
1002989
integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell
1003990
integer, pointer :: nCellsSolve
@@ -1014,6 +1001,10 @@ subroutine remove_small_islands(meshPool, geometryPool)
10141001
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
10151002
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
10161003
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
1004+
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
1005+
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
1006+
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
1007+
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
10171008
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
10181009
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
10191010
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
@@ -1037,6 +1028,12 @@ subroutine remove_small_islands(meshPool, geometryPool)
10371028
! If this is a single cell of ice surrounded by open ocean, kill this location
10381029
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell)
10391030
thickness(iCell) = 0.0_RKIND
1031+
! Update calving budgets
1032+
if (groundedMaskForMassBudget(iCell) .eq. 1) then
1033+
groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell)
1034+
elseif (floatingMaskForMassBudget(iCell) .eq. 1) then
1035+
floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell)
1036+
endif
10401037
elseif (nIceNeighbors == 1) then
10411038
! check if this neighbor has any additional neighbors with ice
10421039
nIceNeighbors2 = 0
@@ -1055,8 +1052,21 @@ subroutine remove_small_islands(meshPool, geometryPool)
10551052
! kill both cells
10561053
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell)
10571054
thickness(iCell) = 0.0_RKIND
1055+
! Update calving budgets
1056+
if (groundedMaskForMassBudget(iCell) .eq. 1) then
1057+
groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell)
1058+
elseif (floatingMaskForMassBudget(iCell) .eq. 1) then
1059+
floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell)
1060+
endif
1061+
10581062
calvingThickness(neighborWithIce) = calvingThickness(neighborWithIce) + thickness(neighborWithIce)
10591063
thickness(neighborWithIce) = 0.0_RKIND
1064+
! Update calving budgets
1065+
if (groundedMaskForMassBudget(neighborWithIce) .eq. 1) then
1066+
groundedCalvingThickness(neighborWithIce) = groundedCalvingThickness(neighborWithIce) + thickness(neighborWithIce)
1067+
elseif (floatingMaskForMassBudget(neighborWithIce) .eq. 1) then
1068+
floatingCalvingThickness(neighborWithIce) = floatingCalvingThickness(neighborWithIce) + thickness(neighborWithIce)
1069+
endif
10601070
endif
10611071
10621072
endif ! check on nIceNeighbors
@@ -1137,6 +1147,7 @@ subroutine topographic_calving(domain, calvingFraction, err)
11371147
! === apply calving ===
11381148
thickness(:) = thickness(:) - calvingThickness(:)
11391149
1150+
call update_calving_budget(domain)
11401151
call remove_small_islands(meshPool, geometryPool)
11411152
11421153
block => block % next
@@ -1316,6 +1327,7 @@ subroutine eigencalving(domain, err)
13161327
call mpas_timer_stop("halo updates")
13171328
! === apply calving ===
13181329
thickness(:) = thickness(:) - calvingThickness(:)
1330+
call update_calving_budget(domain)
13191331

13201332
! update mask
13211333
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
@@ -1340,7 +1352,7 @@ subroutine eigencalving(domain, err)
13401352
endif
13411353
enddo
13421354
! TODO: global reduce & reporting on amount of calving generated in this step
1343-
1355+
call update_calving_budget(domain)
13441356
! update mask
13451357
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
13461358
err = ior(err, err_tmp)
@@ -1365,7 +1377,7 @@ subroutine eigencalving(domain, err)
13651377
endif
13661378
enddo
13671379
! TODO: global reduce & reporting on amount of calving generated in this step
1368-
1380+
call update_calving_budget(domain)
13691381
call remove_small_islands(meshPool, geometryPool)
13701382
13711383
block => block % next
@@ -1494,7 +1506,7 @@ subroutine specified_calving_velocity(domain, err)
14941506
14951507
! === apply calving ===
14961508
thickness(:) = thickness(:) - calvingThickness(:)
1497-
1509+
call update_calving_budget(domain)
14981510
! update mask
14991511
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
15001512
err = ior(err, err_tmp)
@@ -1518,7 +1530,7 @@ subroutine specified_calving_velocity(domain, err)
15181530
endif
15191531
enddo
15201532
! TODO: global reduce & reporting on amount of calving generated in this step
1521-
1533+
call update_calving_budget(domain)
15221534
! update mask
15231535
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
15241536
err = ior(err, err_tmp)
@@ -1543,7 +1555,7 @@ subroutine specified_calving_velocity(domain, err)
15431555
endif
15441556
enddo
15451557
! TODO: global reduce & reporting on amount of calving generated in this step
1546-
1558+
call update_calving_budget(domain)
15471559
call remove_small_islands(meshPool, geometryPool)
15481560

15491561
block => block % next
@@ -1821,13 +1833,13 @@ subroutine von_Mises_calving(domain, err)
18211833

18221834
! === apply calving ===
18231835
thickness(:) = thickness(:) - calvingThickness(:)
1824-
1836+
call update_calving_budget(domain)
18251837
! update mask
18261838
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
18271839
err = ior(err, err_tmp)
18281840

18291841
call remove_small_islands(meshPool, geometryPool)
1830-
1842+
18311843
block => block % next
18321844

18331845
enddo ! associated(block)
@@ -2086,7 +2098,7 @@ subroutine ismip6_retreat(domain, err)
20862098

20872099
! === apply calving ===
20882100
thickness(:) = thickness(:) - calvingThickness(:)
2089-
2101+
call update_calving_budget(domain)
20902102
! update mask
20912103
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
20922104
err = ior(err, err_tmp)
@@ -2999,7 +3011,7 @@ subroutine damage_calving(domain, err)
29993011

30003012
! === apply calving ===
30013013
thickness(:) = thickness(:) - calvingThickness(:)
3002-
3014+
call update_calving_budget(domain)
30033015
! update mask
30043016
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
30053017
err = ior(err, err_tmp)
@@ -3015,7 +3027,7 @@ subroutine damage_calving(domain, err)
30153027
thickness(iCell) = 0.0_RKIND
30163028
endif
30173029
enddo
3018-
3030+
call update_calving_budget(domain)
30193031
! update mask
30203032
call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp)
30213033
err = ior(err, err_tmp)
@@ -3040,9 +3052,8 @@ subroutine damage_calving(domain, err)
30403052
endif
30413053
enddo
30423054
! TODO: global reduce & reporting on amount of calving generated in this step
3043-
3055+
call update_calving_budget(domain)
30443056
call remove_small_islands(meshPool, geometryPool)
3045-
30463057
block => block % next
30473058
30483059
enddo
@@ -3790,9 +3801,13 @@ subroutine remove_icebergs(domain)
37903801
integer, dimension(:), pointer :: seedMask, growMask !masks to pass to flood-fill routine
37913802
!integer, dimension(:), allocatable :: seedMaskOld !mask of where icebergs are removed, for debugging
37923803
3793-
real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine)
3804+
real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine)
3805+
groundedCalvingThickness, &
3806+
floatingCalvingThickness
37943807
real (kind=RKIND), dimension(:), pointer :: thickness
3795-
integer, dimension(:), pointer :: cellMask
3808+
integer, dimension(:), pointer :: cellMask, &
3809+
groundedMaskForMassBudget, &
3810+
floatingMaskForMassBudget
37963811
integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell
37973812
integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell
37983813
@@ -3877,6 +3892,10 @@ subroutine remove_icebergs(domain)
38773892
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
38783893
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
38793894
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
3895+
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
3896+
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
3897+
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
3898+
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
38803899
call mpas_pool_get_dimension(geometryPool, 'nCells', nCells)
38813900
call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve)
38823901
!allocate(seedMaskOld(nCells+1)) ! debug: make this a mask of where icebergs were removed
@@ -3887,6 +3906,11 @@ subroutine remove_icebergs(domain)
38873906
if (seedMask(iCell) == 0 .and. li_mask_is_floating_ice(cellMask(iCell))) then
38883907
calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) ! remove any remaining ice here
38893908
thickness(iCell) = 0.0_RKIND
3909+
if (groundedMaskForMassBudget(iCell) .eq. 1) then
3910+
groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell)
3911+
elseif (floatingMaskForMassBudget(iCell) .eq. 1) then
3912+
floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell)
3913+
endif
38903914
localIcebergCellCount = localIcebergCellCount + 1
38913915
!seedMaskOld(iCell) = 1 ! debug: make this a mask of where icebergs were removed
38923916
endif
@@ -4063,6 +4087,50 @@ subroutine li_flood_fill(seedMask, growMask, domain)
40634087
40644088
end subroutine li_flood_fill
40654089
4090+
4091+
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
4092+
!
4093+
! routine update_calving_budget
4094+
!
4095+
!> \brief Keep a running total of calving applied to grounded and floating
4096+
!> ice, respectively.
4097+
!> \author Trevor Hillebrand
4098+
!> \date May 2022
4099+
!> \details This routine should be called after each time time calvingThickness
4100+
!> is applied, but before masks are updated which often happens multiple times
4101+
!> in a timestep.
4102+
!-----------------------------------------------------------------------
4103+
subroutine update_calving_budget(domain)
4104+
!-----------------------------------------------------------------
4105+
! input/output variables
4106+
!-----------------------------------------------------------------
4107+
type (domain_type), intent(inout) :: domain
4108+
4109+
!-----------------------------------------------------------------
4110+
! local variables
4111+
!-----------------------------------------------------------------
4112+
type (mpas_pool_type), pointer :: meshPool, geometryPool
4113+
integer, dimension(:), pointer :: groundedMaskForMassBudget, & ! binary masks for mass budget
4114+
floatingMaskForMassBudget
4115+
real (kind=RKIND), dimension(:), pointer :: calvingThickness, &
4116+
groundedCalvingThickness, & ! Grounded and floating components for mass budget
4117+
floatingCalvingThickness
4118+
4119+
call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
4120+
call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget)
4121+
call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget)
4122+
call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness)
4123+
call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness)
4124+
call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
4125+
4126+
where (groundedMaskForMassBudget .eq. 1)
4127+
groundedCalvingThickness = calvingThickness
4128+
elsewhere (floatingMaskForMassBudget .eq. 1)
4129+
floatingCalvingThickness = calvingThickness
4130+
end where
4131+
4132+
end subroutine update_calving_budget
4133+
40664134
end module li_calving
40674135
40684136

0 commit comments

Comments
 (0)