From 8d1b472cfac6313c6e5330581bc59bcd4be4b6ff Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 13:58:39 -0600 Subject: [PATCH 01/32] Add specialized masks for mass budget calculations Add grounded and floating masks that are used for the mass budget calculation and are different from the masks in cellMask. Currently, the only difference is that the grounded mask includes floating non-dynamic cells adjacent to grounded ice. The grounded masks will eventually also include floating ice over subglacial lakes, but that requires flood_fill to be moved from li_calving to a shared routine, so that will come in a future commit. In preparation for that, add domain variable in calls to li_calculate_mask. --- .../mpas-albany-landice/src/Registry.xml | 6 ++ .../src/mode_forward/mpas_li_advection.F | 10 ++- .../src/mode_forward/mpas_li_bedtopo.F | 4 +- .../src/mode_forward/mpas_li_calving.F | 28 ++++----- .../src/mode_forward/mpas_li_core.F | 2 +- .../src/mode_forward/mpas_li_iceshelf_melt.F | 4 +- .../mode_forward/mpas_li_subglacial_hydro.F | 2 +- .../src/mode_forward/mpas_li_thermal.F | 2 +- .../mpas_li_time_integration_fe_rk.F | 4 +- .../src/mode_forward/mpas_li_velocity.F | 2 +- .../mode_forward/mpas_li_velocity_external.F | 2 +- .../src/shared/mpas_li_mask.F | 61 ++++++++++++++++++- 12 files changed, 95 insertions(+), 32 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index d41539915214..69309dbe205c 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! + + 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..e468e6a31777 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 @@ -82,6 +82,7 @@ subroutine li_advection_thickness_tracers(& geometryPool, & thermalPool, & scratchPool, & + domain, & err, & advectTracersIn) @@ -111,10 +112,11 @@ subroutine li_advection_thickness_tracers(& !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object + type (domain_type), intent(inout) :: domain + 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 @@ -376,10 +378,6 @@ 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(:) @@ -642,7 +640,7 @@ 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 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..841b9677ba15 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 @@ -200,7 +200,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) @@ -575,7 +575,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 @@ -646,7 +646,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 @@ -1506,7 +1506,7 @@ subroutine eigencalving(domain, err) thickness(:) = thickness(:) - calvingThickness(:) ! 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. @@ -1530,7 +1530,7 @@ subroutine eigencalving(domain, err) ! TODO: global reduce & reporting on amount of calving generated in this step ! 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 @@ -1688,7 +1688,7 @@ subroutine specified_calving_velocity(domain, err) thickness(:) = thickness(:) - calvingThickness(:) ! 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. @@ -1712,7 +1712,7 @@ subroutine specified_calving_velocity(domain, err) ! TODO: global reduce & reporting on amount of calving generated in this step ! 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 @@ -2120,7 +2120,7 @@ subroutine von_Mises_calving(domain, err) thickness(:) = thickness(:) - calvingThickness(:) ! 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 @@ -2390,7 +2390,7 @@ subroutine ismip6_retreat(domain, err) thickness(:) = thickness(:) - calvingThickness(:) ! 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) @@ -3421,7 +3421,7 @@ subroutine damage_calving(domain, err) thickness(:) = thickness(:) - calvingThickness(:) ! 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. @@ -3437,7 +3437,7 @@ subroutine damage_calving(domain, err) enddo ! 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 @@ -4115,7 +4115,7 @@ subroutine mask_calving(domain, err) call mpas_log_write("Mask calving applied. Removed $i floating cells with mask.", intArgs=(/globalMaskCellCount/)) ! 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 @@ -4286,7 +4286,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.") @@ -4425,7 +4425,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.") 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..b74ea4efb18a 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 @@ -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..fedac9bdb1b3 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 @@ -155,7 +155,7 @@ subroutine li_face_melt_grounded_ice(domain, err) thickness(:) = thickness(:) - faceMeltingThickness(:) ! 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 @@ -390,7 +390,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..3c429eb8f762 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 @@ -1161,6 +1161,7 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & + domain, & err_tmp, & advectTracersIn = .true.) @@ -1180,6 +1181,7 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & + domain, & err_tmp, & advectTracersIn = .false.) @@ -1259,7 +1261,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..fbb1bb1fad98 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 @@ -296,7 +298,8 @@ 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 + 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,7 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) real (kind=RKIND) :: thinnestNeighborHeight integer :: iCellNeighbor logical :: openOceanNeighbor + logical :: updateMassBudgetMasks=.true. call mpas_timer_start('calculate mask') @@ -340,6 +344,8 @@ 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, 'groundedMaskForMassBudget', groundedMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) @@ -369,9 +375,14 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) ! 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. 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 +506,52 @@ 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 + endif + cycle ! no need to look at additional neighbors + 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 + ! endif + ! cycle ! no need to look at additional neighbors + ! enddo + 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') ! ==== From 6d25a87aaa4687d9dfdaa931d9494e7120c26eb4 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 14:47:52 -0600 Subject: [PATCH 02/32] Add floatingSfcMassBalApplied variable Also redefine applied surface and basal mass balance fields based on new grounded and floating masks for mass budgets. --- .../mpas-albany-landice/src/Registry.xml | 3 ++ .../src/mode_forward/mpas_li_advection.F | 29 +++++++++++++++---- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 69309dbe205c..d5e3036fed3b 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1245,6 +1245,9 @@ is the value of that variable from the *previous* time level! + 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 e468e6a31777..e1e1b03dc6f1 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 @@ -149,6 +149,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 @@ -206,7 +207,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 @@ -285,6 +288,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) @@ -298,7 +302,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) @@ -605,10 +610,13 @@ subroutine li_advection_thickness_tracers(& sfcMassBal, & sfcMassBalApplied, & groundedSfcMassBalApplied, & + floatingSfcMassBalApplied, & basalMassBal, & basalMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & surfaceTracers, & basalTracers, & layerThickness, & @@ -858,10 +866,13 @@ subroutine apply_mass_balance(& sfcMassBal, & sfcMassBalApplied, & groundedSfcMassBalApplied, & + floatingSfcMassBalApplied, & basalMassBal, & basalMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & surfaceTracers, & basalTracers, & layerThickness, & @@ -875,7 +886,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) @@ -904,7 +917,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) @@ -1060,10 +1074,15 @@ 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 + elsewhere (floatingMaskForMassBudget .eq. 1) + groundedSfcMassBalApplied = 0.0_RKIND + floatingSfcMassBalApplied = sfcMassBalApplied elsewhere groundedSfcMassBalApplied = 0.0_RKIND + floatingSfcMassBalApplied = 0.0_RKIND end where do iCell = 1, nCells From 8391ee0f985c70244b6c8bc2269ae084a68fec1e Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 15:04:34 -0600 Subject: [PATCH 03/32] Add grounded and floating calvingThickness variables Add variables to distinguish between grounded and floating components of calving in mass budgets. --- .../mpas-albany-landice/src/Registry.xml | 5 ++++ .../src/mode_forward/mpas_li_calving.F | 26 ++++++++++++++++--- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index d5e3036fed3b..9516c2d8fb77 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1280,6 +1280,11 @@ is the value of that variable from the *previous* time level! /> + 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) @@ -339,6 +343,8 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) + call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) ! In data calving mode we just calculate what should be calved but don't actually calve it. @@ -368,6 +374,20 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) call li_update_geometry(geometryPool) + ! necessary to get these after masks have been updated + call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + + groundedCalvingThickness(:) = 0.0_RKIND + floatingCalvingThickness(:) = 0.0_RKIND + where (groundedMaskForMassBudget .eq. 1) + groundedCalvingThickness = calvingThickness + elsewhere (floatingMaskForMassBudget .eq. 1) + floatingCalvingThickness = calvingThickness + elsewhere + groundedCalvingThickness = 0.0_RKIND + floatingCalvingThickness = 0.0_RKIND + end where block => block % next end do From 7d11453e8c28d80de3d87bce0acade7177434cd9 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 15:26:31 -0600 Subject: [PATCH 04/32] Add grounded and floating faceMeltingThickness components Add grounded and floating faceMeltingThickness components for calculating mass budgets. Also fix implementation of parsing grounded and floating calvingThickness components, which needs to be done after calvingThickness is applied but before the masks are updated. --- .../mpas-albany-landice/src/Registry.xml | 5 +++ .../src/mode_forward/mpas_li_calving.F | 37 ++++++++++++------- .../src/mode_forward/mpas_li_iceshelf_melt.F | 24 +++++++++++- 3 files changed, 51 insertions(+), 15 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 9516c2d8fb77..fc122bdf13cb 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1444,6 +1444,11 @@ is the value of that variable from the *previous* time level! /> + domain % blocklist + do while (associated(block)) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) + call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness) + call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + + groundedCalvingThickness(:) = 0.0_RKIND + floatingCalvingThickness(:) = 0.0_RKIND + where (groundedMaskForMassBudget .eq. 1) + groundedCalvingThickness = calvingThickness + elsewhere (floatingMaskForMassBudget .eq. 1) + floatingCalvingThickness = calvingThickness + elsewhere + groundedCalvingThickness = 0.0_RKIND + floatingCalvingThickness = 0.0_RKIND + end where + + block => block % next + end do + ! Final operations after calving has been applied, including removal ! of small islands block => domain % blocklist @@ -374,20 +397,6 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) call li_update_geometry(geometryPool) - ! necessary to get these after masks have been updated - call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) - call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) - - groundedCalvingThickness(:) = 0.0_RKIND - floatingCalvingThickness(:) = 0.0_RKIND - where (groundedMaskForMassBudget .eq. 1) - groundedCalvingThickness = calvingThickness - elsewhere (floatingMaskForMassBudget .eq. 1) - floatingCalvingThickness = calvingThickness - elsewhere - groundedCalvingThickness = 0.0_RKIND - floatingCalvingThickness = 0.0_RKIND - end where block => block % next end do 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 fedac9bdb1b3..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,6 +161,16 @@ 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, domain, err_tmp) @@ -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) From ed5f652b68981e7a6762a2e7c2066d0b40ec215b Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 16:08:39 -0600 Subject: [PATCH 05/32] Update global stats calculations using new masks Also account for grounded and floating components of calving and face melting. It's still not clear what to do about VAF. --- .../Registry_global_stats.xml | 12 ++ .../analysis_members/mpas_li_global_stats.F | 135 ++++++++++++------ 2 files changed, 101 insertions(+), 46 deletions(-) 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..19cfa96ae7ba 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 @@ -218,6 +225,10 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) 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 +258,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 +304,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 +359,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 +395,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,17 +420,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 calving (kg yr^{-1}) + blockSumFaceMeltingFlux = blockSumFaceMeltingFlux + faceMeltingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) + blockSumGroundedFaceMeltingFlux = blockSumFaceMeltingFlux + groundedFaceMeltingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) + blockSumFloatingFaceMeltingFlux = blockSumFaceMeltingFlux + floatingFaceMeltingThickness(iCell) * & + areaCell(iCell) * rhoi / ( deltat / scyr ) ! mass loss due to face-melting (kg yr^{-1}) blockSumFaceMeltingFlux = blockSumFaceMeltingFlux + faceMeltingThickness(iCell) * & @@ -531,25 +562,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 +611,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 +645,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 From 16a7b6ef48ad434d9058148caf16401afd20f70b Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 17:02:51 -0600 Subject: [PATCH 06/32] Zero out grounded and floating mass budget masks at start of call Zero out grounded and floating mass budget masks each time li_calculate_mask is called, to avoid inheritting values from previous steps. Also fix the calculation of total grounded and floating face-melting fluxes in global stats. --- .../src/analysis_members/mpas_li_global_stats.F | 4 ++-- components/mpas-albany-landice/src/shared/mpas_li_mask.F | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) 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 19cfa96ae7ba..2fb00154cc15 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 @@ -437,9 +437,9 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) ! mass loss due to calving (kg yr^{-1}) blockSumFaceMeltingFlux = blockSumFaceMeltingFlux + faceMeltingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) - blockSumGroundedFaceMeltingFlux = blockSumFaceMeltingFlux + groundedFaceMeltingThickness(iCell) * & + blockSumGroundedFaceMeltingFlux = blockSumGroundedFaceMeltingFlux + groundedFaceMeltingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) - blockSumFloatingFaceMeltingFlux = blockSumFaceMeltingFlux + floatingFaceMeltingThickness(iCell) * & + blockSumFloatingFaceMeltingFlux = blockSumFloatingFaceMeltingFlux + floatingFaceMeltingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) ! mass loss due to face-melting (kg yr^{-1}) 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 fbb1bb1fad98..1628517a02cc 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -376,6 +376,8 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) ! 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) From 08e8580d73f8eea994244ced3c31bde5656a95f5 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 21:17:54 -0600 Subject: [PATCH 07/32] Account for multiple updates of calvingThickness within a timestep Add a very short subroutine update_calving_budget that is called to recalculate groundedCalvingThickness and floatingCalvingThickness after calvingThickness is applied and before masks are updated. Also treat grounded and floating calving specially within the remove_icebergs and remove_small_islands routines to deal with the face that calvingThickness is updated multiple times per timestep. Thus, it is okay for a cell to have both nonzero groundedCalvingThickness and floatingCalvingThickness. --- .../src/mode_forward/mpas_li_calving.F | 200 ++++++++++++++++-- 1 file changed, 188 insertions(+), 12 deletions(-) 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 ec8fab7faafd..b8a04089e226 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 @@ -247,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 @@ -366,8 +384,6 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) - call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) - call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) ! In data calving mode we just calculate what should be calved but don't actually calve it. @@ -668,6 +684,7 @@ subroutine li_restore_calving_front(domain, err) block => block % next enddo + call update_calving_budget(domain) ! Update mask and geometry block => domain % blocklist do while (associated(block)) @@ -972,6 +989,12 @@ subroutine thickness_calving(domain, calvingFraction, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) +<<<<<<< HEAD +======= + call update_calving_budget(domain) + call remove_small_islands(meshPool, geometryPool) + +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -1046,6 +1069,12 @@ subroutine floating_calving(domain, calvingFraction, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) +<<<<<<< HEAD +======= + call update_calving_budget(domain) + call remove_small_islands(meshPool, geometryPool) + +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -1071,11 +1100,20 @@ 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 +<<<<<<< HEAD real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine) +======= + real(kind=RKIND), pointer :: config_sea_level +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) 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 :: 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 @@ -1105,6 +1143,10 @@ subroutine remove_small_islands(domain, err) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) 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) @@ -1170,6 +1212,7 @@ subroutine remove_small_islands(domain, err) exit endif enddo +<<<<<<< HEAD endif enddo @@ -1219,6 +1262,30 @@ subroutine remove_small_islands(domain, err) (.not. any(connectedCellsList == kCell)) ) then removeIsland = .false. exit +======= + if ((nIceNeighbors == 0) .and. (nOpenOceanNeighbors == nEdgesOnCell(iCell))) then + ! If this is a single cell of ice surrounded by open ocean, kill this location + calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) + thickness(iCell) = 0.0_RKIND + ! Update calving budgets + if (groundedMaskForMassBudget(iCell) .eq. 1) then + groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) + elseif (floatingMaskForMassBudget(iCell) .eq. 1) then + floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) + endif + elseif (nIceNeighbors == 1) then + ! check if this neighbor has any additional neighbors with ice + nIceNeighbors2 = 0 + nOpenOceanNeighbors2 = 0 + do n = 1, nEdgesOnCell(neighborWithIce) + jCell = cellsOnCell(n, neighborWithIce) + if (li_mask_is_ice(cellMask(jCell))) then + nIceNeighbors2 = nIceNeighbors2 + 1 + endif + if (.not. li_mask_is_ice(cellMask(jCell)) .and. bedTopography(jCell) < config_sea_level) then + nOpenOceanNeighbors2 = nOpenOceanNeighbors2 + 1 +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) endif enddo enddo @@ -1268,6 +1335,7 @@ subroutine remove_small_islands(domain, err) calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) thickness(iCell) = 0.0_RKIND +<<<<<<< HEAD endif enddo call mpas_log_write("***Finished cleaning up after removing small islands***") @@ -1276,6 +1344,30 @@ subroutine remove_small_islands(domain, err) islandMask) call mpas_deallocate_scratch_field(seedMaskField, single_block_in=.true.) call mpas_deallocate_scratch_field(growMaskField, single_block_in=.true.) +======= + ! Update calving budgets + if (groundedMaskForMassBudget(iCell) .eq. 1) then + groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) + elseif (floatingMaskForMassBudget(iCell) .eq. 1) then + floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) + endif + + calvingThickness(neighborWithIce) = calvingThickness(neighborWithIce) + thickness(neighborWithIce) + calvingThicknessFromThreshold(neighborWithIce) = calvingThicknessFromThreshold(neighborWithIce) + thickness(neighborWithIce) + thickness(neighborWithIce) = 0.0_RKIND + ! Update calving budgets + if (groundedMaskForMassBudget(neighborWithIce) .eq. 1) then + groundedCalvingThickness(neighborWithIce) = groundedCalvingThickness(neighborWithIce) + thickness(neighborWithIce) + elseif (floatingMaskForMassBudget(neighborWithIce) .eq. 1) then + floatingCalvingThickness(neighborWithIce) = floatingCalvingThickness(neighborWithIce) + thickness(neighborWithIce) + endif + endif + + endif ! check on nIceNeighbors + + endif ! check if iCell has ice + end do ! loop over cells +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) end subroutine remove_small_islands @@ -1353,6 +1445,12 @@ subroutine topographic_calving(domain, calvingFraction, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) +<<<<<<< HEAD +======= + call update_calving_budget(domain) + call remove_small_islands(meshPool, geometryPool) + +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -1533,6 +1631,7 @@ subroutine eigencalving(domain, err) call mpas_timer_stop("halo updates") ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) @@ -1557,7 +1656,7 @@ subroutine eigencalving(domain, err) endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -1582,6 +1681,11 @@ subroutine eigencalving(domain, err) endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step +<<<<<<< HEAD +======= + call update_calving_budget(domain) + call remove_small_islands(meshPool, geometryPool) +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -1715,7 +1819,7 @@ subroutine specified_calving_velocity(domain, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -1739,7 +1843,7 @@ subroutine specified_calving_velocity(domain, err) endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -1764,6 +1868,11 @@ subroutine specified_calving_velocity(domain, err) endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step +<<<<<<< HEAD +======= + call update_calving_budget(domain) + call remove_small_islands(meshPool, geometryPool) +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -2147,11 +2256,16 @@ subroutine von_Mises_calving(domain, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) +<<<<<<< HEAD +======= + call remove_small_islands(meshPool, geometryPool) + +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo ! associated(block) @@ -2417,7 +2531,7 @@ subroutine ismip6_retreat(domain, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -3448,7 +3562,7 @@ subroutine damage_calving(domain, err) ! === apply calving === thickness(:) = thickness(:) - calvingThickness(:) - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -3464,7 +3578,7 @@ subroutine damage_calving(domain, err) thickness(iCell) = 0.0_RKIND endif enddo - + call update_calving_budget(domain) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -3489,7 +3603,12 @@ subroutine damage_calving(domain, err) endif enddo ! TODO: global reduce & reporting on amount of calving generated in this step +<<<<<<< HEAD +======= + call update_calving_budget(domain) + call remove_small_islands(meshPool, geometryPool) +>>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -4257,10 +4376,14 @@ 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 :: 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 @@ -4358,6 +4481,10 @@ subroutine li_remove_icebergs(domain) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) 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 @@ -4369,6 +4496,11 @@ subroutine li_remove_icebergs(domain) calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) ! remove any remaining ice here calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) ! remove any remaining ice here thickness(iCell) = 0.0_RKIND + if (groundedMaskForMassBudget(iCell) .eq. 1) then + groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) + elseif (floatingMaskForMassBudget(iCell) .eq. 1) then + floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) + endif localIcebergCellCount = localIcebergCellCount + 1 !seedMaskOld(iCell) = 1 ! debug: make this a mask of where icebergs were removed endif @@ -4544,6 +4676,50 @@ subroutine li_flood_fill(seedMask, growMask, domain) end subroutine li_flood_fill + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! routine update_calving_budget +! +!> \brief Keep a running total of calving applied to grounded and floating +!> ice, respectively. +!> \author Trevor Hillebrand +!> \date May 2022 +!> \details This routine should be called after each time time calvingThickness +!> is applied, but before masks are updated which often happens multiple times +!> in a timestep. +!----------------------------------------------------------------------- + subroutine update_calving_budget(domain) + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + type (mpas_pool_type), pointer :: meshPool, geometryPool + integer, dimension(:), pointer :: groundedMaskForMassBudget, & ! binary masks for mass budget + floatingMaskForMassBudget + real (kind=RKIND), dimension(:), pointer :: calvingThickness, & + groundedCalvingThickness, & ! Grounded and floating components for mass budget + floatingCalvingThickness + + call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) + call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) + call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + + where (groundedMaskForMassBudget .eq. 1) + groundedCalvingThickness = calvingThickness + elsewhere (floatingMaskForMassBudget .eq. 1) + floatingCalvingThickness = calvingThickness + end where + + end subroutine update_calving_budget + end module li_calving From a18e7f1762a726f7e5632bff18ae0f07e2d10b8e Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 17 May 2022 17:55:27 -0600 Subject: [PATCH 08/32] Use grounded and floating mass budget masks for groundedToFloatingThickness Instead of cellMask before and after advection, use groundedMaskForMassBudget and floatingMaskForMassBudget before and after advection to calculate groundedToFloatingThickness. This prevents the switch between floating non-dynamic cells and grounded cells being registered as ice grounding, which was causing groundedLineMigration to be the wrong sign in testing on the Humboldt domain. --- .../mpas-albany-landice/src/Registry.xml | 8 +++- .../src/mode_forward/mpas_li_advection.F | 48 ++++++++++++++----- 2 files changed, 43 insertions(+), 13 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index fc122bdf13cb..bdbfeb431684 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1472,8 +1472,8 @@ 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/mode_forward/mpas_li_advection.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F index e1e1b03dc6f1..102d0c18a570 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 @@ -198,7 +198,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 @@ -364,8 +366,10 @@ subroutine li_advection_thickness_tracers(& call mpas_allocate_scratch_field(basalTracersField, .true.) basalTracers => basalTracersField % array - call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField) - call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.) + 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, 'edgeMaskTemporary', edgeMaskTemporaryField) call mpas_allocate_scratch_field(edgeMaskTemporaryField, .true.) @@ -388,6 +392,9 @@ subroutine li_advection_thickness_tracers(& vertexMaskTemporaryField % array(:) = vertexMask(:) layerThicknessEdgeFlux(:,:) = 0.0_RKIND + ! save old mass budget masks for determining cells changing from grounded to floating and vice versa + groundedMaskForMassBudgetTemporaryField % array(:) = groundedMaskForMassBudget(:) + floatingMaskForMassBudgetTemporaryField % array(:) = floatingMaskForMassBudget(:) !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers @@ -651,6 +658,16 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + ! Calculate volume converted from grounded to floating + ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state + call grounded_to_floating(groundedMaskForMassBudgetTemporaryField % array, & + floatingMaskForMassBudgetTemporaryField % array, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & + thickness, & + groundedToFloatingThickness, & + nCells) + ! Calculate flux across grounding line ! Do this after new thickness & mask have been calculated, including SMB/BMB fluxAcrossGroundingLine(:) = 0.0_RKIND @@ -756,6 +773,8 @@ subroutine li_advection_thickness_tracers(& call mpas_deallocate_scratch_field(edgeMaskTemporaryField, .true.) call mpas_deallocate_scratch_field(vertexMaskTemporaryField, .true.) call mpas_deallocate_scratch_field(thermalCellMaskField, .true.) + call mpas_deallocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.) + call mpas_deallocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.) ! === error check if (err > 0) then @@ -2129,16 +2148,23 @@ subroutine vertical_remap(thickness, meshPool, layerThickness, tracers, err) end subroutine vertical_remap - subroutine li_grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells) + subroutine grounded_to_floating(groundedMaskForMassBudgetOrig, & + floatingMaskForMassBudgetOrig, & + groundedMaskForMassBudgetNew, & + floatingMaskForMassBudgetNew, & + thicknessNew, & + groundedToFloatingThickness, & + nCells) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- integer, dimension(:), intent(in) :: & - cellMaskOrig !< Input: mask for cells before advection + groundedMaskForMassBudgetOrig, & !< Input: mask for cells before advection + floatingMaskForMassBudgetOrig integer, dimension(:), intent(in) :: & - cellMaskNew !< Input: mask for cells after advection - + groundedMaskForMassBudgetNew, & !< Input: mask for cells after advection + floatingMaskForMassBudgetNew real(kind=RKIND), dimension(:), intent(in) :: & thicknessNew !< Input: ice thickness after advection @@ -2163,12 +2189,12 @@ subroutine li_grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, grou 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 + if ( (groundedMaskForMassBudgetOrig(iCell) .eq. 1) .and. & + (floatingMaskForMassBudgetNew(iCell) .eq. 1) ) then groundedToFloatingThickness(iCell) = thicknessNew(iCell) ! find locations that had floating ice before but now have grounded ice - else if (li_mask_is_floating_ice(cellMaskOrig(iCell)) .and. & - li_mask_is_grounded_ice(cellMaskNew(iCell)) ) then + else if ( (floatingMaskForMassBudgetOrig(iCell) .eq. 1) .and. & + (groundedMaskForMassBudgetNew(iCell) .eq. 1) ) then groundedToFloatingThickness(iCell) = -1.0_RKIND * thicknessNew(iCell) endif enddo From ad9f96d6c83f0c2fe29e9eb936096eca4e0fd7aa Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 May 2022 17:21:31 -0600 Subject: [PATCH 09/32] Calculate grounding line flux before advection, smb, and bmb Calculating grounding line flux after advection led to double (or triple?) counting between GL Flux and GL Migration Flux (and maybe SMB as well). Calculate grounding line flux before advection, smb, and bmb instead. --- .../src/mode_forward/mpas_li_advection.F | 76 ++++++++----------- 1 file changed, 32 insertions(+), 44 deletions(-) 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 102d0c18a570..f7709d37a6ed 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 @@ -403,6 +403,38 @@ subroutine li_advection_thickness_tracers(& if ( (trim(config_thickness_advection) .eq. 'fo') .or. & (trim(config_thickness_advection) .eq. 'fct') ) then + ! Calculate flux across grounding line before applying + ! sfc/basal mass balance and advection + 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 + do k = 1, nVertLevels + thicknessFluxEdge = layerNormalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge) + fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * thicknessFluxEdge / dvEdge(iEdge) + enddo + ! 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 + if (config_print_thickness_advection_info) then call mpas_log_write('Using first-order upwind for thickness advection') call mpas_log_write('advectTracers = $l', logicArgs=(/advectTracers/)) @@ -667,50 +699,6 @@ subroutine li_advection_thickness_tracers(& thickness, & groundedToFloatingThickness, & nCells) - - ! 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. From f8a266250565537a544644f7c10c81b5be1bc051 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 May 2022 18:38:19 -0600 Subject: [PATCH 10/32] Update groundedToFloatingThickness when masks are updated Update groundedToFloatingThickness whenever masks are updated to account for cells changing between masks when face-melting and calving are applied, and not just when sfcMassBal, basalMassBal, and advection are applied, as was the case previously. --- .../src/mode_forward/mpas_li_advection.F | 73 +------------------ .../mpas_li_time_integration_fe_rk.F | 6 ++ .../src/shared/mpas_li_mask.F | 28 ++++++- 3 files changed, 34 insertions(+), 73 deletions(-) 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 f7709d37a6ed..2d4bbbfbd72b 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 !-------------------------------------------------------------------- ! @@ -392,10 +391,6 @@ subroutine li_advection_thickness_tracers(& vertexMaskTemporaryField % array(:) = vertexMask(:) layerThicknessEdgeFlux(:,:) = 0.0_RKIND - ! save old mass budget masks for determining cells changing from grounded to floating and vice versa - groundedMaskForMassBudgetTemporaryField % array(:) = groundedMaskForMassBudget(:) - floatingMaskForMassBudgetTemporaryField % array(:) = floatingMaskForMassBudget(:) - !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers !----------------------------------------------------------------- @@ -599,8 +594,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 !----------------------------------------------------------------- @@ -690,15 +683,6 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! Calculate volume converted from grounded to floating - ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state - call grounded_to_floating(groundedMaskForMassBudgetTemporaryField % array, & - floatingMaskForMassBudgetTemporaryField % array, & - groundedMaskForMassBudget, & - floatingMaskForMassBudget, & - thickness, & - groundedToFloatingThickness, & - nCells) ! Remap tracers to the standard vertical sigma coordinate ! Note: If tracers are not being advected, then this subroutine simply restores the ! layer thickness to sigma coordinate values. @@ -761,8 +745,6 @@ subroutine li_advection_thickness_tracers(& call mpas_deallocate_scratch_field(edgeMaskTemporaryField, .true.) call mpas_deallocate_scratch_field(vertexMaskTemporaryField, .true.) call mpas_deallocate_scratch_field(thermalCellMaskField, .true.) - call mpas_deallocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.) - call mpas_deallocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.) ! === error check if (err > 0) then @@ -2136,59 +2118,6 @@ subroutine vertical_remap(thickness, meshPool, layerThickness, tracers, err) end subroutine vertical_remap - subroutine grounded_to_floating(groundedMaskForMassBudgetOrig, & - floatingMaskForMassBudgetOrig, & - groundedMaskForMassBudgetNew, & - floatingMaskForMassBudgetNew, & - thicknessNew, & - groundedToFloatingThickness, & - nCells) - !----------------------------------------------------------------- - ! input variables - !----------------------------------------------------------------- - integer, dimension(:), intent(in) :: & - groundedMaskForMassBudgetOrig, & !< Input: mask for cells before advection - floatingMaskForMassBudgetOrig - - integer, dimension(:), intent(in) :: & - groundedMaskForMassBudgetNew, & !< Input: mask for cells after advection - floatingMaskForMassBudgetNew - real(kind=RKIND), dimension(:), intent(in) :: & - thicknessNew !< Input: ice thickness after advection - - integer, pointer, intent(in) :: & - nCells - !----------------------------------------------------------------- - ! input/output variables - !----------------------------------------------------------------- - - !----------------------------------------------------------------- - ! output variables - !----------------------------------------------------------------- - real(kind=RKIND), dimension(:), intent(out) :: & - groundedToFloatingThickness !< Input: ice thickness - - !----------------------------------------------------------------- - ! local variables - !----------------------------------------------------------------- - integer :: iCell - - groundedToFloatingThickness(:) = 0.0_RKIND - - do iCell = 1, nCells - ! find locations that had grounded ice before but now have floating ice - if ( (groundedMaskForMassBudgetOrig(iCell) .eq. 1) .and. & - (floatingMaskForMassBudgetNew(iCell) .eq. 1) ) then - groundedToFloatingThickness(iCell) = thicknessNew(iCell) - ! find locations that had floating ice before but now have grounded ice - else if ( (floatingMaskForMassBudgetOrig(iCell) .eq. 1) .and. & - (groundedMaskForMassBudgetNew(iCell) .eq. 1) ) then - groundedToFloatingThickness(iCell) = -1.0_RKIND * thicknessNew(iCell) - endif - enddo - - end subroutine li_grounded_to_floating - !*********************************************************************** end module li_advection 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 3c429eb8f762..01d7540530a3 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 @@ -697,6 +697,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 +791,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( & 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 1628517a02cc..05eaefc9670e 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -297,7 +297,7 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) !----------------------------------------------------------------- integer, pointer :: nCells, nVertices, nEdges, vertexDegree integer, pointer :: nVertInterfaces - real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography + real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography, groundedToFloatingThickness integer, dimension(:), pointer :: nEdgesOnCell, cellMask, vertexMask, edgeMask, & groundedMaskForMassBudget, floatingMaskForMassBudget integer, dimension(:,:), pointer :: cellsOnCell, cellsOnVertex, cellsOnEdge, dirichletVelocityMask @@ -322,6 +322,9 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) integer :: iCellNeighbor logical :: openOceanNeighbor logical :: updateMassBudgetMasks=.true. + type (field1DInteger), pointer :: & ! used to calculate groundedToFloatingThickness + groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget + floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget call mpas_timer_start('calculate mask') @@ -344,9 +347,15 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + call mpas_pool_get_field(geometryPool, 'groundedMaskForMassBudgetTemporary', groundedMaskForMassBudgetTemporaryField) + call mpas_pool_get_field(geometryPool, 'floatingMaskForMassBudgetTemporary', floatingMaskForMassBudgetTemporaryField) + call mpas_allocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.) + call mpas_allocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.) + call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) @@ -372,6 +381,11 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) cellMask = ior(cellMask, li_mask_ValueIce) end where + + ! save old mass budget masks for determining cells changing from grounded + ! to floating and vice versa + groundedMaskForMassBudgetTemporaryField % array(:) = groundedMaskForMassBudget(:) + floatingMaskForMassBudgetTemporaryField % array(:) = floatingMaskForMassBudget(:) ! Is it floating? (ice thickness equal to floatation is considered floating) ! For now floating ice and grounded ice are mutually exclusive. ! This may change if a grounding line parameterization is added. @@ -540,6 +554,18 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) endif enddo + ! Update groundedToFloatingThickness based on new masks + do iCell = 1, nCells + ! find locations that had grounded ice before but now have floating ice + if ( (groundedMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (floatingMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) + thickness(iCell) + ! find locations that had floating ice before but now have grounded ice + else if ( (floatingMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (groundedMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) - thickness(iCell) + endif + enddo ! Now use flood fill to identify all floating ice connected to the open ! ocean. Subglacial lakes are any floating ice that are not in the ! resulting flood fill mask. Commented code below does not yet work From 2cc6f924410cdbac5377859d4834fb249d310c29 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 3 Aug 2022 18:09:22 -0600 Subject: [PATCH 11/32] Correct use of 'cycle' when when calculating masks The cycle statement used to exit a loop over neighbors when calculating the grounding line and adding non-dynamic floating fringe to the grounded mask was inadvertently outside the relevant if-statement. --- components/mpas-albany-landice/src/shared/mpas_li_mask.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 05eaefc9670e..3fd0b7215213 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -535,8 +535,8 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) if (li_mask_is_grounded_ice(cellMask(iCellNeighbor))) then groundedMaskForMassBudget(i) = 1 floatingMaskForMassBudget(i) = 0 + cycle ! no need to look at additional neighbors endif - cycle ! no need to look at additional neighbors enddo endif @@ -548,8 +548,8 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) ! 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 - ! cycle ! no need to look at additional neighbors ! enddo endif enddo From de98bb07c9895817a68d3037ae9382179180cedf Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 5 Oct 2022 11:40:34 -0600 Subject: [PATCH 12/32] Do not apply float-kill to non-dynamic fringe Use floatingMaskForMassBudget to define which cells are calving when using config_calving = 'floating'. This is because floating non-dynamic cells adjacent to grounded cells should not really be considered part of the ice shelf. --- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) 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 b8a04089e226..c782347f4954 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 @@ -1043,6 +1043,7 @@ subroutine floating_calving(domain, calvingFraction, err) real (kind=RKIND), dimension(:), pointer :: calvingThicknessFromThreshold ! thickness of ice that calves (computed in this subroutine) real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask + integer, dimension(:), pointer :: floatingMaskForMassBudget err = 0 @@ -1057,11 +1058,13 @@ 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) + call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) calvingThickness = 0.0_RKIND - ! Note: The floating_ice mask includes all floating ice, both inactive and active - where (li_mask_is_floating_ice(cellMask)) + ! Note: The floating_ice mask does not include floating + ! non-dynamic cells adjacent to grounded cells. + where ( floatingMaskForMassBudget .eq. 1 ) calvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness From 9cfc6996392f3a6a77b08ac3cdf09e4d1637230e Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 5 Oct 2022 12:56:48 -0600 Subject: [PATCH 13/32] Use mass budget masks for basal mass balance Use floatingMaskForMassBudget and groundedMaskForMassBudget to partition basal mass balance between grounded and floating. --- .../src/mode_forward/mpas_li_advection.F | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) 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 2d4bbbfbd72b..560ba3a16332 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 @@ -603,18 +603,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. From c5030542700ed7450e3a2190b4cbf70f8630fca7 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 1 Nov 2022 19:30:05 -0600 Subject: [PATCH 14/32] Update calving budget incrementally Use a local variable dCalvingThickness to update the calving budget each time calving is applied. The variable calvingThickness is now not used directly to modify thickness, but calculated by summing dCalvingThickness throughout a time step. This might lead to issues with halo updates because dCalvingThickness is always an allocatable array. If so, it will need to either be made a Registry variable or an MPAS allocatable scratch array. --- .../src/mode_forward/mpas_li_calving.F | 349 +++++++----------- 1 file changed, 133 insertions(+), 216 deletions(-) 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 c782347f4954..d59a8b3e4147 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 @@ -516,6 +516,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(:), allocatable:: & + dCalvingThickness ! incremental calving thickness real (kind=RKIND) :: & restoreThicknessMin ! small thickness to which ice is restored should it fall below this thickness @@ -661,6 +663,9 @@ subroutine li_restore_calving_front(domain, err) endif ! if preventing retreat ! Now check for marine regions that have advanced + allocate(dCalvingThickness(nCellsSolve+1)) + + ! 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. @@ -673,7 +678,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) + dCalvingThickness(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) @@ -681,10 +686,12 @@ subroutine li_restore_calving_front(domain, err) endif ! bedTopography < config_sea_level enddo ! iCell + call update_calving_budget(geometryPool, dCalvingThickness) + + deallocate(dCalvingThickness) block => block % next enddo - call update_calving_budget(domain) ! Update mask and geometry block => domain % blocklist do while (associated(block)) @@ -825,6 +832,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(:), allocatable :: & + dCalvingThickness integer, pointer :: nCells integer :: iCell, iCellOnCell, iCellNeighbor @@ -892,7 +901,8 @@ subroutine thickness_calving(domain, calvingFraction, err) ! Initialize calvingThickness = 0.0_RKIND - + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 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 @@ -982,19 +992,15 @@ 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 + dCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - dCalvingThickness(:) -<<<<<<< HEAD -======= - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_small_islands(meshPool, geometryPool) - ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -1002,6 +1008,7 @@ subroutine thickness_calving(domain, calvingFraction, err) call mpas_deallocate_scratch_field(activeForCalvingMaskField, .true.) call mpas_deallocate_scratch_field(inactiveMarginMaskField, .true.) call mpas_deallocate_scratch_field(oceanMaskField, .true.) + deallocate(dCalvingThickness) end subroutine thickness_calving @@ -1041,9 +1048,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(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask integer, dimension(:), pointer :: floatingMaskForMassBudget + integer, pointer :: nCells err = 0 @@ -1059,25 +1068,24 @@ subroutine floating_calving(domain, calvingFraction, err) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + call mpas_pool_get_array(meshPool, 'nCells', nCells) - calvingThickness = 0.0_RKIND - + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND ! Note: The floating_ice mask does not include floating ! non-dynamic cells adjacent to grounded cells. where ( floatingMaskForMassBudget .eq. 1 ) - calvingThickness = thickness * calvingFraction + dCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - dCalvingThickness(:) -<<<<<<< HEAD -======= - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_small_islands(meshPool, geometryPool) ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) + deallocate(dCalvingThickness) block => block % next enddo @@ -1103,15 +1111,12 @@ 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 -<<<<<<< HEAD - real (kind=RKIND), dimension(:), pointer :: calvingThickness ! thickness of ice that calves (computed in this subroutine) -======= real(kind=RKIND), pointer :: config_sea_level ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) 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(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: bedTopography integer, dimension(:), pointer :: cellMask, & @@ -1155,7 +1160,8 @@ subroutine remove_small_islands(domain, err) allocate(connectedCellsList((maxEdges+1)**2), & islandMask(nCells+1)) - + allocate(dCalvingThickness(nCellsSolve+1)) + dCalvingThickness(:) = 0.0_RKIND islandMask(:) = 0 nIslandCellsLocal = 0 nIslandCellsGlobal = 0 @@ -1215,7 +1221,6 @@ subroutine remove_small_islands(domain, err) exit endif enddo -<<<<<<< HEAD endif enddo @@ -1265,30 +1270,6 @@ subroutine remove_small_islands(domain, err) (.not. any(connectedCellsList == kCell)) ) then removeIsland = .false. exit -======= - if ((nIceNeighbors == 0) .and. (nOpenOceanNeighbors == nEdgesOnCell(iCell))) then - ! If this is a single cell of ice surrounded by open ocean, kill this location - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) - calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) - thickness(iCell) = 0.0_RKIND - ! Update calving budgets - if (groundedMaskForMassBudget(iCell) .eq. 1) then - groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) - elseif (floatingMaskForMassBudget(iCell) .eq. 1) then - floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) - endif - elseif (nIceNeighbors == 1) then - ! check if this neighbor has any additional neighbors with ice - nIceNeighbors2 = 0 - nOpenOceanNeighbors2 = 0 - do n = 1, nEdgesOnCell(neighborWithIce) - jCell = cellsOnCell(n, neighborWithIce) - if (li_mask_is_ice(cellMask(jCell))) then - nIceNeighbors2 = nIceNeighbors2 + 1 - endif - if (.not. li_mask_is_ice(cellMask(jCell)) .and. bedTopography(jCell) < config_sea_level) then - nOpenOceanNeighbors2 = nOpenOceanNeighbors2 + 1 ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) endif enddo enddo @@ -1335,10 +1316,9 @@ 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) + dCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND -<<<<<<< HEAD endif enddo call mpas_log_write("***Finished cleaning up after removing small islands***") @@ -1347,31 +1327,10 @@ subroutine remove_small_islands(domain, err) islandMask) call mpas_deallocate_scratch_field(seedMaskField, single_block_in=.true.) call mpas_deallocate_scratch_field(growMaskField, single_block_in=.true.) -======= - ! Update calving budgets - if (groundedMaskForMassBudget(iCell) .eq. 1) then - groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) - elseif (floatingMaskForMassBudget(iCell) .eq. 1) then - floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) - endif - - calvingThickness(neighborWithIce) = calvingThickness(neighborWithIce) + thickness(neighborWithIce) - calvingThicknessFromThreshold(neighborWithIce) = calvingThicknessFromThreshold(neighborWithIce) + thickness(neighborWithIce) - thickness(neighborWithIce) = 0.0_RKIND - ! Update calving budgets - if (groundedMaskForMassBudget(neighborWithIce) .eq. 1) then - groundedCalvingThickness(neighborWithIce) = groundedCalvingThickness(neighborWithIce) + thickness(neighborWithIce) - elseif (floatingMaskForMassBudget(neighborWithIce) .eq. 1) then - floatingCalvingThickness(neighborWithIce) = floatingCalvingThickness(neighborWithIce) + thickness(neighborWithIce) - endif - endif - - endif ! check on nIceNeighbors - - endif ! check if iCell has ice - end do ! loop over cells ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) + call update_calving_budget(geometryPool, dCalvingThickness) + + deallocate(dCalvingThickness) end subroutine remove_small_islands @@ -1409,11 +1368,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(:), allocatable :: dCalvingThickness 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 @@ -1437,23 +1398,24 @@ subroutine topographic_calving(domain, calvingFraction, err) 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 + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND where ( (li_mask_is_floating_ice(cellMask)) .and. (bedTopography < config_calving_topography + config_sea_level) ) - calvingThickness = thickness * calvingFraction + dCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - dCalvingThickness(:) -<<<<<<< HEAD -======= - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_small_islands(meshPool, geometryPool) ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) + deallocate(dCalvingThickness) + block => block % next enddo @@ -1511,6 +1473,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(:), allocatable :: dCalvingThickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell @@ -1592,6 +1555,8 @@ subroutine eigencalving(domain, err) realArgs=(/minval(eigencalvingParameter), maxval(eigencalvingParameter)/)) endif + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 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 @@ -1599,7 +1564,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, & + dCalvingThickness, calvingVelocity, applyToGrounded, & applyToFloating, applyToGroundingLine, domain, calvingCFLdt, dtCalvingCFLratio, & totalRatebasedCalvedVolume, totalRatebasedUncalvedVolume, err_tmp) err = ior(err, err_tmp) @@ -1611,7 +1576,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, dCalvingThickness, 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 & @@ -1633,8 +1598,8 @@ subroutine eigencalving(domain, err) call mpas_dmpar_field_halo_exch(domain, 'calvingThickness') call mpas_timer_stop("halo updates") ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - call update_calving_budget(domain) + thickness(:) = thickness(:) - dCalvingThickness(:) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) @@ -1654,42 +1619,21 @@ 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) + dCalvingThickness(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(domain) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask 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 -<<<<<<< HEAD -======= call update_calving_budget(domain) + call remove_icebergs(domain) call remove_small_islands(meshPool, geometryPool) ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) + deallocate(dCalvingThickness) block => block % next enddo @@ -1741,6 +1685,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(:), allocatable :: dCalvingThickness ! Incremental calving thickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell @@ -1794,6 +1739,8 @@ subroutine specified_calving_velocity(domain, err) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND ! get parameter value if (trim(config_calving_specified_source) == 'const') then @@ -1812,7 +1759,7 @@ 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, dCalvingThickness, calvingVelocity, & applyToGrounded=.true., applyToFloating=.true., applyToGroundingLine=.false., & domain=domain, maxDt=calvingCFLdt, CFLratio=dtCalvingCFLratio, & totalAblatedVolume=totalRatebasedCalvedVolume, & @@ -1821,8 +1768,8 @@ subroutine specified_calving_velocity(domain, err) err = ior(err, err_tmp) ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - call update_calving_budget(domain) + thickness(:) = thickness(:) - dCalvingThickness(:) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -1841,43 +1788,22 @@ 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) + dCalvingThickness(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(domain) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask 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 -<<<<<<< HEAD -======= call update_calving_budget(domain) - call remove_small_islands(meshPool, geometryPool) ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next + + deallocate(dCalvingThickness) enddo end subroutine specified_calving_velocity @@ -1938,6 +1864,7 @@ subroutine von_Mises_calving(domain, err) floatingVonMisesThresholdStress, & groundedVonMisesThresholdStress, & surfaceSpeed + real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:,:), pointer :: flowParamA, & temperature, layerThickness real (kind=RKIND), pointer :: config_default_flowParamA @@ -2050,7 +1977,9 @@ subroutine von_Mises_calving(domain, err) endif vonMisesStress(:) = 0.0_RKIND - + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 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 @@ -2186,7 +2115,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, & + dCalvingThickness, calvingVelocity, applyToGrounded, & applyToFloating, applyToGroundingLine, domain, & calvingCFLdt, dtCalvingCFLratio, & totalRatebasedCalvedVolume, totalRatebasedUncalvedVolume, err_tmp) @@ -2210,7 +2139,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) + dCalvingThickness(iCell) = thickness(iCell) elseif ( li_mask_is_floating_ice(cellMask(icell)) .and. & (.not. li_mask_is_dynamic_ice(cellMask(iCell))) ) then nGroundedNeighbors = 0 @@ -2221,12 +2150,19 @@ subroutine von_Mises_calving(domain, err) endif enddo if ( nGroundedNeighbors == 0 ) then - calvingThickness(iCell) = thickness(iCell) + dCalvingThickness(iCell) = thickness(iCell) endif endif enddo endif + ! === apply calving velocity before damage threshold calving === + thickness(:) = thickness(:) - dCalvingThickness(:) + call update_calving_budget(geometryPool, dCalvingThickness) + ! 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 @@ -2234,7 +2170,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, dCalvingThickness, 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 & @@ -2247,28 +2183,14 @@ 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(:) - call update_calving_budget(domain) + ! === apply calving velocity before damage threshold calving === + thickness(:) = thickness(:) - dCalvingThickness(:) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) -<<<<<<< HEAD -======= - call remove_small_islands(meshPool, geometryPool) - ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) + deallocate(dCalvingThickness) block => block % next enddo ! associated(block) @@ -2335,6 +2257,7 @@ subroutine ismip6_retreat(domain, err) calvingVelocity, thickness, & xvelmean, yvelmean, calvingThickness, & bedTopography + real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, pointer :: nCells, timestepNumber @@ -2395,6 +2318,9 @@ subroutine ismip6_retreat(domain, err) ! submergedArea used in runoff unit conversion allocate(submergedArea(nCells+1)) + submergedArea(:) = 0.0_RKIND + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND streamFound = .false. ! Changed to true if ismip6-gis stream is found, otherwise throws error stream_cursor => domain % streamManager % streams % head @@ -2518,7 +2444,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., & + dCalvingThickness, calvingVelocity, applyToGrounded=.true., & applyToFloating=.true., applyToGroundingLine=.false., & domain=domain, maxDt=calvingCFLdt, CFLratio=dtCalvingCFLratio, & totalAblatedVolume=totalRatebasedCalvedVolume, & @@ -2533,14 +2459,14 @@ subroutine ismip6_retreat(domain, err) call mpas_timer_stop("halo updates") ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - call update_calving_budget(domain) + thickness(:) = thickness(:) - dCalvingThickness(:) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) deallocate(submergedArea) - + deallocate(dCalvingThickness) end subroutine ismip6_retreat !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -3473,6 +3399,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(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:), pointer :: calvingVelocity real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio @@ -3536,6 +3463,8 @@ subroutine damage_calving(domain, err) call mpas_pool_get_array(geometryPool, 'totalRatebasedCalvedVolume', totalRatebasedCalvedVolume) call mpas_pool_get_array(geometryPool, 'totalRatebasedUncalvedVolume', totalRatebasedUncalvedVolume) + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND call calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask) @@ -3557,15 +3486,15 @@ 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, dCalvingThickness, 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 ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) - call update_calving_budget(domain) + thickness(:) = thickness(:) - dCalvingThickness(:) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -3577,41 +3506,15 @@ 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) + dCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) ! update mask 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 -<<<<<<< HEAD - -======= - call update_calving_budget(domain) - call remove_small_islands(meshPool, geometryPool) ->>>>>>> d85d0dc7a7 (Account for multiple updates of calvingThickness within a timestep) block => block % next enddo @@ -4085,7 +3988,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, dCalvingThickness, err) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- @@ -4097,12 +4000,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) :: dCalvingThickness ! incremental calving thickness !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -4143,6 +4045,8 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d growMask => growMaskField % array growMask(:) = 0 + dCalvingThickness(:) = 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 @@ -4162,13 +4066,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) + dCalvingThickness(jCell) = thickness(jCell) calvingThicknessFromThreshold(jCell) = calvingThicknessFromThreshold(jCell) + thickness(jCell) endif enddo - calvingThickness(iCell) = thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) + dCalvingThickness(iCell) = thickness(iCell) endif ! if cell is calving margin enddo ! cell loop @@ -4223,6 +4127,7 @@ subroutine mask_calving(domain, err) integer, pointer :: nCells, nCellsSolve integer :: iCell integer :: localMaskCellCount, globalMaskCellCount + real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness integer :: err_tmp err = 0 @@ -4244,6 +4149,10 @@ subroutine mask_calving(domain, err) 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) + + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) @@ -4255,7 +4164,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) + dCalvingThickness(iCell) = thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) thickness(iCell) = 0.0_RKIND if (iCell <= nCellsSolve) localMaskCellCount = localMaskCellCount + 1 @@ -4265,10 +4174,13 @@ 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, dCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + + deallocate(dCalvingThickness) block => block % next enddo @@ -4383,6 +4295,7 @@ subroutine li_remove_icebergs(domain) real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine) groundedCalvingThickness, & floatingCalvingThickness + real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask, & groundedMaskForMassBudget, & @@ -4423,6 +4336,9 @@ subroutine li_remove_icebergs(domain) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND + seedMask(:) = 0 growMask(:) = 0 @@ -4496,19 +4412,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 + dCalvingThickness(iCell) = thickness(iCell) ! remove any remaining ice here thickness(iCell) = 0.0_RKIND - if (groundedMaskForMassBudget(iCell) .eq. 1) then - groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) - elseif (floatingMaskForMassBudget(iCell) .eq. 1) then - floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) - endif localIcebergCellCount = localIcebergCellCount + 1 !seedMaskOld(iCell) = 1 ! debug: make this a mask of where icebergs were removed endif enddo + call update_calving_budget(geometryPool, dCalvingThickness) block => block % next end do @@ -4526,7 +4438,10 @@ 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 + + deallocate(dCalvingThickness) + end subroutine remove_icebergs +>>>>>>> a5bee82292 (Update calving budget incrementally) !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -4692,23 +4607,22 @@ end subroutine li_flood_fill !> is applied, but before masks are updated which often happens multiple times !> in a timestep. !----------------------------------------------------------------------- - subroutine update_calving_budget(domain) + subroutine update_calving_budget(geometryPool, dCalvingThickness) !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- - type (domain_type), intent(inout) :: domain + type (mpas_pool_type), pointer, intent(inout) :: geometryPool + real (kind=RKIND), dimension(:), intent(inout) :: dCalvingThickness !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- - type (mpas_pool_type), pointer :: meshPool, geometryPool integer, dimension(:), pointer :: groundedMaskForMassBudget, & ! binary masks for mass budget floatingMaskForMassBudget real (kind=RKIND), dimension(:), pointer :: calvingThickness, & groundedCalvingThickness, & ! Grounded and floating components for mass budget floatingCalvingThickness - call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) @@ -4716,10 +4630,13 @@ subroutine update_calving_budget(domain) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) where (groundedMaskForMassBudget .eq. 1) - groundedCalvingThickness = calvingThickness + groundedCalvingThickness = dCalvingThickness elsewhere (floatingMaskForMassBudget .eq. 1) - floatingCalvingThickness = calvingThickness - end where + floatingCalvingThickness = dCalvingThickness + end where + + calvingThickness(:) = calvingThickness(:) + dCalvingThickness(:) + dCalvingThickness(:) = 0.0_RKIND end subroutine update_calving_budget From b9e6f3f99f55d6427ff9f36d4b03dcf42bbfdc36 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 1 Nov 2022 20:21:56 -0600 Subject: [PATCH 15/32] Cleanup after rebase onto MALI-Dev/develop Clean up a few items after rebasing, including changing a few instances of 'flood_fill' to 'li_flood_fill', an instance of 'nCellsSolve' to 'nCells', and removing cellMaskTemporaryField from li_advection_thickness_tracers. --- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 d59a8b3e4147..1b4b6a2b75b2 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 @@ -663,7 +663,7 @@ subroutine li_restore_calving_front(domain, err) endif ! if preventing retreat ! Now check for marine regions that have advanced - allocate(dCalvingThickness(nCellsSolve+1)) + allocate(dCalvingThickness(nCells+1)) ! loop over locally owned cells do iCell = 1, nCells From 898444ce6b06b55f881724c486509b6369860ca6 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 3 Nov 2022 10:21:51 -0600 Subject: [PATCH 16/32] Update grounded and floating calving budgets incrementally Update grounded and floating calving budgets incrementally, which is required to close budget. --- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 1b4b6a2b75b2..51db001b9093 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 @@ -4630,9 +4630,9 @@ subroutine update_calving_budget(geometryPool, dCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) where (groundedMaskForMassBudget .eq. 1) - groundedCalvingThickness = dCalvingThickness + groundedCalvingThickness = groundedCalvingThickness + dCalvingThickness elsewhere (floatingMaskForMassBudget .eq. 1) - floatingCalvingThickness = dCalvingThickness + floatingCalvingThickness = floatingCalvingThickness + dCalvingThickness end where calvingThickness(:) = calvingThickness(:) + dCalvingThickness(:) From b8236133b9a0192f9138f34d6cad0b72269cb3f1 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 23 Feb 2023 16:10:59 -0700 Subject: [PATCH 17/32] Clean-up after rebase Fix a few small issues after rebasing onto MALI-Dev/develop --- components/mpas-albany-landice/src/Registry.xml | 3 ++- .../src/analysis_members/mpas_li_global_stats.F | 1 - .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index bdbfeb431684..e264e5ca7ebf 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1280,6 +1280,7 @@ is the value of that variable from the *previous* time level! /> @@ -1368,7 +1369,6 @@ is the value of that variable from the *previous* time level! - /> @@ -1444,6 +1444,7 @@ is the value of that variable from the *previous* time level! /> 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 2fb00154cc15..abb29aa0b4e0 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 @@ -223,7 +223,6 @@ 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, & 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 51db001b9093..748de7981767 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 @@ -4155,7 +4155,7 @@ subroutine mask_calving(domain, err) dCalvingThickness(:) = 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 From 8a0e9f4a00ea4a7c067d29650613040552508bdd Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 17 Oct 2023 14:16:17 -0600 Subject: [PATCH 18/32] Clean up after rebase Clean up after rebase, again --- components/mpas-albany-landice/src/mode_forward/mpas_li_core.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 b74ea4efb18a..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) From 7ac08614516e2fdd9303a270c3eb5c69b0cb35dc Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 17 Oct 2023 16:26:14 -0600 Subject: [PATCH 19/32] Fix dCalvingThickness dimension in remove_small_islands Fix dCalvingThickness dimension in remove_small_islands. Must be nCells rather than nCellsSolve in size. --- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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 748de7981767..10b4726c49ce 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 @@ -664,6 +664,7 @@ subroutine li_restore_calving_front(domain, err) ! Now check for marine regions that have advanced allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND ! loop over locally owned cells do iCell = 1, nCells @@ -1147,7 +1148,6 @@ subroutine remove_small_islands(domain, err) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges) - 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, 'calvingThicknessFromThreshold', calvingThicknessFromThreshold) @@ -1160,7 +1160,6 @@ subroutine remove_small_islands(domain, err) allocate(connectedCellsList((maxEdges+1)**2), & islandMask(nCells+1)) - allocate(dCalvingThickness(nCellsSolve+1)) dCalvingThickness(:) = 0.0_RKIND islandMask(:) = 0 nIslandCellsLocal = 0 From f41567730f2ab77fae05c3ddd30b7f99d1a3688c Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 Oct 2023 16:00:19 -0600 Subject: [PATCH 20/32] Make treatment of basal mass bal components consistent Grounded and floating components of sfcMassBal were using the mass balance mask fields, while basalMassBal components used cellMask. Make both of them use mass balance masks. --- .../src/mode_forward/mpas_li_advection.F | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) 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 560ba3a16332..22fe1a5b5eec 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 @@ -1070,27 +1070,20 @@ subroutine apply_mass_balance(& 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 From 364a78ac3e3f4efb760c30dc0381e8249112ed4b Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 Oct 2023 17:05:25 -0600 Subject: [PATCH 21/32] Fix accidental doubling of faceMeltingFlux Fix accidental doubling of faceMeltingFlux when calculating global stats. --- .../src/analysis_members/mpas_li_global_stats.F | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) 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 abb29aa0b4e0..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 @@ -433,7 +433,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) blockSumFloatingCalvingFlux = blockSumFloatingCalvingFlux + floatingCalvingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) - ! mass loss due to calving (kg yr^{-1}) + ! mass loss due to face-melting (kg yr^{-1}) blockSumFaceMeltingFlux = blockSumFaceMeltingFlux + faceMeltingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) blockSumGroundedFaceMeltingFlux = blockSumGroundedFaceMeltingFlux + groundedFaceMeltingThickness(iCell) * & @@ -441,10 +441,6 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) blockSumFloatingFaceMeltingFlux = blockSumFloatingFaceMeltingFlux + floatingFaceMeltingThickness(iCell) * & areaCell(iCell) * rhoi / ( deltat / scyr ) - ! mass loss due to face-melting (kg yr^{-1}) - blockSumFaceMeltingFlux = blockSumFaceMeltingFlux + faceMeltingThickness(iCell) * & - areaCell(iCell) * rhoi / ( deltat / scyr ) - ! max surface speed if (surfaceSpeed(iCell) > blockMaxSurfaceSpeed) then blockMaxSurfaceSpeed = surfaceSpeed(iCell) From 2919a3bbce77619bb65076a7e0aa54647c362021 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 19 Oct 2023 20:34:19 -0600 Subject: [PATCH 22/32] Update layerThickness halos after advection Update layerThickness halos after advection. This reduces the error in a test with 500 m/yr face-melt, so it seems to be necessary, but it does not solve the issue entirely. --- .../mpas-albany-landice/src/mode_forward/mpas_li_advection.F | 4 ++++ 1 file changed, 4 insertions(+) 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 22fe1a5b5eec..984996bd8f6a 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 @@ -584,6 +584,10 @@ 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 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 From 72e4b89208ddeccec456988e0492293973561bf6 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 20 Oct 2023 10:49:40 -0600 Subject: [PATCH 23/32] Make grounding line flux calculation consistent with FCT Make grounding line flux calculation consistent with FCT advection. Also update edgeMask halos before calculating grounding line flux. --- .../src/mode_forward/mpas_li_advection.F | 37 +++++++++++-------- 1 file changed, 22 insertions(+), 15 deletions(-) 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 984996bd8f6a..5ee3aa5dde33 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 @@ -389,38 +389,45 @@ subroutine li_advection_thickness_tracers(& cellMaskTemporaryField % array(:) = cellMask(:) edgeMaskTemporaryField % array(:) = edgeMask(:) vertexMaskTemporaryField % array(:) = vertexMask(:) + ! 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_timer_stop("halo updates") - layerThicknessEdgeFlux(:,:) = 0.0_RKIND !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers !----------------------------------------------------------------- if ( (trim(config_thickness_advection) .eq. 'fo') .or. & (trim(config_thickness_advection) .eq. 'fct') ) then - - ! Calculate flux across grounding line before applying - ! sfc/basal mass balance and advection + ! 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 + 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 + 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 + else GLfluxSign = -1.0_RKIND theGroundedCell = iCell2 endif - do k = 1, nVertLevels - thicknessFluxEdge = layerNormalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge) - fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * thicknessFluxEdge / dvEdge(iEdge) - enddo + 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 + 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 @@ -494,7 +501,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. & From 67351eb5a53730af81c0f11ed8c3e024c4b4017f Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 20 Oct 2023 10:51:00 -0600 Subject: [PATCH 24/32] Move face-melting after last RK velocity solve Move face-melting after last velocity solve in RK loop. This is necessary to close mass budgets when using face-melting. --- .../mode_forward/mpas_li_time_integration_fe_rk.F | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) 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 01d7540530a3..e111c7422a26 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 @@ -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) @@ -594,6 +587,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 From de6825dbb4da60717e0ee540939a587497722b98 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 16 Nov 2023 11:46:14 -0800 Subject: [PATCH 25/32] Clean up after rebase following Runge-Kutta merge Clean up after rebase following Runge-Kutta merge --- .../src/mode_forward/mpas_li_advection.F | 3 --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 8 +------- 2 files changed, 1 insertion(+), 10 deletions(-) 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 5ee3aa5dde33..55abcdff5268 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 @@ -81,7 +81,6 @@ subroutine li_advection_thickness_tracers(& geometryPool, & thermalPool, & scratchPool, & - domain, & err, & advectTracersIn) @@ -111,8 +110,6 @@ subroutine li_advection_thickness_tracers(& !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object - type (domain_type), intent(inout) :: domain - type (mpas_pool_type), intent(inout) :: & velocityPool !< Input/output: velocity information ! (needs to be inout for li_calculate_mask call 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 e111c7422a26..6765cfbb4c59 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 @@ -514,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) @@ -526,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(:) @@ -1166,7 +1162,6 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & - domain, & err_tmp, & advectTracersIn = .true.) @@ -1186,7 +1181,6 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & - domain, & err_tmp, & advectTracersIn = .false.) From 372287ee10cf903f2c721f2d7ae552516656f2a5 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 16 Nov 2023 14:23:32 -0700 Subject: [PATCH 26/32] Add cellMaskTemporary back to Registry Add cellMaskTemporary back to Registry for use in calving. --- components/mpas-albany-landice/src/Registry.xml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index e264e5ca7ebf..227e4d76b3a2 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1473,7 +1473,11 @@ 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." /> - + From 8aff98f939389156bc6aa00bfc20b7f761805b75 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 16 Nov 2023 15:10:03 -0700 Subject: [PATCH 27/32] Loop over nEdgeSolve instead of nEdges when calculating GL flux Loop over nEdgeSolve instead of nEdges when calculating GL flux. Looping over nEdges will give an error resulting from grounding line thickness <= 0 when the grounding line is at the domain boundary. --- .../src/mode_forward/mpas_li_advection.F | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) 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 55abcdff5268..9d9a757c10e7 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 @@ -244,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 @@ -328,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) @@ -389,6 +390,7 @@ subroutine li_advection_thickness_tracers(& ! 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") !----------------------------------------------------------------- @@ -401,7 +403,7 @@ subroutine li_advection_thickness_tracers(& ! Do this after new thickness & mask have been calculated, including SMB/BMB fluxAcrossGroundingLine(:) = 0.0_RKIND fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND - do iEdge = 1, nEdges + do iEdge = 1, nEdgesSolve if (li_mask_is_grounding_line(edgeMask(iEdge))) then iCell1 = cellsOnEdge(1,iEdge) iCell2 = cellsOnEdge(2,iEdge) @@ -434,6 +436,11 @@ subroutine li_advection_thickness_tracers(& 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") + if (config_print_thickness_advection_info) then call mpas_log_write('Using first-order upwind for thickness advection') call mpas_log_write('advectTracers = $l', logicArgs=(/advectTracers/)) From e3d145148df303f5a3a852c8778724751a23f041 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 16 Nov 2023 19:16:58 -0700 Subject: [PATCH 28/32] Calculate grounding line flux after advection but before updating thickness Calculate grounding line flux after advection but before updating thickness, which is required to close mass budgets. --- .../src/mode_forward/mpas_li_advection.F | 83 ++++++++++--------- 1 file changed, 42 insertions(+), 41 deletions(-) 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 9d9a757c10e7..da147eb5b60e 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 @@ -399,47 +399,6 @@ subroutine li_advection_thickness_tracers(& if ( (trim(config_thickness_advection) .eq. 'fo') .or. & (trim(config_thickness_advection) .eq. 'fct') ) then - ! 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, 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") if (config_print_thickness_advection_info) then call mpas_log_write('Using first-order upwind for thickness advection') @@ -600,6 +559,48 @@ subroutine li_advection_thickness_tracers(& 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 From 49ce8b2d6c3fad9d5a9d5b45588197584ea6f043 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 17 Apr 2024 14:41:03 -0600 Subject: [PATCH 29/32] Clean up after rebase Fix a few minor issues that were not caught during the rebase. --- .../src/mode_forward/mpas_li_calving.F | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) 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 10b4726c49ce..3f140f5e5fa4 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 @@ -332,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 @@ -341,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) ! Parse grounded vs floating calving for budgets block => domain % blocklist @@ -410,7 +410,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) @@ -1001,7 +1001,6 @@ subroutine thickness_calving(domain, calvingFraction, err) thickness(:) = thickness(:) - dCalvingThickness(:) call update_calving_budget(geometryPool, dCalvingThickness) - call remove_small_islands(meshPool, geometryPool) block => block % next enddo @@ -1084,7 +1083,6 @@ subroutine floating_calving(domain, calvingFraction, err) thickness(:) = thickness(:) - dCalvingThickness(:) call update_calving_budget(geometryPool, dCalvingThickness) - call remove_small_islands(meshPool, geometryPool) deallocate(dCalvingThickness) block => block % next @@ -1296,7 +1294,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) @@ -1411,7 +1409,6 @@ subroutine topographic_calving(domain, calvingFraction, err) thickness(:) = thickness(:) - dCalvingThickness(:) call update_calving_budget(geometryPool, dCalvingThickness) - call remove_small_islands(meshPool, geometryPool) deallocate(dCalvingThickness) @@ -1630,7 +1627,6 @@ subroutine eigencalving(domain, err) call update_calving_budget(domain) call remove_icebergs(domain) - call remove_small_islands(meshPool, geometryPool) deallocate(dCalvingThickness) block => block % next @@ -1797,9 +1793,6 @@ subroutine specified_calving_velocity(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! TODO: global reduce & reporting on amount of calving generated in this step - call update_calving_budget(domain) - block => block % next deallocate(dCalvingThickness) @@ -2189,6 +2182,8 @@ subroutine von_Mises_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call remove_icebergs(domain) + deallocate(dCalvingThickness) block => block % next @@ -3514,6 +3509,10 @@ subroutine damage_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call remove_icebergs(domain) + + deallocate(dCalvingThickness) + block => block % next enddo @@ -4178,7 +4177,6 @@ subroutine mask_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - deallocate(dCalvingThickness) block => block % next enddo From 5a892242b4abbae2fcbe58e9e99021625fbbe41b Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 10 Jul 2024 13:39:06 -0700 Subject: [PATCH 30/32] Remove code that was inadvertently overwriting grounded and floating calving thickness Remove block of old code that was inadvertently overwriting grounded and floating calving thickness after they were calculated incrementally. --- .../src/mode_forward/mpas_li_calving.F | 23 ------------------- 1 file changed, 23 deletions(-) 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 3f140f5e5fa4..3c2fb62f02b2 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 @@ -352,29 +352,6 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - ! Parse grounded vs floating calving for budgets - block => domain % blocklist - do while (associated(block)) - call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) - call mpas_pool_get_array(geometryPool, 'groundedCalvingThickness', groundedCalvingThickness) - call mpas_pool_get_array(geometryPool, 'floatingCalvingThickness', floatingCalvingThickness) - call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) - call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) - - groundedCalvingThickness(:) = 0.0_RKIND - floatingCalvingThickness(:) = 0.0_RKIND - where (groundedMaskForMassBudget .eq. 1) - groundedCalvingThickness = calvingThickness - elsewhere (floatingMaskForMassBudget .eq. 1) - floatingCalvingThickness = calvingThickness - elsewhere - groundedCalvingThickness = 0.0_RKIND - floatingCalvingThickness = 0.0_RKIND - end where - - block => block % next - end do - ! Final operations after calving has been applied, including removal ! of small islands block => domain % blocklist From b0581d1893849b11a4c812e93229cf1b4af57691 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 10 Jul 2024 13:41:01 -0700 Subject: [PATCH 31/32] Update halos on incrementalCalvingThickness Replace the allocatable array dCalvingThickness with a Registry variable incrementalCalvingThickness, and update halos after calling li_apply_front_ablation_velocity and before subtracting from ice thickness. --- .../mpas-albany-landice/src/Registry.xml | 3 + .../src/mode_forward/mpas_li_calving.F | 215 +++++++++--------- 2 files changed, 105 insertions(+), 113 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 227e4d76b3a2..7523249ffe27 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1275,6 +1275,9 @@ is the value of that variable from the *previous* time level! + 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 3c2fb62f02b2..f21cf0c12d15 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 @@ -493,8 +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(:), allocatable:: & - dCalvingThickness ! incremental calving thickness + 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 @@ -548,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) @@ -640,8 +641,7 @@ subroutine li_restore_calving_front(domain, err) endif ! if preventing retreat ! Now check for marine regions that have advanced - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND ! loop over locally owned cells do iCell = 1, nCells @@ -656,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 - dCalvingThickness(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) @@ -664,9 +664,8 @@ subroutine li_restore_calving_front(domain, err) endif ! bedTopography < config_sea_level enddo ! iCell - call update_calving_budget(geometryPool, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) - deallocate(dCalvingThickness) block => block % next enddo @@ -709,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 @@ -810,8 +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(:), allocatable :: & - dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: & + incrementalCalvingThickness integer, pointer :: nCells integer :: iCell, iCellOnCell, iCellNeighbor @@ -879,8 +878,7 @@ subroutine thickness_calving(domain, calvingFraction, err) ! Initialize calvingThickness = 0.0_RKIND - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 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 @@ -970,14 +968,15 @@ 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) - dCalvingThickness = thickness * calvingFraction + incrementalCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - dCalvingThickness(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + + call update_calving_budget(geometryPool, incrementalCalvingThickness) - call update_calving_budget(geometryPool, dCalvingThickness) block => block % next enddo @@ -985,7 +984,6 @@ subroutine thickness_calving(domain, calvingFraction, err) call mpas_deallocate_scratch_field(activeForCalvingMaskField, .true.) call mpas_deallocate_scratch_field(inactiveMarginMaskField, .true.) call mpas_deallocate_scratch_field(oceanMaskField, .true.) - deallocate(dCalvingThickness) end subroutine thickness_calving @@ -1025,7 +1023,7 @@ 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(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask integer, dimension(:), pointer :: floatingMaskForMassBudget @@ -1047,21 +1045,19 @@ subroutine floating_calving(domain, calvingFraction, err) call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) call mpas_pool_get_array(meshPool, 'nCells', nCells) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND ! Note: The floating_ice mask does not include floating ! non-dynamic cells adjacent to grounded cells. where ( floatingMaskForMassBudget .eq. 1 ) - dCalvingThickness = thickness * calvingFraction + incrementalCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - dCalvingThickness(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) - deallocate(dCalvingThickness) block => block % next enddo @@ -1092,7 +1088,7 @@ subroutine remove_small_islands(domain, err) real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine) groundedCalvingThickness, & floatingCalvingThickness - real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: bedTopography integer, dimension(:), pointer :: cellMask, & @@ -1125,6 +1121,7 @@ subroutine remove_small_islands(domain, err) call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges) 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) @@ -1135,7 +1132,7 @@ subroutine remove_small_islands(domain, err) allocate(connectedCellsList((maxEdges+1)**2), & islandMask(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND islandMask(:) = 0 nIslandCellsLocal = 0 nIslandCellsGlobal = 0 @@ -1291,7 +1288,7 @@ subroutine remove_small_islands(domain, err) do iCell = 1, nCells if (li_mask_is_floating_ice(cellMask(iCell)) .and. seedMask(iCell) == 0) then calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) - dCalvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo @@ -1302,9 +1299,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, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) - deallocate(dCalvingThickness) end subroutine remove_small_islands @@ -1342,7 +1338,7 @@ 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(:), allocatable :: dCalvingThickness + 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 @@ -1368,26 +1364,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) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND where ( (li_mask_is_floating_ice(cellMask)) .and. (bedTopography < config_calving_topography + config_sea_level) ) - dCalvingThickness = thickness * calvingFraction + incrementalCalvingThickness = thickness * calvingFraction endwhere - calvingThicknessFromThreshold = calvingThickness + calvingThicknessFromThreshold = incrementalCalvingThickness ! === apply calving === - thickness(:) = thickness(:) - dCalvingThickness(:) - - call update_calving_budget(geometryPool, dCalvingThickness) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) - deallocate(dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next enddo @@ -1446,7 +1440,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(:), allocatable :: dCalvingThickness + 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 @@ -1504,6 +1498,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) @@ -1528,8 +1523,7 @@ subroutine eigencalving(domain, err) realArgs=(/minval(eigencalvingParameter), maxval(eigencalvingParameter)/)) endif - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + 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 @@ -1537,7 +1531,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, & - dCalvingThickness, calvingVelocity, applyToGrounded, & + incrementalCalvingThickness, calvingVelocity, applyToGrounded, & applyToFloating, applyToGroundingLine, domain, calvingCFLdt, dtCalvingCFLratio, & totalRatebasedCalvedVolume, totalRatebasedUncalvedVolume, err_tmp) err = ior(err, err_tmp) @@ -1549,7 +1543,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, dCalvingThickness, 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 & @@ -1568,11 +1562,11 @@ 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(:) - dCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) @@ -1592,12 +1586,12 @@ 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 - dCalvingThickness(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, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) @@ -1605,7 +1599,6 @@ subroutine eigencalving(domain, err) call update_calving_budget(domain) call remove_icebergs(domain) - deallocate(dCalvingThickness) block => block % next enddo @@ -1657,7 +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(:), allocatable :: dCalvingThickness ! Incremental calving thickness + 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 @@ -1703,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) @@ -1711,8 +1705,7 @@ subroutine specified_calving_velocity(domain, err) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND ! get parameter value if (trim(config_calving_specified_source) == 'const') then @@ -1731,17 +1724,19 @@ subroutine specified_calving_velocity(domain, err) endif ! Convert calvingVelocity to calvingThickness - call li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool, dCalvingThickness, 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(:) - dCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + 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) @@ -1760,19 +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 - dCalvingThickness(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, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) block => block % next - deallocate(dCalvingThickness) enddo end subroutine specified_calving_velocity @@ -1833,7 +1827,7 @@ subroutine von_Mises_calving(domain, err) floatingVonMisesThresholdStress, & groundedVonMisesThresholdStress, & surfaceSpeed - real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:,:), pointer :: flowParamA, & temperature, layerThickness real (kind=RKIND), pointer :: config_default_flowParamA @@ -1898,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) @@ -1946,8 +1941,7 @@ subroutine von_Mises_calving(domain, err) endif vonMisesStress(:) = 0.0_RKIND - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 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 @@ -2084,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, & - dCalvingThickness, calvingVelocity, applyToGrounded, & + incrementalCalvingThickness, calvingVelocity, applyToGrounded, & applyToFloating, applyToGroundingLine, domain, & calvingCFLdt, dtCalvingCFLratio, & totalRatebasedCalvedVolume, totalRatebasedUncalvedVolume, err_tmp) @@ -2096,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 @@ -2108,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 - dCalvingThickness(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 @@ -2119,15 +2113,15 @@ subroutine von_Mises_calving(domain, err) endif enddo if ( nGroundedNeighbors == 0 ) then - dCalvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) endif endif enddo endif ! === apply calving velocity before damage threshold calving === - thickness(:) = thickness(:) - dCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + 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) @@ -2139,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, dCalvingThickness, 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 & @@ -2153,15 +2147,14 @@ subroutine von_Mises_calving(domain, err) endif ! === apply calving velocity before damage threshold calving === - thickness(:) = thickness(:) - dCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + 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) call remove_icebergs(domain) - deallocate(dCalvingThickness) block => block % next enddo ! associated(block) @@ -2228,7 +2221,7 @@ subroutine ismip6_retreat(domain, err) calvingVelocity, thickness, & xvelmean, yvelmean, calvingThickness, & bedTopography - real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio integer, pointer :: nCells, timestepNumber @@ -2272,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) @@ -2290,8 +2284,7 @@ subroutine ismip6_retreat(domain, err) ! submergedArea used in runoff unit conversion allocate(submergedArea(nCells+1)) submergedArea(:) = 0.0_RKIND - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 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 @@ -2415,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, & - dCalvingThickness, calvingVelocity, applyToGrounded=.true., & + incrementalCalvingThickness, calvingVelocity, applyToGrounded=.true., & applyToFloating=.true., applyToGroundingLine=.false., & domain=domain, maxDt=calvingCFLdt, CFLratio=dtCalvingCFLratio, & totalAblatedVolume=totalRatebasedCalvedVolume, & @@ -2426,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(:) - dCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + 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) deallocate(submergedArea) - deallocate(dCalvingThickness) end subroutine ismip6_retreat !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -3370,7 +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(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:), pointer :: calvingVelocity real (kind=RKIND), pointer :: calvingCFLdt real (kind=RKIND), pointer :: dtCalvingCFLratio @@ -3423,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) @@ -3434,8 +3427,7 @@ subroutine damage_calving(domain, err) call mpas_pool_get_array(geometryPool, 'totalRatebasedCalvedVolume', totalRatebasedCalvedVolume) call mpas_pool_get_array(geometryPool, 'totalRatebasedUncalvedVolume', totalRatebasedUncalvedVolume) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND call calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask) @@ -3457,15 +3449,18 @@ 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, dCalvingThickness, 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(:) - dCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) + 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) @@ -3477,19 +3472,17 @@ 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 - dCalvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo - call update_calving_budget(geometryPool, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) call remove_icebergs(domain) - deallocate(dCalvingThickness) - block => block % next enddo @@ -3963,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, dCalvingThickness, err) + subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, domain, incrementalCalvingThickness, err) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- @@ -3979,7 +3972,7 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - real (kind=RKIND), dimension(:), intent(out) :: dCalvingThickness ! incremental calving thickness + real (kind=RKIND), dimension(:), intent(out) :: incrementalCalvingThickness ! incremental calving thickness !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -4020,7 +4013,7 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d growMask => growMaskField % array growMask(:) = 0 - dCalvingThickness(:) = 0.0_RKIND + 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) ) @@ -4041,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 - dCalvingThickness(jCell) = thickness(jCell) + incrementalCalvingThickness(jCell) = thickness(jCell) calvingThicknessFromThreshold(jCell) = calvingThicknessFromThreshold(jCell) + thickness(jCell) endif enddo calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) - dCalvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) endif ! if cell is calving margin enddo ! cell loop @@ -4102,7 +4095,7 @@ subroutine mask_calving(domain, err) integer, pointer :: nCells, nCellsSolve integer :: iCell integer :: localMaskCellCount, globalMaskCellCount - real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness integer :: err_tmp err = 0 @@ -4120,14 +4113,14 @@ 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) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND + incrementalCalvingThickness(:) = 0.0_RKIND ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) @@ -4139,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/)) - dCalvingThickness(iCell) = thickness(iCell) + incrementalCalvingThickness(iCell) = thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) thickness(iCell) = 0.0_RKIND if (iCell <= nCellsSolve) localMaskCellCount = localMaskCellCount + 1 @@ -4149,12 +4142,11 @@ 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, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - deallocate(dCalvingThickness) block => block % next enddo @@ -4269,7 +4261,7 @@ subroutine li_remove_icebergs(domain) real (kind=RKIND), dimension(:), pointer :: calvingThickness, & ! thickness of ice that calves (computed in this subroutine) groundedCalvingThickness, & floatingCalvingThickness - real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness + real (kind=RKIND), dimension(:), pointer :: incrementalCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask, & groundedMaskForMassBudget, & @@ -4310,9 +4302,6 @@ subroutine li_remove_icebergs(domain) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND - seedMask(:) = 0 growMask(:) = 0 @@ -4373,6 +4362,7 @@ 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) @@ -4387,14 +4377,14 @@ subroutine li_remove_icebergs(domain) do iCell = 1, nCellsSolve if (seedMask(iCell) == 0 .and. li_mask_is_floating_ice(cellMask(iCell))) then calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) ! remove any remaining ice here - dCalvingThickness(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, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next end do @@ -4413,7 +4403,6 @@ subroutine li_remove_icebergs(domain) call mpas_log_write("Iceberg-detection flood-fill complete. Removed $i iceberg cells.", intArgs=(/globalIcebergCellCount/)) call mpas_timer_stop("iceberg detection") - deallocate(dCalvingThickness) end subroutine remove_icebergs >>>>>>> a5bee82292 (Update calving budget incrementally) @@ -4581,12 +4570,12 @@ end subroutine li_flood_fill !> is applied, but before masks are updated which often happens multiple times !> in a timestep. !----------------------------------------------------------------------- - subroutine update_calving_budget(geometryPool, dCalvingThickness) + subroutine update_calving_budget(geometryPool, incrementalCalvingThickness) !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- type (mpas_pool_type), pointer, intent(inout) :: geometryPool - real (kind=RKIND), dimension(:), intent(inout) :: dCalvingThickness + real (kind=RKIND), dimension(:), intent(inout) :: incrementalCalvingThickness !----------------------------------------------------------------- ! local variables @@ -4604,13 +4593,13 @@ subroutine update_calving_budget(geometryPool, dCalvingThickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) where (groundedMaskForMassBudget .eq. 1) - groundedCalvingThickness = groundedCalvingThickness + dCalvingThickness + groundedCalvingThickness = groundedCalvingThickness + incrementalCalvingThickness elsewhere (floatingMaskForMassBudget .eq. 1) - floatingCalvingThickness = floatingCalvingThickness + dCalvingThickness + floatingCalvingThickness = floatingCalvingThickness + incrementalCalvingThickness end where - calvingThickness(:) = calvingThickness(:) + dCalvingThickness(:) - dCalvingThickness(:) = 0.0_RKIND + calvingThickness(:) = calvingThickness(:) + incrementalCalvingThickness(:) + incrementalCalvingThickness(:) = 0.0_RKIND end subroutine update_calving_budget From 6597665558da4527027112bdaf14d048aca2eb77 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 25 Mar 2025 13:07:16 -0700 Subject: [PATCH 32/32] Clean up after rebase Clean up some issues that remained after rebasing onto develop. --- .../src/mode_forward/mpas_li_advection.F | 3 +++ .../src/mode_forward/mpas_li_calving.F | 11 +++++------ .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 4 ++-- 3 files changed, 10 insertions(+), 8 deletions(-) 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 da147eb5b60e..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 @@ -368,6 +368,9 @@ subroutine li_advection_thickness_tracers(& 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.) + call mpas_pool_get_field(geometryPool, 'edgeMaskTemporary', edgeMaskTemporaryField) call mpas_allocate_scratch_field(edgeMaskTemporaryField, .true.) 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 f21cf0c12d15..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 @@ -1119,6 +1119,7 @@ subroutine remove_small_islands(domain, err) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges) + 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) @@ -1596,8 +1597,7 @@ subroutine eigencalving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call update_calving_budget(domain) - call remove_icebergs(domain) + call li_remove_icebergs(domain) block => block % next enddo @@ -2153,7 +2153,7 @@ subroutine von_Mises_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call remove_icebergs(domain) + call li_remove_icebergs(domain) block => block % next @@ -3481,7 +3481,7 @@ subroutine damage_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call remove_icebergs(domain) + call li_remove_icebergs(domain) block => block % next @@ -4403,8 +4403,7 @@ subroutine li_remove_icebergs(domain) call mpas_log_write("Iceberg-detection flood-fill complete. Removed $i iceberg cells.", intArgs=(/globalIcebergCellCount/)) call mpas_timer_stop("iceberg detection") - end subroutine remove_icebergs ->>>>>>> a5bee82292 (Update calving budget incrementally) + end subroutine li_remove_icebergs !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! 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 6765cfbb4c59..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 !----------------------------------------------------------------- @@ -366,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 ***