diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index d41539915214..7523249ffe27 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1222,6 +1222,12 @@ is the value of that variable from the *previous* time level! + + @@ -1239,6 +1245,9 @@ is the value of that variable from the *previous* time level! + @@ -1266,12 +1275,21 @@ is the value of that variable from the *previous* time level! + + + @@ -1354,7 +1372,6 @@ is the value of that variable from the *previous* time level! - /> @@ -1431,6 +1448,12 @@ is the value of that variable from the *previous* time level! + + @@ -1453,8 +1476,12 @@ is the value of that variable from the *previous* time level! units="m" description="amount of change in bedTopography. This field is calculated by the sea-level model over its time step. Positive values indicate bed uplift relative to sea level." /> - + + diff --git a/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml b/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml index 50ad79158354..d16ed3eb8f43 100755 --- a/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml +++ b/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml @@ -85,9 +85,21 @@ + + + + diff --git a/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F b/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F index 10350fa9d235..76999b107fe5 100755 --- a/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F +++ b/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F @@ -174,11 +174,16 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) real (kind=RKIND), dimension(:), pointer :: bedTopography real (kind=RKIND), dimension(:), pointer :: sfcMassBalApplied real (kind=RKIND), dimension(:), pointer :: groundedSfcMassBalApplied + real (kind=RKIND), dimension(:), pointer :: floatingSfcMassBalApplied real (kind=RKIND), dimension(:), pointer :: basalMassBalApplied real (kind=RKIND), dimension(:), pointer :: groundedBasalMassBalApplied real (kind=RKIND), dimension(:), pointer :: floatingBasalMassBalApplied real (kind=RKIND), dimension(:), pointer :: calvingThickness - real (kind=RKIND), dimension(:), pointer :: faceMeltingThickness + real (kind=RKIND), dimension(:), pointer :: groundedCalvingThickness + real (kind=RKIND), dimension(:), pointer :: floatingCalvingThickness + real (kind=RKIND), dimension(:), pointer :: faceMeltingThickness, & + groundedFaceMeltingThickness, & + floatingFaceMeltingThickness real (kind=RKIND), dimension(:), pointer :: surfaceSpeed real (kind=RKIND), dimension(:), pointer :: basalSpeed real (kind=RKIND), dimension(:), pointer :: fluxAcrossGroundingLine @@ -195,6 +200,8 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) integer, dimension(:), pointer :: cellMask integer, dimension(:), pointer :: hydroMarineMarginMask integer, dimension(:), pointer :: hydroTerrestrialMarginMask + integer, dimension(:), pointer :: groundedMaskForMassBudget, & + floatingMaskForMassBudget integer, pointer :: nCellsSolve integer, pointer :: nEdgesSolve @@ -216,8 +223,11 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) real (kind=RKIND), pointer :: iceThicknessMax, iceThicknessMin, iceThicknessMean real (kind=RKIND), pointer :: totalSfcMassBal, totalBasalMassBal real (kind=RKIND), pointer :: totalGroundedSfcMassBal, totalFloatingSfcMassBal - real (kind=RKIND), pointer :: totalFaceMeltingFlux real (kind=RKIND), pointer :: totalGroundedBasalMassBal, totalFloatingBasalMassBal + real (kind=RKIND), pointer :: totalGroundedCalvingFlux, totalFloatingCalvingFlux + real (kind=RKIND), pointer :: totalFaceMeltingFlux, & + totalGroundedFaceMeltingFlux, & + totalFloatingFaceMeltingFlux real (kind=RKIND), pointer :: avgNetAccumulation real (kind=RKIND), pointer :: avgGroundedBasalMelt real (kind=RKIND), pointer :: avgSubshelfMelt @@ -247,8 +257,9 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) real (kind=RKIND) :: blockSumSfcMassBal, blockSumBasalMassBal real (kind=RKIND) :: blockSumGroundedSfcMassBal, blockSumFloatingSfcMassBal real (kind=RKIND) :: blockSumGroundedBasalMassBal, blockSumFloatingBasalMassBal - real (kind=RKIND) :: blockSumCalvingFlux - real (kind=RKIND) :: blockSumFaceMeltingFlux + real (kind=RKIND) :: blockSumCalvingFlux, blockSumGroundedCalvingFlux, blockSumFloatingCalvingFlux + real (kind=RKIND) :: blockSumFaceMeltingFlux, blockSumGroundedFaceMeltingFlux, & + blockSumFloatingFaceMeltingFlux real (kind=RKIND) :: blockMaxSurfaceSpeed real (kind=RKIND) :: blockMaxBasalSpeed real (kind=RKIND) :: blockGLflux @@ -292,7 +303,11 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) blockSumGroundedBasalMassBal = 0.0_RKIND blockSumFloatingBasalMassBal = 0.0_RKIND blockSumCalvingFlux = 0.0_RKIND + blockSumGroundedCalvingFlux = 0.0_RKIND + blockSumFloatingCalvingFlux = 0.0_RKIND blockSumFaceMeltingFlux = 0.0_RKIND + blockSumGroundedFaceMeltingFlux = 0.0_RKIND + blockSumFloatingFaceMeltingFlux = 0.0_RKIND blockGLflux = 0.0_RKIND blockGLMigrationFlux = 0.0_RKIND blockSumSubglacialWaterVolume = 0.0_RKIND @@ -343,12 +358,19 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'floatingSfcMassBalApplied', floatingSfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'faceMeltingThickness', faceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) + call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness) + call mpas_pool_get_array(geometryPool, 'groundedFaceMeltingThickness', groundedFaceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'floatingFaceMeltingThickness', floatingFaceMeltingThickness) 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_array(velocityPool, 'surfaceSpeed', surfaceSpeed) call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) @@ -372,20 +394,18 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) * areaCell(iCell) blockSumIceVolume = blockSumIceVolume + real(li_mask_is_ice_int(cellMask(iCell)),RKIND) & * areaCell(iCell) * thickness(iCell) - + ! can we somehow incorporate non-dynamic floating fringe into VAF? blockSumVAF = blockSumVAF + real(li_mask_is_grounded_ice_int(cellMask(iCell)),RKIND) * areaCell(iCell) * & ( thickness(iCell) + (rhow / rhoi) * min(0.0_RKIND, (bedTopography(iCell) - config_sea_level)) ) - blockSumGroundedIceArea = blockSumGroundedIceArea + real(li_mask_is_grounded_ice_int(cellMask(iCell)),RKIND) & + blockSumGroundedIceArea = blockSumGroundedIceArea + real(groundedMaskForMassBudget(iCell),RKIND) & * areaCell(iCell) - - blockSumGroundedIceVolume = blockSumGroundedIceVolume + real(li_mask_is_grounded_ice_int(cellMask(iCell)),RKIND) & + blockSumGroundedIceVolume = blockSumGroundedIceVolume + real(groundedMaskForMassBudget(iCell),RKIND) & * areaCell(iCell) * thickness(iCell) - blockSumFloatingIceArea = blockSumFloatingIceArea + real(li_mask_is_floating_ice_int(cellMask(iCell)),RKIND) & + blockSumFloatingIceArea = blockSumFloatingIceArea + real(floatingMaskForMassBudget(iCell),RKIND) & * areaCell(iCell) - - blockSumFloatingIceVolume = blockSumFloatingIceVolume + real(li_mask_is_floating_ice_int(cellMask(iCell)),RKIND) & + blockSumFloatingIceVolume = blockSumFloatingIceVolume + real(floatingMaskForMassBudget(iCell),RKIND) & * areaCell(iCell) * thickness(iCell) ! max, min thickness values (m) @@ -399,21 +419,27 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) ! SMB (kg yr^{-1}) blockSumSfcMassBal = blockSumSfcMassBal + areaCell(iCell) * sfcMassBalApplied(iCell) * scyr blockSumGroundedSfcMassBal = blockSumGroundedSfcMassBal + areaCell(iCell) * groundedSfcMassBalApplied(iCell) * scyr - blockSumFloatingSfcMassBal = blockSumFloatingSfcMassBal + & - (sfcMassBalApplied(iCell) - groundedSfcMassBalApplied(iCell)) * areaCell(iCell) * scyr - + blockSumFloatingSfcMassBal = blockSumFloatingSfcMassBal + areaCell(iCell) * floatingSfcMassBalApplied(iCell) * scyr ! BMB (kg yr-1) blockSumBasalMassBal = blockSumBasalMassBal + areaCell(iCell) * basalMassBalApplied(iCell) * scyr blockSumGroundedBasalMassBal = blockSumGroundedBasalMassBal + areaCell(iCell) * groundedBasalMassBalApplied(iCell) * scyr blockSumFloatingBasalMassBal = blockSumFloatingBasalMassBal + areaCell(iCell) * floatingBasalMassBalApplied(iCell) * scyr - ! mass loss due do calving (kg yr^{-1}) + ! mass loss due to calving (kg yr^{-1}) blockSumCalvingFlux = blockSumCalvingFlux + calvingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) + blockSumGroundedCalvingFlux = blockSumGroundedCalvingFlux + groundedCalvingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) + blockSumFloatingCalvingFlux = blockSumFloatingCalvingFlux + floatingCalvingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) ! mass loss due to face-melting (kg yr^{-1}) blockSumFaceMeltingFlux = blockSumFaceMeltingFlux + faceMeltingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) + blockSumGroundedFaceMeltingFlux = blockSumGroundedFaceMeltingFlux + groundedFaceMeltingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) + blockSumFloatingFaceMeltingFlux = blockSumFloatingFaceMeltingFlux + floatingFaceMeltingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) ! max surface speed if (surfaceSpeed(iCell) > blockMaxSurfaceSpeed) then @@ -531,25 +557,29 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) sums(11) = blockSumGroundedBasalMassBal sums(12) = blockSumFloatingBasalMassBal sums(13) = blockSumCalvingFlux - sums(14) = blockSumFaceMeltingFlux - sums(15) = blockSumVAF - sums(16) = blockGLflux - sums(17) = blockGLMigrationflux + sums(14) = blockSumGroundedCalvingFlux + sums(15) = blockSumFloatingCalvingFlux + sums(16) = blockSumFaceMeltingFlux + sums(17) = blockSumGroundedFaceMeltingFlux + sums(18) = blockSumFloatingFaceMeltingFlux + sums(19) = blockSumVAF + sums(20) = blockGLflux + sums(21) = blockGLMigrationflux if (config_SGH) then - sums(18) = blockSumSubglacialWaterVolume - sums(19) = blockSumBasalMeltInput - sums(20) = blockSumExternalWaterInput - sums(21) = blockSumChannelMelt - sums(22) = blockSumLakeVolume - sums(23) = blockSumLakeArea - sums(24) = blockSumGLMeltFlux - sums(25) = blockSumTerrestrialMeltFlux - sums(26) = blockSumChannelGLMeltFlux - sums(27) = blockSumChannelTerrestrialMeltFlux - sums(28) = blockSumFlotationFraction - nVars = 28 + sums(22) = blockSumSubglacialWaterVolume + sums(23) = blockSumBasalMeltInput + sums(24) = blockSumExternalWaterInput + sums(25) = blockSumChannelMelt + sums(26) = blockSumLakeVolume + sums(27) = blockSumLakeArea + sums(28) = blockSumGLMeltFlux + sums(29) = blockSumTerrestrialMeltFlux + sums(30) = blockSumChannelGLMeltFlux + sums(31) = blockSumChannelTerrestrialMeltFlux + sums(32) = blockSumFlotationFraction + nVars = 32 else - nVars = 17 + nVars = 21 endif call mpas_dmpar_sum_real_array(dminfo, nVars, sums(1:nVars), reductions(1:nVars)) @@ -576,7 +606,11 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) call mpas_pool_get_array(globalStatsAMPool, 'totalFloatingBasalMassBal', totalFloatingBasalMassBal) call mpas_pool_get_array(globalStatsAMPool, 'avgSubshelfMelt', avgSubshelfMelt) call mpas_pool_get_array(globalStatsAMPool, 'totalCalvingFlux', totalCalvingFlux) + call mpas_pool_get_array(globalStatsAMPool, 'totalGroundedCalvingFlux', totalGroundedCalvingFlux) + call mpas_pool_get_array(globalStatsAMPool, 'totalFloatingCalvingFlux', totalFloatingCalvingFlux) call mpas_pool_get_array(globalStatsAMPool, 'totalFaceMeltingFlux', totalFaceMeltingFlux) + call mpas_pool_get_array(globalStatsAMPool, 'totalGroundedFaceMeltingFlux', totalGroundedFaceMeltingFlux) + call mpas_pool_get_array(globalStatsAMPool, 'totalFloatingFaceMeltingFlux', totalFloatingFaceMeltingFlux) call mpas_pool_get_array(globalStatsAMPool, 'groundingLineFlux', groundingLineFlux) call mpas_pool_get_array(globalStatsAMPool, 'groundingLineMigrationFlux', groundingLineMigrationFlux) if (config_SGH) then @@ -606,22 +640,26 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) totalGroundedBasalMassBal = reductions(11) totalFloatingBasalMassBal = reductions(12) totalCalvingFlux = reductions(13) - totalFaceMeltingFlux = reductions(14) - volumeAboveFloatation = reductions(15) - groundingLineFlux = reductions(16) - groundingLineMigrationFlux = reductions(17) + totalGroundedCalvingFlux = reductions(14) + totalFloatingCalvingFlux = reductions(15) + totalFaceMeltingFlux = reductions(16) + totalGroundedFaceMeltingFlux = reductions(17) + totalFloatingFaceMeltingFlux = reductions(18) + volumeAboveFloatation = reductions(19) + groundingLineFlux = reductions(20) + groundingLineMigrationFlux = reductions(21) if (config_SGH) then - totalSubglacialWaterVolume = reductions(18) - totalBasalMeltInput = reductions(19) - totalExternalWaterInput = reductions(20) - totalChannelMelt = reductions(21) - totalSubglacialLakeVolume = reductions(22) - totalSubglacialLakeArea = reductions(23) - totalDistWaterFluxMarineMargin = reductions(24) - totalDistWaterFluxTerrestrialMargin = reductions(25) - totalChnlWaterFluxMarineMargin = reductions(26) - totalChnlWaterFluxTerrestrialMargin = reductions(27) - totalFlotationFraction = reductions(28) + totalSubglacialWaterVolume = reductions(22) + totalBasalMeltInput = reductions(23) + totalExternalWaterInput = reductions(24) + totalChannelMelt = reductions(25) + totalSubglacialLakeVolume = reductions(26) + totalSubglacialLakeArea = reductions(27) + totalDistWaterFluxMarineMargin = reductions(28) + totalDistWaterFluxTerrestrialMargin = reductions(29) + totalChnlWaterFluxMarineMargin = reductions(30) + totalChnlWaterFluxTerrestrialMargin = reductions(31) + totalFlotationFraction = reductions(32) endif if (totalIceArea > 0.0_RKIND) then 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 4cb0f3e132db..a91b843c0f27 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 @@ -46,8 +46,7 @@ module li_advection !-------------------------------------------------------------------- public :: li_advection_thickness_tracers, & li_layer_normal_velocity, & - li_update_geometry, & - li_grounded_to_floating + li_update_geometry !-------------------------------------------------------------------- ! @@ -114,7 +113,6 @@ subroutine li_advection_thickness_tracers(& type (mpas_pool_type), intent(inout) :: & velocityPool !< Input/output: velocity information ! (needs to be inout for li_calculate_mask call - type (mpas_pool_type), intent(inout) :: & geometryPool !< Input/output: geometry information to be updated @@ -147,6 +145,7 @@ subroutine li_advection_thickness_tracers(& sfcMassBal, & ! surface mass balance (potential forcing) sfcMassBalApplied, & ! surface mass balance (actually applied) groundedSfcMassBalApplied, & ! surface mass balance on grounded locations (actually applied) + floatingSfcMassBalApplied, & ! surface mass balance on floating locations (actually applied) basalMassBal, & ! basal mass balance groundedBasalMassBal, & ! basal mass balance for grounded ice floatingBasalMassBal, & ! basal mass balance for floating ice @@ -195,7 +194,9 @@ subroutine li_advection_thickness_tracers(& cellMaskTemporaryField, & ! scratch field containing old values of cellMask edgeMaskTemporaryField, & vertexMaskTemporaryField, & - thermalCellMaskField + thermalCellMaskField, & + groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget + floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget ! Allocatable arrays need for flux-corrected transport advection real (kind=RKIND), dimension(:,:,:), allocatable :: tend @@ -204,7 +205,9 @@ subroutine li_advection_thickness_tracers(& cellMask, & ! integer bitmask for cells edgeMask, & ! integer bitmask for edges vertexMask, & - thermalCellMask + thermalCellMask, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget integer, dimension(:,:), pointer :: cellsOnEdge @@ -241,7 +244,7 @@ subroutine li_advection_thickness_tracers(& !WHL - debug integer, dimension(:), pointer :: indexToCellID, indexToEdgeID - integer, pointer :: nCells, nEdges + integer, pointer :: nCells, nEdges, nEdgesSolve integer, pointer :: config_stats_cell_ID integer :: iEdge, iEdgeOnCell integer, dimension(:), pointer :: nEdgesOnCell @@ -283,6 +286,7 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'floatingSfcMassBalApplied', floatingSfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'basalMassBal', basalMassBal) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBal', groundedBasalMassBal) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) @@ -296,7 +300,8 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening) call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) - + call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) ! get arrays from the velocity pool call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity) call mpas_pool_get_array(velocityPool, 'layerThicknessEdgeFlux', layerThicknessEdgeFlux) @@ -323,6 +328,7 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) + call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'maxEdges2', maxEdges2) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) @@ -357,6 +363,11 @@ 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.) + call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField) call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.) @@ -376,15 +387,14 @@ subroutine li_advection_thickness_tracers(& ! given the old thickness, compute the thickness in each layer call li_calculate_layerThickness(meshPool, thickness, layerThickness) - ! Save copies of masks because we need to preserve mask - ! states prior to advection for accurate time integration. - ! A mask update is necessary to calculate grounding line flux, - ! after which we will reset the masks to their previous states. cellMaskTemporaryField % array(:) = cellMask(:) edgeMaskTemporaryField % array(:) = edgeMask(:) vertexMaskTemporaryField % array(:) = vertexMask(:) - - layerThicknessEdgeFlux(:,:) = 0.0_RKIND + ! Update halo on edgeMask for calculating flux across grounding line. + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'edgeMask') + call mpas_dmpar_field_halo_exch(domain, 'cellMask') + call mpas_timer_stop("halo updates") !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers @@ -457,7 +467,7 @@ subroutine li_advection_thickness_tracers(& ! copy old (input) values of layer thickness and tracers to scratch arrays layerThicknessOld(:,:) = layerThickness(:,:) advectedTracersOld(:,:,:) = advectedTracers(:,:,:) - + layerThicknessEdgeFlux(:,:) = 0.0_RKIND ! compute new values of layer thickness and tracers if ((trim(config_thickness_advection) .eq. 'fo') .and. & ( (trim(config_tracer_advection) .eq. 'fo') .or. & @@ -547,6 +557,52 @@ subroutine li_advection_thickness_tracers(& enddo endif + ! Update halos after advection. + call mpas_timer_start("halo updates") + 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 @@ -557,8 +613,6 @@ subroutine li_advection_thickness_tracers(& ! before advection took place. thickness = sum(layerThickness, 1) - - !----------------------------------------------------------------- ! Add the surface and basal mass balance to the layer thickness !----------------------------------------------------------------- @@ -568,18 +622,22 @@ subroutine li_advection_thickness_tracers(& ! TODO: more complicated treatment at GL? - do iCell = 1, nCells - if (li_mask_is_grounded_ice(cellMask(iCell))) then - basalMassBal(iCell) = groundedBasalMassBal(iCell) - elseif (li_mask_is_floating_ice(cellMask(iCell))) then - ! Currently, floating and grounded ice are mutually exclusive. - ! This could change if the GL is parameterized, in which case this logic may need adjustment. - basalMassBal(iCell) = floatingBasalMassBal(iCell) - elseif ( .not. (li_mask_is_ice(cellMask(iCell)))) then - ! We don't allow a positive basal mass balance where ice is not already present. - basalMassBal(iCell) = 0.0_RKIND - endif - enddo + where ( groundedMaskForMassBudget .eq. 1 ) + + basalMassBal = groundedBasalMassBal + + elsewhere ( floatingMaskForMassBudget .eq. 1) + + ! Currently, floating and grounded ice are mutually exclusive. + ! This could change if the GL is parameterized, in which case this logic may need adjustment. + basalMassBal = floatingBasalMassBal + + elsewhere ( .not. (li_mask_is_ice(cellMask) ) ) + + ! We don't allow a positive basal mass balance where ice is not already present. + basalMassBal = 0.0_RKIND + + end where ! It is possible that excess internal melting was computed and assigned ! to the drainedInternalMeltRate array in mpas_li_thermal.F. If so, then add it to basalMassBal. @@ -607,10 +665,13 @@ subroutine li_advection_thickness_tracers(& sfcMassBal, & sfcMassBalApplied, & groundedSfcMassBalApplied, & + floatingSfcMassBalApplied, & basalMassBal, & basalMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & surfaceTracers, & basalTracers, & layerThickness, & @@ -642,52 +703,9 @@ subroutine li_advection_thickness_tracers(& ! We need an updated set of masks to calculate fluxAcrossGroundingLine, ! but we will reset this to the previous state below for accuracy of the ! time integration scheme. - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! 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, nEdges - 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 - - ! Reset masks to state before advection and mass balance for - ! accuracy of time integration scheme. - cellMask(:) = cellMaskTemporaryField % array(:) - edgeMask(:) = edgeMaskTemporaryField % array(:) - vertexMask(:) = vertexMaskTemporaryField % array(:) - ! 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. @@ -860,10 +878,13 @@ subroutine apply_mass_balance(& sfcMassBal, & sfcMassBalApplied, & groundedSfcMassBalApplied, & + floatingSfcMassBalApplied, & basalMassBal, & basalMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & surfaceTracers, & basalTracers, & layerThickness, & @@ -877,7 +898,9 @@ subroutine apply_mass_balance(& dt, & !< Input: time step (s) rhoi !< Input: ice density (kg/m^3) - integer, dimension(:), intent(in) :: cellMask !< Input: mask on cells + integer, dimension(:), intent(in) :: cellMask, & !< Input: mask on cells + groundedMaskForMassBudget, & + floatingMaskForMassBudget real (kind=RKIND), dimension(:), intent(in) :: & bedTopography, & !< Input: bed elevation (m) @@ -906,7 +929,8 @@ subroutine apply_mass_balance(& sfcMassBalApplied !< Output: surface mass balance actually applied on this time step (kg/m^2/s) real(kind=RKIND), dimension(:), intent(out) :: & - groundedSfcMassBalApplied !< Output: surface mass balance actually applied to grounded ice on this time step (kg/m^2/s) + groundedSfcMassBalApplied, & !< Output: surface mass balance actually applied to grounded ice on this time step (kg/m^2/s) + floatingSfcMassBalApplied !< Output: surface mass balance actually applied to floating ice on this time step (kg/m^2/s) real(kind=RKIND), dimension(:), intent(out) :: & basalMassBalApplied, & !< Output: basal mass balance actually applied on this time step (kg/m^2/s) @@ -1062,25 +1086,23 @@ subroutine apply_mass_balance(& enddo ! iCell ! Separate grounded and floating components, as necessary. - where (li_mask_is_grounded_ice(cellMask) .or. bedTopography > config_sea_level) + where ( (groundedMaskForMassBudget .eq. 1) .or. (bedTopography > config_sea_level) ) groundedSfcMassBalApplied = sfcMassBalApplied + floatingSfcMassBalApplied = 0.0_RKIND + groundedBasalMassBalApplied = basalMassBalApplied + floatingBasalMassBalApplied = 0.0_RKIND + elsewhere (floatingMaskForMassBudget .eq. 1) + groundedSfcMassBalApplied = 0.0_RKIND + floatingSfcMassBalApplied = sfcMassBalApplied + groundedBasalMassBalApplied = 0.0_RKIND + floatingBasalMassBalApplied = basalMassBalApplied elsewhere groundedSfcMassBalApplied = 0.0_RKIND + floatingSfcMassBalApplied = 0.0_RKIND + groundedBasalMassBalApplied = 0.0_RKIND + floatingBasalMassBalApplied = 0.0_RKIND end where - do iCell = 1, nCells - if (li_mask_is_grounded_ice(cellMask(iCell))) then - groundedBasalMassBalApplied(iCell) = basalMassBalApplied(iCell) - floatingBasalMassBalApplied(iCell) = 0.0_RKIND - elseif (li_mask_is_floating_ice(cellMask(iCell))) then - floatingBasalMassBalApplied(iCell) = basalMassBalApplied(iCell) - groundedBasalMassBalApplied(iCell) = 0.0_RKIND - else - groundedBasalMassBalApplied(iCell) = 0.0_RKIND - floatingBasalMassBalApplied(iCell) = 0.0_RKIND - endif - enddo - deallocate(thckTracerProducts) end subroutine apply_mass_balance @@ -2112,52 +2134,6 @@ subroutine vertical_remap(thickness, meshPool, layerThickness, tracers, err) end subroutine vertical_remap - subroutine li_grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells) - !----------------------------------------------------------------- - ! input variables - !----------------------------------------------------------------- - integer, dimension(:), intent(in) :: & - cellMaskOrig !< Input: mask for cells before advection - - integer, dimension(:), intent(in) :: & - cellMaskNew !< Input: mask for cells after advection - - 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 (li_mask_is_grounded_ice(cellMaskOrig(iCell)) .and. & - li_mask_is_floating_ice(cellMaskNew(iCell)) ) then - groundedToFloatingThickness(iCell) = thicknessNew(iCell) - ! find locations that had floating ice before but now have grounded ice - else if (li_mask_is_floating_ice(cellMaskOrig(iCell)) .and. & - li_mask_is_grounded_ice(cellMaskNew(iCell)) ) then - groundedToFloatingThickness(iCell) = -1.0_RKIND * thicknessNew(iCell) - endif - enddo - - end subroutine li_grounded_to_floating - !*********************************************************************** end module li_advection diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_bedtopo.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_bedtopo.F index 51c9a8f69e23..be5ea3281e8c 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_bedtopo.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_bedtopo.F @@ -251,7 +251,7 @@ subroutine li_bedtopo_solve(domain, err) bedTopography(:) = bedTopography(:) + upliftRate(:) * deltat call li_update_geometry(geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) block => block % next end do @@ -761,7 +761,7 @@ subroutine slmodel_solve(slmTimeStep, domain) ! Perform Halo exchange update call mpas_dmpar_field_halo_exch(domain,'bedTopography') call li_update_geometry(geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! deallocate memory diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 1637c7d697df..ba5073f2c424 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -129,9 +129,13 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) bedTopography ! bed topography (negative below sea level) real (kind=RKIND), dimension(:), pointer :: & - calvingThickness ! thickness of ice that calves (computed in this subroutine) - ! typically the entire ice thickness, but will be a fraction of the thickness - ! if calving_timescale > dt + calvingThickness, & ! thickness of ice that calves (computed in this subroutine) + ! typically the entire ice thickness, but will be a fraction of the thickness + ! if calving_timescale > dt + groundedCalvingThickness, & ! Grounded and floating components for mass budget + floatingCalvingThickness + integer, dimension(:), pointer :: groundedMaskForMassBudget, & ! binary masks for mass budget + floatingMaskForMassBudget real (kind=RKIND), dimension(:), pointer :: & calvingThicknessFromThreshold ! thickness of ice that calves from threshold processes (computed in this subroutine) @@ -200,7 +204,7 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) ! Update mask and geometry before calling any calving routines. May not be necessary, but best be safe. call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call li_update_geometry(geometryPool) @@ -243,6 +247,24 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) end do endif + ! Zero calvingThickness here instead of or in addition to in individual subroutines. + ! This is necessary when using damage threshold calving with other calving + ! routines. Some individual routines still set calvingThickness to zero, but + ! this is redundant. li_apply_front_ablation_velocity will zero + ! calvingThickness, so any routines that use that should not be applied + ! after other calving routines. + block => domain % blocklist + 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 ! compute calvingThickness based on the calving_config option if (trim(config_calving) == 'none') then @@ -310,7 +332,7 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) ! Consider mask calving as a possible additional step ! Mask calving can occur by itself or in conjunction with a physical calving law @@ -319,16 +341,16 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) err = ior(err, err_tmp) endif - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call remove_small_islands(domain, err_tmp) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) ! now also remove any icebergs call li_remove_icebergs(domain) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) ! Final operations after calving has been applied, including removal ! of small islands @@ -365,7 +387,7 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) endif ! config_print_calving_info ! Update mask and geometry - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call li_update_geometry(geometryPool) @@ -471,6 +493,8 @@ subroutine li_restore_calving_front(domain, err) ! > 0 for cells below sea level that were initially ice-free and now have ice restoreThickness ! thickness of ice that is added to restore the calving front to its initial position ! > 0 for cells below sea level that were initially ice-covered and now have very thin or no ice + real (kind=RKIND), dimension(:), pointer :: & + incrementalCalvingThickness ! incremental calving thickness real (kind=RKIND) :: & restoreThicknessMin ! small thickness to which ice is restored should it fall below this thickness @@ -524,6 +548,7 @@ subroutine li_restore_calving_front(domain, err) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'restoreThickness', restoreThickness) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'calvingVelocity', calvingVelocity) @@ -575,7 +600,7 @@ subroutine li_restore_calving_front(domain, err) restoreThicknessMin = 0.1_RKIND * config_dynamic_thickness ! calculate masks - so we know where the calving front was located initially - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! initialize @@ -616,6 +641,9 @@ subroutine li_restore_calving_front(domain, err) endif ! if preventing retreat ! Now check for marine regions that have advanced + incrementalCalvingThickness(:) = 0.0_RKIND + + ! loop over locally owned cells do iCell = 1, nCells if (bedTopography(iCell) < config_sea_level) then ! The bed is below sea level; test for calving-front advance and retreat. @@ -628,7 +656,7 @@ subroutine li_restore_calving_front(domain, err) call mpas_log_write('Remove ice: indexToCellID=$i, thickness=$r', intArgs=(/indexToCellID(iCell)/), realArgs=(/thickness(iCell)/)) endif - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND ! Assign flow speed to calving speed - this allows calculation of licalvf in ISMIP6 postprocessing calvingVelocity(iCell) = surfaceSpeed(iCell) @@ -636,6 +664,8 @@ subroutine li_restore_calving_front(domain, err) endif ! bedTopography < config_sea_level enddo ! iCell + call update_calving_budget(geometryPool, incrementalCalvingThickness) + block => block % next enddo @@ -646,7 +676,7 @@ subroutine li_restore_calving_front(domain, err) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call li_update_geometry(geometryPool) block => block % next @@ -678,7 +708,7 @@ subroutine li_restore_calving_front(domain, err) if (err > 0) then call mpas_log_write("An error has occurred in li_restore_calving_front.", MPAS_LOG_ERR) endif - + call mpas_log_write("finished with li_restore_calving_front") end subroutine li_restore_calving_front @@ -779,6 +809,8 @@ subroutine thickness_calving(domain, calvingFraction, err) real (kind=RKIND), dimension(:), pointer :: & calvingThicknessFromThreshold ! thickness of ice that calves (computed in this subroutine) + real (kind=RKIND), dimension(:), pointer :: & + incrementalCalvingThickness integer, pointer :: nCells integer :: iCell, iCellOnCell, iCellNeighbor @@ -846,7 +878,7 @@ subroutine thickness_calving(domain, calvingFraction, err) ! Initialize calvingThickness = 0.0_RKIND - + incrementalCalvingThickness(:) = 0.0_RKIND ! Identify cells that are active-for-calving: ! (1) Grounded ice with thickness > config_dynamic_thickness ! (2) Floating ice with thickness > config_calving_thickness @@ -936,12 +968,14 @@ subroutine thickness_calving(domain, calvingFraction, err) ! Calve ice in ocean cells that are not on the protected inactive margin where (oceanMask == 1 .and. inactiveMarginMask == 0 .and. thickness > 0.0_RKIND) - calvingThickness = thickness * calvingFraction + incrementalCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next enddo @@ -989,8 +1023,11 @@ subroutine floating_calving(domain, calvingFraction, err) type (mpas_pool_type), pointer :: meshPool real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine) real (kind=RKIND), dimension(:), pointer :: calvingThicknessFromThreshold ! thickness of ice that calves (computed in this subroutine) + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask + integer, dimension(:), pointer :: floatingMaskForMassBudget + integer, pointer :: nCells err = 0 @@ -1005,17 +1042,21 @@ subroutine floating_calving(domain, calvingFraction, err) call mpas_pool_get_array(geometryPool, 'calvingThicknessFromThreshold', calvingThicknessFromThreshold) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) - - calvingThickness = 0.0_RKIND - - ! Note: The floating_ice mask includes all floating ice, both inactive and active - where (li_mask_is_floating_ice(cellMask)) - calvingThickness = thickness * calvingFraction + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + call mpas_pool_get_array(meshPool, 'nCells', nCells) + + incrementalCalvingThickness(:) = 0.0_RKIND + ! Note: The floating_ice mask does not include floating + ! non-dynamic cells adjacent to grounded cells. + where ( floatingMaskForMassBudget .eq. 1 ) + incrementalCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next enddo @@ -1042,11 +1083,17 @@ subroutine remove_small_islands(domain, err) integer, intent(inout) :: err type (mpas_pool_type), pointer :: scratchPool, meshPool, geometryPool, velocityPool logical, pointer :: config_remove_small_islands - real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine) + real(kind=RKIND), pointer :: config_sea_level real (kind=RKIND), dimension(:), pointer :: calvingThicknessFromThreshold ! 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 :: incrementalCalvingThickness 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 :: nCells, maxEdges @@ -1075,13 +1122,18 @@ subroutine remove_small_islands(domain, err) 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, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingThicknessFromThreshold', calvingThicknessFromThreshold) + 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) allocate(connectedCellsList((maxEdges+1)**2), & islandMask(nCells+1)) - + incrementalCalvingThickness(:) = 0.0_RKIND islandMask(:) = 0 nIslandCellsLocal = 0 nIslandCellsGlobal = 0 @@ -1217,7 +1269,7 @@ subroutine remove_small_islands(domain, err) call mpas_dmpar_field_halo_exch(domain, 'thickness') call mpas_dmpar_field_halo_exch(domain, 'calvingThickness') call mpas_timer_stop("halo updates") - call li_calculate_mask(meshPool, velocityPool, geometryPool, err) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) call mpas_dmpar_sum_int(domain % dminfo, nIslandCellsLocal, nIslandCellsGlobal) @@ -1236,8 +1288,8 @@ subroutine remove_small_islands(domain, err) call li_flood_fill(seedMask, growMask, domain) do iCell = 1, nCells if (li_mask_is_floating_ice(cellMask(iCell)) .and. seedMask(iCell) == 0) then - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo @@ -1248,6 +1300,8 @@ subroutine remove_small_islands(domain, err) call mpas_deallocate_scratch_field(seedMaskField, single_block_in=.true.) call mpas_deallocate_scratch_field(growMaskField, single_block_in=.true.) + call update_calving_budget(geometryPool, incrementalCalvingThickness) + end subroutine remove_small_islands @@ -1285,11 +1339,13 @@ subroutine topographic_calving(domain, calvingFraction, err) type (mpas_pool_type), pointer :: geometryPool, meshPool real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine) real (kind=RKIND), dimension(:), pointer :: calvingThicknessFromThreshold + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real(kind=RKIND), pointer :: config_calving_topography real(kind=RKIND), pointer :: config_sea_level logical, pointer :: config_print_calving_info real (kind=RKIND), dimension(:), pointer :: bedTopography, thickness integer, dimension(:), pointer :: cellMask + integer, pointer :: nCells err = 0 @@ -1309,20 +1365,24 @@ subroutine topographic_calving(domain, calvingFraction, err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingThicknessFromThreshold', calvingThicknessFromThreshold) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(meshPool, 'nCells', nCells) - calvingThickness = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND where ( (li_mask_is_floating_ice(cellMask)) .and. (bedTopography < config_calving_topography + config_sea_level) ) - calvingThickness = thickness * calvingFraction + incrementalCalvingThickness = thickness * calvingFraction endwhere - calvingThicknessFromThreshold = calvingThickness + calvingThicknessFromThreshold = incrementalCalvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next enddo @@ -1381,6 +1441,7 @@ subroutine eigencalving(domain, err) real (kind=RKIND), dimension(:), pointer :: angleEdge real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: calvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell @@ -1438,6 +1499,7 @@ subroutine eigencalving(domain, err) call mpas_pool_get_array(velocityPool, 'eMin', eMin) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingCFLdt', calvingCFLdt) call mpas_pool_get_array(geometryPool, 'dtCalvingCFLratio', dtCalvingCFLratio) call mpas_pool_get_array(geometryPool, 'totalRatebasedCalvedVolume', totalRatebasedCalvedVolume) @@ -1462,6 +1524,7 @@ subroutine eigencalving(domain, err) realArgs=(/minval(eigencalvingParameter), maxval(eigencalvingParameter)/)) endif + incrementalCalvingThickness(:) = 0.0_RKIND calvingVelocity(:) = 0.0_RKIND ! First calculate the front retreat rate (Levermann eq. 1) calvingVelocity(:) = eigencalvingParameter(:) * max(0.0_RKIND, eMax(:)) * max(0.0_RKIND, eMin(:)) ! m/s @@ -1469,7 +1532,7 @@ subroutine eigencalving(domain, err) call mpas_log_write("calling li_apply_front_ablation_velocity from eigencalving") ! Convert calvingVelocity to calvingThickness call li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool, & - calvingThickness, calvingVelocity, applyToGrounded, & + incrementalCalvingThickness, calvingVelocity, applyToGrounded, & applyToFloating, applyToGroundingLine, domain, calvingCFLdt, dtCalvingCFLratio, & totalRatebasedCalvedVolume, totalRatebasedUncalvedVolume, err_tmp) err = ior(err, err_tmp) @@ -1481,7 +1544,7 @@ subroutine eigencalving(domain, err) call mpas_log_write('config_damage_calving_method == threshold; & removing ice with damage > $r', realArgs=(/config_damage_calving_threshold/)) - call apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, err_tmp) + call apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, incrementalCalvingThickness, err_tmp) err = ior(err, err_tmp) elseif ( trim(config_damage_calving_method) == 'calving_rate' ) then call mpas_log_write('config_damage_calving_method == calving_rate & @@ -1500,13 +1563,14 @@ subroutine eigencalving(domain, err) ! why not, so leaving it. ! NOTE: THIS WILL NOT WORK ON MULTIPLE BLOCKS PER PROCESSOR call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'calvingThickness') + call mpas_dmpar_field_halo_exch(domain, 'incrementalCalvingThickness') call mpas_timer_stop("halo updates") ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. @@ -1523,36 +1587,17 @@ subroutine eigencalving(domain, err) ! Tests of the current implementation show reasonable behavior. do iCell = 1, nCells if (calvingFrontMask(iCell) == 1 .and. thickness(iCell) < config_calving_thickness) then - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step - + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! remove abandoned floating ice (i.e. icebergs) and add it to the calving flux - ! Defined as: floating ice (dynamic or non-dynamic) that is not adjacent to dynamic ice (floating or grounded) - ! This won't necessarily find all abandoned ice, but in practice it does a pretty good job at general cleanup - calvingSubtotal = 0.0_RKIND - do iCell = 1, nCells - if (li_mask_is_floating_ice(cellMask(iCell))) then - ! check neighbors for dynamic ice (floating or grounded) - dynamicNeighbor = .false. - do iNeighbor = 1, nEdgesOnCell(iCell) - jCell = cellsOnCell(iNeighbor, iCell) - if (li_mask_is_dynamic_ice(cellMask(jCell))) dynamicNeighbor = .true. - enddo - if (.not. dynamicNeighbor) then ! calve this ice - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) - thickness(iCell) = 0.0_RKIND - calvingSubtotal = calvingSubtotal + calvingThickness(iCell) * areaCell(iCell) - endif - endif - enddo - ! TODO: global reduce & reporting on amount of calving generated in this step + call li_remove_icebergs(domain) block => block % next enddo @@ -1605,6 +1650,7 @@ subroutine specified_calving_velocity(domain, err) real (kind=RKIND), dimension(:), pointer :: angleEdge real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: calvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell @@ -1650,6 +1696,7 @@ subroutine specified_calving_velocity(domain, err) call mpas_pool_get_array(geometryPool, 'calvingVelocityData', calvingVelocityData) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingCFLdt', calvingCFLdt) call mpas_pool_get_array(geometryPool, 'dtCalvingCFLratio', dtCalvingCFLratio) call mpas_pool_get_array(geometryPool, 'totalRatebasedCalvedVolume', totalRatebasedCalvedVolume) @@ -1658,6 +1705,7 @@ subroutine specified_calving_velocity(domain, err) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) + incrementalCalvingThickness(:) = 0.0_RKIND ! get parameter value if (trim(config_calving_specified_source) == 'const') then @@ -1676,19 +1724,21 @@ subroutine specified_calving_velocity(domain, err) endif ! Convert calvingVelocity to calvingThickness - call li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool, calvingThickness, calvingVelocity, & + call li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool, incrementalCalvingThickness, calvingVelocity, & applyToGrounded=.true., applyToFloating=.true., applyToGroundingLine=.false., & domain=domain, maxDt=calvingCFLdt, CFLratio=dtCalvingCFLratio, & totalAblatedVolume=totalRatebasedCalvedVolume, & totalUnablatedVolume=totalRatebasedUncalvedVolume, & err=err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'incrementalCalvingThickness') + call mpas_timer_stop("halo updates") ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. @@ -1705,38 +1755,18 @@ subroutine specified_calving_velocity(domain, err) ! Tests of the current implementation show reasonable behavior. do iCell = 1, nCells if (calvingFrontMask(iCell) == 1 .and. thickness(iCell) < config_calving_thickness) then - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step - + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! remove abandoned floating ice (i.e. icebergs) and add it to the calving flux - ! Defined as: floating ice (dynamic or non-dynamic) that is not adjacent to dynamic ice (floating or grounded) - ! This won't necessarily find all abandoned ice, but in practice it does a pretty good job at general cleanup - calvingSubtotal = 0.0_RKIND - do iCell = 1, nCells - if (li_mask_is_floating_ice(cellMask(iCell))) then - ! check neighbors for dynamic ice (floating or grounded) - dynamicNeighbor = .false. - do iNeighbor = 1, nEdgesOnCell(iCell) - jCell = cellsOnCell(iNeighbor, iCell) - if (li_mask_is_dynamic_ice(cellMask(jCell))) dynamicNeighbor = .true. - enddo - if (.not. dynamicNeighbor) then ! calve this ice - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) - thickness(iCell) = 0.0_RKIND - calvingSubtotal = calvingSubtotal + calvingThickness(iCell) * areaCell(iCell) - endif - endif - enddo - ! TODO: global reduce & reporting on amount of calving generated in this step - block => block % next + enddo end subroutine specified_calving_velocity @@ -1797,6 +1827,7 @@ subroutine von_Mises_calving(domain, err) floatingVonMisesThresholdStress, & groundedVonMisesThresholdStress, & surfaceSpeed + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:,:), pointer :: flowParamA, & temperature, layerThickness real (kind=RKIND), pointer :: config_default_flowParamA @@ -1861,6 +1892,7 @@ subroutine von_Mises_calving(domain, err) call mpas_pool_get_array(velocityPool, 'yvelmean', yvelmean) call mpas_pool_get_array(velocityPool, 'surfaceSpeed', surfaceSpeed) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingVelocity', calvingVelocity) call mpas_pool_get_array(geometryPool, 'calvingCFLdt', calvingCFLdt) call mpas_pool_get_array(geometryPool, 'dtCalvingCFLratio', dtCalvingCFLratio) @@ -1909,7 +1941,8 @@ subroutine von_Mises_calving(domain, err) endif vonMisesStress(:) = 0.0_RKIND - + incrementalCalvingThickness(:) = 0.0_RKIND + ! get flowParamA from MPAS or use Albany-like equation if ( config_use_Albany_flowA_eqn_for_vM ) then !calculate Albany-type flowParamA @@ -2045,7 +2078,7 @@ subroutine von_Mises_calving(domain, err) call mpas_log_write("calling li_apply_front_ablation_velocity from von Mises stress calving routine") ! Convert calvingVelocity to calvingThickness call li_apply_front_ablation_velocity(meshPool, geometryPool,velocityPool, & - calvingThickness, calvingVelocity, applyToGrounded, & + incrementalCalvingThickness, calvingVelocity, applyToGrounded, & applyToFloating, applyToGroundingLine, domain, & calvingCFLdt, dtCalvingCFLratio, & totalRatebasedCalvedVolume, totalRatebasedUncalvedVolume, err_tmp) @@ -2057,7 +2090,7 @@ subroutine von_Mises_calving(domain, err) ! why not, so leaving it. ! NOTE: THIS WILL NOT WORK ON MULTIPLE BLOCKS PER PROCESSOR call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'calvingThickness') + call mpas_dmpar_field_halo_exch(domain, 'incrementalCalvingThickness') call mpas_timer_stop("halo updates") ! If floating VM threshold is zero, remove floating ice. Remove @@ -2069,7 +2102,7 @@ subroutine von_Mises_calving(domain, err) do iCell = 1,nCells if ( li_mask_is_floating_ice(cellMask(iCell)) .and. & li_mask_is_dynamic_ice(cellMask(iCell)) ) then - calvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) elseif ( li_mask_is_floating_ice(cellMask(icell)) .and. & (.not. li_mask_is_dynamic_ice(cellMask(iCell))) ) then nGroundedNeighbors = 0 @@ -2080,12 +2113,19 @@ subroutine von_Mises_calving(domain, err) endif enddo if ( nGroundedNeighbors == 0 ) then - calvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) endif endif enddo endif + ! === apply calving velocity before damage threshold calving === + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) + ! update mask + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) + err = ior(err, err_tmp) + if ( trim(config_damage_calving_method) == 'none' ) then ! do nothing elseif ( trim(config_damage_calving_method) == 'threshold' ) then @@ -2093,7 +2133,7 @@ subroutine von_Mises_calving(domain, err) call mpas_log_write('config_damage_calving_method == threshold; & removing ice with damage > $r', realArgs=(/config_damage_calving_threshold/)) - call apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, err_tmp) + call apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, incrementalCalvingThickness, err_tmp) err = ior(err, err_tmp) elseif ( trim(config_damage_calving_method) == 'calving_rate' ) then call mpas_log_write('config_damage_calving_method == calving_rate & @@ -2106,23 +2146,15 @@ subroutine von_Mises_calving(domain, err) return endif - ! Update halos on calvingThickness or faceMeltingThickness before - ! applying it. - ! Testing seemed to indicate this is not necessary, but I don't - ! understand - ! why not, so leaving it. - ! NOTE: THIS WILL NOT WORK ON MULTIPLE BLOCKS PER PROCESSOR - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'calvingThickness') - call mpas_timer_stop("halo updates") - - ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - + ! === apply calving velocity before damage threshold calving === + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call li_remove_icebergs(domain) + block => block % next enddo ! associated(block) @@ -2189,6 +2221,7 @@ subroutine ismip6_retreat(domain, err) calvingVelocity, thickness, & xvelmean, yvelmean, calvingThickness, & bedTopography + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, pointer :: nCells, timestepNumber @@ -2232,6 +2265,7 @@ subroutine ismip6_retreat(domain, err) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingVelocity', calvingVelocity) call mpas_pool_get_array(geometryPool, 'calvingCFLdt', calvingCFLdt) call mpas_pool_get_array(geometryPool, 'dtCalvingCFLratio', dtCalvingCFLratio) @@ -2249,6 +2283,8 @@ subroutine ismip6_retreat(domain, err) ! submergedArea used in runoff unit conversion allocate(submergedArea(nCells+1)) + submergedArea(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND streamFound = .false. ! Changed to true if ismip6-gis stream is found, otherwise throws error stream_cursor => domain % streamManager % streams % head @@ -2372,7 +2408,7 @@ subroutine ismip6_retreat(domain, err) call mpas_log_write("calling li_apply_front_ablation_velocity from ismip6 retreat routine") ! Convert calvingVelocity to calvingThickness call li_apply_front_ablation_velocity(meshPool, geometryPool,velocityPool, & - calvingThickness, calvingVelocity, applyToGrounded=.true., & + incrementalCalvingThickness, calvingVelocity, applyToGrounded=.true., & applyToFloating=.true., applyToGroundingLine=.false., & domain=domain, maxDt=calvingCFLdt, CFLratio=dtCalvingCFLratio, & totalAblatedVolume=totalRatebasedCalvedVolume, & @@ -2383,18 +2419,17 @@ subroutine ismip6_retreat(domain, err) ! Update halos on calvingThickness before applying it. ! NOTE: THIS WILL NOT WORK ON MULTIPLE BLOCKS PER PROCESSOR call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'calvingThickness') + call mpas_dmpar_field_halo_exch(domain, 'incrementalCalvingThickness') call mpas_timer_stop("halo updates") ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) deallocate(submergedArea) - end subroutine ismip6_retreat !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -3327,6 +3362,7 @@ subroutine damage_calving(domain, err) real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: bedTopography real (kind=RKIND), dimension(:), pointer :: calvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:), pointer :: calvingVelocity real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio @@ -3379,6 +3415,7 @@ subroutine damage_calving(domain, err) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_array(meshPool, 'areaCell', areaCell) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) @@ -3390,6 +3427,7 @@ subroutine damage_calving(domain, err) call mpas_pool_get_array(geometryPool, 'totalRatebasedCalvedVolume', totalRatebasedCalvedVolume) call mpas_pool_get_array(geometryPool, 'totalRatebasedUncalvedVolume', totalRatebasedUncalvedVolume) + incrementalCalvingThickness(:) = 0.0_RKIND call calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask) @@ -3411,17 +3449,20 @@ subroutine damage_calving(domain, err) err=err_tmp) err = ior(err, err_tmp) elseif (trim(config_damage_calving_method) == 'threshold') then - call apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, err_tmp) + call apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, incrementalCalvingThickness, err_tmp) err = ior(err, err_tmp) else call mpas_log_write("Unknown value for config_damage_calving_method was specified!", MPAS_LOG_ERR) endif + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'incrementalCalvingThickness') + call mpas_timer_stop("halo updates") ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. @@ -3431,35 +3472,16 @@ subroutine damage_calving(domain, err) ! Tests of the current implementation show reasonable behavior. do iCell = 1, nCells if (calvingFrontMask(iCell) == 1 .and. thickness(iCell) < config_calving_thickness) then - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo - + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! remove abandoned floating ice (i.e. icebergs) and add it to the calving flux - ! Defined as: floating ice (dynamic or non-dynamic) that is not adjacent to dynamic ice (floating or grounded) - ! This won't necessarily find all abandoned ice, but in practice it does a pretty good job at general cleanup - calvingSubtotal = 0.0_RKIND - do iCell = 1, nCells - if (li_mask_is_floating_ice(cellMask(iCell))) then - ! check neighbors for dynamic ice (floating or grounded) - dynamicNeighbor = .false. - do iNeighbor = 1, nEdgesOnCell(iCell) - jCell = cellsOnCell(iNeighbor, iCell) - if (li_mask_is_dynamic_ice(cellMask(jCell))) dynamicNeighbor = .true. - enddo - if (.not. dynamicNeighbor) then ! calve this ice - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) - thickness(iCell) = 0.0_RKIND - calvingSubtotal = calvingSubtotal + calvingThickness(iCell) * areaCell(iCell) - endif - endif - enddo - ! TODO: global reduce & reporting on amount of calving generated in this step + call li_remove_icebergs(domain) block => block % next @@ -3934,7 +3956,7 @@ end subroutine li_finalize_damage_after_advection !> value exceeds a specified threshold, assuming the ice is connected to the calving !> front. !----------------------------------------------------------------------- - subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, err) + subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, incrementalCalvingThickness, err) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- @@ -3946,12 +3968,11 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object type (mpas_pool_type), pointer, intent(inout) :: geometryPool !< Input: geometry pool - !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - + real (kind=RKIND), dimension(:), intent(out) :: incrementalCalvingThickness ! incremental calving thickness !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -3992,6 +4013,8 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d growMask => growMaskField % array growMask(:) = 0 + incrementalCalvingThickness(:) = 0.0_RKIND + ! define seed and grow masks for flood fill. where ( li_mask_is_dynamic_margin(cellMask) .and. (damage .ge. config_damage_calving_threshold) ) seedMask = 1 @@ -4011,13 +4034,13 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d jCell = cellsOnCell(iNeighbor, iCell) if (li_mask_is_floating_ice(cellMask(jCell)) .and. .not. li_mask_is_dynamic_ice(cellMask(jCell))) then ! this is a thin neighbor - remove the whole cell volume - calvingThickness(jCell) = thickness(jCell) + incrementalCalvingThickness(jCell) = thickness(jCell) calvingThicknessFromThreshold(jCell) = calvingThicknessFromThreshold(jCell) + thickness(jCell) endif enddo - calvingThickness(iCell) = thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) endif ! if cell is calving margin enddo ! cell loop @@ -4072,6 +4095,7 @@ subroutine mask_calving(domain, err) integer, pointer :: nCells, nCellsSolve integer :: iCell integer :: localMaskCellCount, globalMaskCellCount + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness integer :: err_tmp err = 0 @@ -4089,13 +4113,17 @@ subroutine mask_calving(domain, err) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingThicknessFromThreshold', calvingThicknessFromThreshold) call mpas_pool_get_array(geometryPool, 'calvingMask', calvingMask) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(meshPool, 'nCells', nCells) + + incrementalCalvingThickness(:) = 0.0_RKIND ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) localMaskCellCount = 0 @@ -4104,7 +4132,7 @@ subroutine mask_calving(domain, err) do iCell = 1, nCells if (li_mask_is_floating_ice(cellMask(iCell)) .and. (calvingMask(iCell) >= 1)) then !call mpas_log_write("Found masked cell at $i", intArgs=(/iCell/)) - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) thickness(iCell) = 0.0_RKIND if (iCell <= nCellsSolve) localMaskCellCount = localMaskCellCount + 1 @@ -4114,8 +4142,9 @@ subroutine mask_calving(domain, err) call mpas_dmpar_sum_int(domain % dminfo, localMaskCellCount, globalMaskCellCount) call mpas_log_write("Mask calving applied. Removed $i floating cells with mask.", intArgs=(/globalMaskCellCount/)) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) block => block % next @@ -4228,10 +4257,15 @@ subroutine li_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 :: calvingThicknessFromThreshold + real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine) + groundedCalvingThickness, & + floatingCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness 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 @@ -4286,7 +4320,7 @@ subroutine li_remove_icebergs(domain) ! make sure masks are up to date. May not be necessary, but safer to ! do anyway. - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call mpas_log_write("Iceberg-detection flood-fill: updated masks.") @@ -4328,7 +4362,12 @@ subroutine li_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, 'incrementalCalvingThickness', incrementalCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingThicknessFromThreshold', calvingThicknessFromThreshold) + 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 @@ -4337,14 +4376,15 @@ subroutine li_remove_icebergs(domain) localIcebergCellCount = 0 do iCell = 1, nCellsSolve if (seedMask(iCell) == 0 .and. li_mask_is_floating_ice(cellMask(iCell))) then - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) ! remove any remaining ice here calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) ! remove any remaining ice here + incrementalCalvingThickness(iCell) = thickness(iCell) ! remove any remaining ice here thickness(iCell) = 0.0_RKIND localIcebergCellCount = localIcebergCellCount + 1 !seedMaskOld(iCell) = 1 ! debug: make this a mask of where icebergs were removed endif enddo + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next end do @@ -4362,6 +4402,7 @@ subroutine li_remove_icebergs(domain) call mpas_deallocate_scratch_field(growMaskField, single_block_in=.false.) call mpas_log_write("Iceberg-detection flood-fill complete. Removed $i iceberg cells.", intArgs=(/globalIcebergCellCount/)) call mpas_timer_stop("iceberg detection") + end subroutine li_remove_icebergs !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| @@ -4425,7 +4466,7 @@ subroutine li_flood_fill(seedMask, growMask, domain) ! make sure masks are up to date. May not be necessary, but safer to ! do anyway. - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) !call mpas_log_write("Flood-fill: updated masks.") @@ -4515,6 +4556,52 @@ 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(geometryPool, incrementalCalvingThickness) + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + type (mpas_pool_type), pointer, intent(inout) :: geometryPool + real (kind=RKIND), dimension(:), intent(inout) :: incrementalCalvingThickness + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + 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_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 = groundedCalvingThickness + incrementalCalvingThickness + elsewhere (floatingMaskForMassBudget .eq. 1) + floatingCalvingThickness = floatingCalvingThickness + incrementalCalvingThickness + end where + + calvingThickness(:) = calvingThickness(:) + incrementalCalvingThickness(:) + incrementalCalvingThickness(:) = 0.0_RKIND + + end subroutine update_calving_budget + end module li_calving diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index 0e6ff08df22d..3ae5e3c308ab 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -300,7 +300,7 @@ function li_core_init(domain, startTimeStamp) result(err) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call li_update_geometry(geometryPool) @@ -803,7 +803,7 @@ function li_core_initial_solve(domain) result(err) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call li_update_geometry(geometryPool) block => block % next diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F index daa1d5cca87d..00e381d275a6 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F @@ -106,8 +106,12 @@ subroutine li_face_melt_grounded_ice(domain, err) real (kind=RKIND), dimension(:), pointer :: & faceMeltSpeed, & faceMeltingThickness, & + groundedFaceMeltingThickness, & + floatingFaceMeltingThickness, & thickness - integer, dimension(:), pointer :: cellMask + integer, dimension(:), pointer :: cellMask, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget integer :: err_tmp logical :: applyToFloating, applyToGrounded, applyToGroundingLine @@ -141,6 +145,10 @@ subroutine li_face_melt_grounded_ice(domain, err) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'faceMeltingThickness', faceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'groundedFaceMeltingThickness', groundedFaceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'floatingFaceMeltingThickness', floatingFaceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) call mpas_pool_get_array(geometryPool, 'thickness', thickness) ! Update halos on calvingThickness or faceMeltingThickness before applying it. @@ -153,9 +161,19 @@ subroutine li_face_melt_grounded_ice(domain, err) ! Apply facemelt: open the Ark thickness(:) = thickness(:) - faceMeltingThickness(:) + where (groundedMaskForMassBudget .eq. 1) + groundedFaceMeltingThickness = faceMeltingThickness + floatingFaceMeltingThickness = 0.0_RKIND + elsewhere (floatingMaskForMassBudget .eq. 1) + groundedFaceMeltingThickness = 0.0_RKIND + floatingFaceMeltingThickness = faceMeltingThickness + elsewhere + groundedFaceMeltingThickness = 0.0_RKIND + floatingFaceMeltingThickness = 0.0_RKIND + end where ! update mask - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) block => block % next @@ -170,9 +188,13 @@ subroutine li_face_melt_grounded_ice(domain, err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_array(geometryPool, 'faceMeltSpeed', faceMeltSpeed) call mpas_pool_get_array(geometryPool, 'faceMeltingThickness', faceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'groundedFaceMeltingThickness', groundedFaceMeltingThickness) + call mpas_pool_get_array(geometryPool, 'floatingFaceMeltingThickness', floatingFaceMeltingThickness) faceMeltSpeed(:) = 0.0_RKIND faceMeltingThickness(:) = 0.0_RKIND + groundedFaceMeltingThickness(:) = 0.0_RKIND + floatingfaceMeltingThickness(:) = 0.0_RKIND block => block % next enddo ! associated(block) @@ -390,7 +412,7 @@ subroutine li_basal_melt_floating_ice(domain, err) thermalCellMask => thermalCellMaskField % array ! calculate masks - so we know where the ice is floating - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! calculate a mask to identify ice that is thick enough to be thermally active diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index bf9616b4339b..ce4f0d727d42 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -172,7 +172,7 @@ subroutine li_SGH_init(domain, err) endif ! Mask needs to be initialized for pressure calcs to be correct - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call calc_hydro_mask(domain) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F index acc5048e3a4e..7fec0b821932 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F @@ -858,7 +858,7 @@ subroutine li_thermal_solver(domain, err) endif ! calculate masks - so we know where the ice is floating - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! calculate a mask to identify ice that is thick enough to be thermally active diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F index e6051a3ec753..7e98bb4c3819 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F @@ -87,7 +87,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) use li_velocity use li_bedtopo use li_mask - use li_advection, only: li_grounded_to_floating + use li_advection use li_ocean_extrap !----------------------------------------------------------------- @@ -244,19 +244,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advancing clock") -!TODO: Determine whether grounded melting should in fact be called first ! === Ocean forcing extrapolation into ice-shelf cavities =========== call mpas_timer_start("ocean forcing extrapolation") call li_ocean_extrap_solve(domain, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("ocean forcing extrapolation") -! === Face melting for grounded ice =========== - call mpas_timer_start("face melting for grounded ice") - call li_face_melt_grounded_ice(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("face melting for grounded ice") - ! === Basal melting for floating ice =========== call mpas_timer_start("basal melting for floating ice") call li_basal_melt_floating_ice(domain, err_tmp) @@ -373,7 +366,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Calculate masks prior to RK loop, but do not update masks within the loop ! to preserve the accuracy of time integration. - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! *** Start RK loop *** @@ -521,7 +514,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Finalize budget updates ! Update masks after RK integration - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) @@ -533,10 +526,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") - ! 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 li_grounded_to_floating(cellMaskPrev, cellMask, thickness, groundedToFloatingThickness, nCells) - sfcMassBalApplied(:) = sfcMassBalAccum(:) groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) basalMassBalApplied(:) = basalMassBalAccum(:) @@ -594,6 +583,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") +! === Face melting for grounded ice =========== + call mpas_timer_start("face melting for grounded ice") + call li_face_melt_grounded_ice(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("face melting for grounded ice") + ! === Update bed topo ===================== ! It's not clear when the best time to do this is. ! Seems cleaner to do it either before or after all of the time evolution of the ice @@ -697,6 +692,7 @@ subroutine prepare_advection(domain, err) real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity real (kind=RKIND), pointer :: calvingCFLdt, faceMeltingCFLdt integer, pointer :: processLimitingTimestep, config_rk_order, config_rk3_stages + real (kind=RKIND), dimension(:), pointer :: groundedToFloatingThickness integer, dimension(:), pointer :: edgeMask logical, pointer :: config_print_thickness_advection_info @@ -790,8 +786,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( & @@ -1259,7 +1260,7 @@ subroutine advection_solver(domain, err) call li_calculate_layerThickness(meshPool, thickness, layerThickness) endif - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call li_update_geometry(geometryPool) block => block % next diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F index cf80b2497ebe..acffc35328ae 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F @@ -333,7 +333,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call li_update_geometry(geometryPool) block => block % next diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F index 077ff4c1773f..29f7a6fab809 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F @@ -875,7 +875,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! We could call diagnostic_solve_before_velocity to do all that, but the way ! code is currently organized, that would be a circular dependency.a - call li_calculate_mask(meshPool, velocityPool, geometryPool, err) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) ! Lower surface is based on floatation for floating ice. For grounded ice (and non-ice areas) it is the bed. where ( li_mask_is_floating_ice(cellMask) ) diff --git a/components/mpas-albany-landice/src/shared/mpas_li_mask.F b/components/mpas-albany-landice/src/shared/mpas_li_mask.F index 64050cb83adf..3fd0b7215213 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -258,7 +258,7 @@ end subroutine li_calculate_mask_init ! !----------------------------------------------------------------------- - subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) + subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) !----------------------------------------------------------------- ! @@ -277,6 +277,8 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) ! input/output variables ! !----------------------------------------------------------------- + type (domain_type), intent(inout) :: domain + type (mpas_pool_type), intent(inout) :: & geometryPool !< Input/Output: geometry information @@ -295,8 +297,9 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) !----------------------------------------------------------------- integer, pointer :: nCells, nVertices, nEdges, vertexDegree integer, pointer :: nVertInterfaces - real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography - integer, dimension(:), pointer :: nEdgesOnCell, cellMask, vertexMask, edgeMask + real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography, groundedToFloatingThickness + integer, dimension(:), pointer :: nEdgesOnCell, cellMask, vertexMask, edgeMask, & + groundedMaskForMassBudget, floatingMaskForMassBudget integer, dimension(:,:), pointer :: cellsOnCell, cellsOnVertex, cellsOnEdge, dirichletVelocityMask real (kind=RKIND), pointer :: config_ice_density, config_ocean_density, & config_sea_level, config_dynamic_thickness @@ -318,6 +321,10 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) real (kind=RKIND) :: thinnestNeighborHeight 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') @@ -340,6 +347,14 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, 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) @@ -366,12 +381,24 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, 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. + ! Initialize grounded and floating masks for mass budget as well. + floatingMaskForMassBudget(:) = 0 + groundedMaskForMassBudget(:) = 0 where ( li_mask_is_ice(cellMask) .and. (config_ice_density / config_ocean_density * thickness) <= & (config_sea_level - bedTopography) ) cellMask = ior(cellMask, li_mask_ValueFloating) + floatingMaskForMassBudget = 1 + elsewhere ( li_mask_is_ice(cellMask) .and. (config_ice_density / config_ocean_density * thickness) > & + (config_sea_level - bedTopography) ) + groundedMaskForMassBudget = 1 end where ! Identify the margin @@ -495,6 +522,64 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) enddo enddo + ! For mass budget calculations, add floating ice over subglacial lakes and + ! floating non-dynamic cells bordering grounded ice to the grounded mask + ! and remove them from floating mask. Need the following if statement to + ! avoid circularity between li_flood_fill and this routine. + if (updateMassBudgetMasks) then + do i=1,nCells + if (li_mask_is_floating_ice(cellMask(i))) then + if (.not. li_mask_is_dynamic_ice(cellMask(i))) then + do j=1,nEdgesOnCell(i) ! Check if any neighbors are grounded + iCellNeighbor = cellsOnCell(j,i) + if (li_mask_is_grounded_ice(cellMask(iCellNeighbor))) then + groundedMaskForMassBudget(i) = 1 + floatingMaskForMassBudget(i) = 0 + cycle ! no need to look at additional neighbors + endif + enddo + endif + + ! Seed mask with floating cells with open ocean neighbors. + ! Currently commented because flood fill is not yet available + ! here. + ! do j=1,nEdgesOnCell(i) ! Check if any neighbors are open ocean + ! iCellNeighbor = cellsOnCell(j,i) + ! if ( (.not. li_mask_is_ice(cellMask(iCellNeighbor))) .and. & + ! (bedTopography(iCellNeighbor) < config_sea_level) ) then + ! seedMask = 1 + ! cycle ! no need to look at additional neighbors + ! endif + ! enddo + 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 + ! because flood fill needs to be moved from li_calving to a shared + ! routine. + ! call li_flood_fill(seedMask=seedMask, growMask=li_mask_is_floating_ice(cellMask), domain) + + !where (li_mask_is_floating(cellMask) and seedMask .eq. 0) + ! groundedMaskForBudget = 1 + ! gloatingMaskForMassBudget = 0 + !end where + + endif + !call mpas_timer_stop('calculate mask cell') ! ====