From dc058a21be3deec2afc99c3d3aece4b30dba965c Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 13:58:39 -0600 Subject: [PATCH 01/33] 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 | 30 +++++----- .../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 | 60 ++++++++++++++++++- 12 files changed, 97 insertions(+), 31 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 3aa73768f04f..eab540db34a1 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1201,6 +1201,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 56b018cdd399..e02f81c89000 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 @@ -367,7 +369,7 @@ subroutine li_advection_thickness_tracers(& call li_calculate_layerThickness(meshPool, thickness, layerThickness) ! update masks - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) ! save old copycellMask for determining cells changing from grounded to floating and vice versa @@ -543,7 +545,7 @@ subroutine li_advection_thickness_tracers(& ! Update the thickness and cellMask before applying the mass balance. ! The update is needed because the SMB and BMB depend on whether ice is present. thickness = sum(layerThickness, 1) - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) @@ -627,7 +629,7 @@ subroutine li_advection_thickness_tracers(& enddo endif - 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 de2a523e9339..5bf46f0644a5 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 @@ -199,7 +199,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) @@ -574,7 +574,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 @@ -645,7 +645,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 @@ -1505,7 +1505,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. @@ -1529,7 +1529,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 @@ -1687,7 +1687,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. @@ -1711,7 +1711,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 @@ -2119,7 +2119,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 @@ -2389,7 +2389,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) @@ -3420,7 +3420,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. @@ -3436,7 +3436,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 @@ -3828,7 +3828,7 @@ subroutine li_finalize_damage_after_advection(domain, err) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) ! 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) if (config_preserve_damage) then @@ -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 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 5a92ecba377d..a7a236492fcc 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 60cd6e5ca795..4be457cb50c9 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 @@ -167,7 +167,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 6ca13f8deb0d..bb76e6a3ab36 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 @@ -1099,6 +1099,7 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & + domain, & err_tmp, & advectTracersIn = .true.) @@ -1118,6 +1119,7 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & + domain, & err_tmp, & advectTracersIn = .false.) @@ -1197,7 +1199,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 673ed8186e47..a14be4b43806 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -246,7 +246,7 @@ end subroutine li_calculate_mask_init ! !----------------------------------------------------------------------- - subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) + subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) !----------------------------------------------------------------- ! @@ -265,6 +265,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 @@ -284,7 +286,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 @@ -306,6 +309,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') @@ -328,6 +332,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) @@ -357,9 +363,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 @@ -473,6 +484,51 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) endif 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 093f93eefa2856513bf1c0ca94a3b127ba8f2296 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 14:47:52 -0600 Subject: [PATCH 02/33] 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 eab540db34a1..faa6544750af 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1224,6 +1224,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 e02f81c89000..04ac17da3656 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 @@ -203,7 +204,9 @@ subroutine li_advection_thickness_tracers(& integer, dimension(:), pointer :: & cellMask, & ! integer bitmask for cells edgeMask, & ! integer bitmask for edges - thermalCellMask + thermalCellMask, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget integer, dimension(:,:), pointer :: cellsOnEdge @@ -282,6 +285,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) @@ -294,7 +298,8 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) 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) @@ -597,10 +602,13 @@ subroutine li_advection_thickness_tracers(& sfcMassBal, & sfcMassBalApplied, & groundedSfcMassBalApplied, & + floatingSfcMassBalApplied, & basalMassBal, & basalMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & surfaceTracers, & basalTracers, & layerThickness, & @@ -846,10 +854,13 @@ subroutine apply_mass_balance(& sfcMassBal, & sfcMassBalApplied, & groundedSfcMassBalApplied, & + floatingSfcMassBalApplied, & basalMassBal, & basalMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + groundedMaskForMassBudget, & + floatingMaskForMassBudget, & surfaceTracers, & basalTracers, & layerThickness, & @@ -863,7 +874,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) @@ -892,7 +905,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) @@ -1048,10 +1062,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 4a130e95b663ed3cbdd1e9e4626fd08fcf15b20b Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 15:04:34 -0600 Subject: [PATCH 03/33] 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 faa6544750af..10b7bffe85d9 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1259,6 +1259,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) @@ -338,6 +342,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. @@ -367,6 +373,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 915be80a58dc72dfc1c299218df349c59146c666 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 15:26:31 -0600 Subject: [PATCH 04/33] 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 10b7bffe85d9..21ce81608c1d 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1420,6 +1420,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 @@ -373,20 +396,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 a7a236492fcc..9d74defa2201 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 11f4e1ea6d0a63f1b77c52e4ad0c1e782ad75043 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 16:08:39 -0600 Subject: [PATCH 05/33] 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 37727a42dfa6c5675ffaeda64fbc8e6fa5e6f5a2 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 17:02:51 -0600 Subject: [PATCH 06/33] 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 a14be4b43806..2d96c1c76699 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -364,6 +364,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 82dc2662f104bd0c1eae3f0e4dd7d58b3b4cf2c3 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 16 May 2022 21:17:54 -0600 Subject: [PATCH 07/33] 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 14c9d6996baa..21904fce4c60 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 @@ -246,6 +246,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 @@ -365,8 +383,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. @@ -667,6 +683,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)) @@ -971,6 +988,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 @@ -1045,6 +1068,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 @@ -1070,11 +1099,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 @@ -1104,6 +1142,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) @@ -1169,6 +1211,7 @@ subroutine remove_small_islands(domain, err) exit endif enddo +<<<<<<< HEAD endif enddo @@ -1218,6 +1261,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 @@ -1267,6 +1334,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***") @@ -1275,6 +1343,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 @@ -1352,6 +1444,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 @@ -1532,6 +1630,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) @@ -1556,7 +1655,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) @@ -1581,6 +1680,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 @@ -1714,7 +1818,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) @@ -1738,7 +1842,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) @@ -1763,6 +1867,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 @@ -2146,11 +2255,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) @@ -2416,7 +2530,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) @@ -3447,7 +3561,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) @@ -3463,7 +3577,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) @@ -3488,7 +3602,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 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 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 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 2f48a21f50bdb2d6d63cc09834019820ab47aabf Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 17 May 2022 17:55:27 -0600 Subject: [PATCH 08/33] 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 | 50 ++++++++++++++----- 2 files changed, 43 insertions(+), 15 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 21ce81608c1d..3518fc09c394 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1448,10 +1448,14 @@ 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 04ac17da3656..b2600c12c933 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 @@ -196,7 +196,9 @@ subroutine li_advection_thickness_tracers(& type (field1DInteger), pointer :: & cellMaskTemporaryField, & ! scratch field containing old values of cellMask - 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 @@ -360,8 +362,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(scratchPool, 'iceCellMask', thermalCellMaskField) call mpas_allocate_scratch_field(thermalCellMaskField, .true.) @@ -377,9 +381,10 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! save old copycellMask for determining cells changing from grounded to floating and vice versa - cellMaskTemporaryField % array(:) = cellMask(:) 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 @@ -640,6 +645,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 @@ -744,6 +759,8 @@ subroutine li_advection_thickness_tracers(& call mpas_deallocate_scratch_field(surfaceTracersField, .true.) call mpas_deallocate_scratch_field(cellMaskTemporaryField, .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 @@ -2125,16 +2142,23 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers 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 @@ -2159,12 +2183,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 be7bdf148ac594534962cb5a5f430f1e52f7b778 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 May 2022 17:21:31 -0600 Subject: [PATCH 09/33] 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 | 70 +++++++++---------- 1 file changed, 32 insertions(+), 38 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 b2600c12c933..9f8d847959d3 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 @@ -393,6 +393,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/)) @@ -655,44 +687,6 @@ subroutine li_advection_thickness_tracers(& 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 - - ! 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 61282498608901977b88ce7196d2c9a293b561c9 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 May 2022 18:38:19 -0600 Subject: [PATCH 10/33] 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 | 83 +------------------ .../mpas_li_time_integration_fe_rk.F | 6 ++ .../src/shared/mpas_li_mask.F | 28 ++++++- 3 files changed, 35 insertions(+), 82 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 9f8d847959d3..c5abcbba500c 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 !-------------------------------------------------------------------- ! @@ -196,9 +195,7 @@ subroutine li_advection_thickness_tracers(& type (field1DInteger), pointer :: & cellMaskTemporaryField, & ! scratch field containing old values of cellMask - thermalCellMaskField, & - groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget - floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget + thermalCellMaskField ! Allocatable arrays need for flux-corrected transport advection real (kind=RKIND), dimension(:,:,:), allocatable :: tend @@ -362,11 +359,6 @@ subroutine li_advection_thickness_tracers(& call mpas_allocate_scratch_field(basalTracersField, .true.) basalTracers => basalTracersField % array - call mpas_pool_get_field(geometryPool, 'groundedMaskForMassBudgetTemporary', groundedMaskForMassBudgetTemporaryField) - call mpas_pool_get_field(geometryPool, 'floatingMaskForMassBudgetTemporary', floatingMaskForMassBudgetTemporaryField) - call mpas_allocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.) - call mpas_allocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.) - call mpas_pool_get_field(scratchPool, 'iceCellMask', thermalCellMaskField) call mpas_allocate_scratch_field(thermalCellMaskField, .true.) thermalCellMask => thermalCellMaskField % array @@ -382,10 +374,6 @@ subroutine li_advection_thickness_tracers(& err = ior(err, err_tmp) 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 !----------------------------------------------------------------- @@ -589,8 +577,6 @@ subroutine li_advection_thickness_tracers(& thickness = sum(layerThickness, 1) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - - !----------------------------------------------------------------- ! Add the surface and basal mass balance to the layer thickness !----------------------------------------------------------------- @@ -677,16 +663,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. @@ -753,8 +729,6 @@ subroutine li_advection_thickness_tracers(& call mpas_deallocate_scratch_field(surfaceTracersField, .true.) call mpas_deallocate_scratch_field(cellMaskTemporaryField, .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 +2110,6 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers end subroutine vertical_remap - subroutine grounded_to_floating(groundedMaskForMassBudgetOrig, & - floatingMaskForMassBudgetOrig, & - groundedMaskForMassBudgetNew, & - floatingMaskForMassBudgetNew, & - thicknessNew, & - groundedToFloatingThickness, & - nCells) - !----------------------------------------------------------------- - ! input variables - !----------------------------------------------------------------- - integer, dimension(:), intent(in) :: & - groundedMaskForMassBudgetOrig, & !< Input: mask for cells before advection - floatingMaskForMassBudgetOrig - - integer, dimension(:), intent(in) :: & - groundedMaskForMassBudgetNew, & !< Input: mask for cells after advection - floatingMaskForMassBudgetNew - real(kind=RKIND), dimension(:), intent(in) :: & - thicknessNew !< Input: ice thickness after advection - - integer, pointer, intent(in) :: & - nCells - !----------------------------------------------------------------- - ! input/output variables - !----------------------------------------------------------------- - - !----------------------------------------------------------------- - ! output variables - !----------------------------------------------------------------- - real(kind=RKIND), dimension(:), intent(out) :: & - groundedToFloatingThickness !< Input: ice thickness - - !----------------------------------------------------------------- - ! local variables - !----------------------------------------------------------------- - integer :: iCell - - groundedToFloatingThickness(:) = 0.0_RKIND - - do iCell = 1, nCells - ! find locations that had grounded ice before but now have floating ice - if ( (groundedMaskForMassBudgetOrig(iCell) .eq. 1) .and. & - (floatingMaskForMassBudgetNew(iCell) .eq. 1) ) then - groundedToFloatingThickness(iCell) = thicknessNew(iCell) - ! find locations that had floating ice before but now have grounded ice - else if ( (floatingMaskForMassBudgetOrig(iCell) .eq. 1) .and. & - (groundedMaskForMassBudgetNew(iCell) .eq. 1) ) then - groundedToFloatingThickness(iCell) = -1.0_RKIND * thicknessNew(iCell) - endif - enddo - - end subroutine 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 bb76e6a3ab36..dad87f4d62b1 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 @@ -655,6 +655,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 @@ -748,8 +749,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 2d96c1c76699..19cf9c1ef6ee 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -285,7 +285,7 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) !----------------------------------------------------------------- integer, pointer :: nCells, nVertices, nEdges, vertexDegree integer, pointer :: nVertInterfaces - real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography + real(KIND=RKIND), dimension(:), pointer :: thickness, bedTopography, groundedToFloatingThickness integer, dimension(:), pointer :: nEdgesOnCell, cellMask, vertexMask, edgeMask, & groundedMaskForMassBudget, floatingMaskForMassBudget integer, dimension(:,:), pointer :: cellsOnCell, cellsOnVertex, cellsOnEdge, dirichletVelocityMask @@ -310,6 +310,9 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) integer :: iCellNeighbor logical :: openOceanNeighbor logical :: updateMassBudgetMasks=.true. + type (field1DInteger), pointer :: & ! used to calculate groundedToFloatingThickness + groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget + floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget call mpas_timer_start('calculate mask') @@ -332,9 +335,15 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) call mpas_pool_get_array(geometryPool, 'groundedMaskForMassBudget', groundedMaskForMassBudget) call mpas_pool_get_array(geometryPool, 'floatingMaskForMassBudget', floatingMaskForMassBudget) + call mpas_pool_get_field(geometryPool, 'groundedMaskForMassBudgetTemporary', groundedMaskForMassBudgetTemporaryField) + call mpas_pool_get_field(geometryPool, 'floatingMaskForMassBudgetTemporary', floatingMaskForMassBudgetTemporaryField) + call mpas_allocate_scratch_field(groundedMaskForMassBudgetTemporaryField, .true.) + call mpas_allocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.) + call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) @@ -360,6 +369,11 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) cellMask = ior(cellMask, li_mask_ValueIce) end where + + ! save old mass budget masks for determining cells changing from grounded + ! to floating and vice versa + groundedMaskForMassBudgetTemporaryField % array(:) = groundedMaskForMassBudget(:) + floatingMaskForMassBudgetTemporaryField % array(:) = floatingMaskForMassBudget(:) ! Is it floating? (ice thickness equal to floatation is considered floating) ! For now floating ice and grounded ice are mutually exclusive. ! This may change if a grounding line parameterization is added. @@ -518,6 +532,18 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) endif enddo + ! Update groundedToFloatingThickness based on new masks + do iCell = 1, nCells + ! find locations that had grounded ice before but now have floating ice + if ( (groundedMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (floatingMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) + thickness(iCell) + ! find locations that had floating ice before but now have grounded ice + else if ( (floatingMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (groundedMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) - thickness(iCell) + endif + enddo ! Now use flood fill to identify all floating ice connected to the open ! ocean. Subglacial lakes are any floating ice that are not in the ! resulting flood fill mask. Commented code below does not yet work From 4f566f0ad571d330564b701235a29d0d9c4d8c70 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 15 Jul 2022 15:39:48 -0600 Subject: [PATCH 11/33] Add halo updates after calling li_calculate_mask Add halo updates for mass budget masks after calling li_calculate_mask. This did not help close the budget, but in case it is correct I don't want to lose the work. Committing and then reverting. --- .../src/mode_forward/mpas_li_advection.F | 12 ++ .../src/mode_forward/mpas_li_bedtopo.F | 8 + .../src/mode_forward/mpas_li_calving.F | 152 ++++++++---------- .../src/mode_forward/mpas_li_core.F | 16 +- .../src/mode_forward/mpas_li_iceshelf_melt.F | 10 +- .../mode_forward/mpas_li_subglacial_hydro.F | 4 + .../src/mode_forward/mpas_li_thermal.F | 5 +- .../mpas_li_time_integration_fe_rk.F | 4 + .../src/mode_forward/mpas_li_velocity.F | 4 + .../mode_forward/mpas_li_velocity_external.F | 5 +- .../src/shared/mpas_li_mask.F | 46 +++--- 11 files changed, 150 insertions(+), 116 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 c5abcbba500c..ffa72d5c73a6 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 @@ -374,6 +374,10 @@ subroutine li_advection_thickness_tracers(& err = ior(err, err_tmp) layerThicknessEdgeFlux(:,:) = 0.0_RKIND + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers !----------------------------------------------------------------- @@ -577,6 +581,10 @@ subroutine li_advection_thickness_tracers(& thickness = sum(layerThickness, 1) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") !----------------------------------------------------------------- ! Add the surface and basal mass balance to the layer thickness !----------------------------------------------------------------- @@ -663,6 +671,10 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! 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. 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 be5ea3281e8c..4f9ec5d6915c 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 @@ -253,6 +253,10 @@ subroutine li_bedtopo_solve(domain, err) call li_update_geometry(geometryPool) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") block => block % next end do @@ -764,6 +768,10 @@ subroutine slmodel_solve(slmTimeStep, domain) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! deallocate memory deallocate(globalArrayThickness) deallocate(gatheredArrayThickness) 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 21904fce4c60..e2835597f8a0 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 @@ -207,6 +207,10 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) err = ior(err, err_tmp) call li_update_geometry(geometryPool) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! based on the calving timescale, set the fraction of ice that calves if (config_calving_timescale > 0.0_RKIND) then calvingFraction = min(deltat/config_calving_timescale, 1.0_RKIND) @@ -412,7 +416,10 @@ 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) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") block => block % next end do if (present(solveVeloAfterCalving)) then @@ -622,6 +629,10 @@ subroutine li_restore_calving_front(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! initialize restoreThickness = 0.0_RKIND @@ -694,6 +705,10 @@ subroutine li_restore_calving_front(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call li_update_geometry(geometryPool) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") block => block % next end do @@ -988,12 +1003,9 @@ 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 @@ -1068,12 +1080,9 @@ 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 @@ -1099,11 +1108,7 @@ 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, & @@ -1211,7 +1216,6 @@ subroutine remove_small_islands(domain, err) exit endif enddo -<<<<<<< HEAD endif enddo @@ -1261,30 +1265,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 @@ -1334,7 +1314,6 @@ 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***") @@ -1343,30 +1322,6 @@ 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 @@ -1444,12 +1399,9 @@ 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 @@ -1636,6 +1588,10 @@ subroutine eigencalving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. !where ((li_mask_is_floating_ice(cellMask) .and. li_mask_is_dynamic_ice(cellMask) .and. & @@ -1659,7 +1615,10 @@ subroutine eigencalving(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! 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 @@ -1680,11 +1639,8 @@ 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 @@ -1823,6 +1779,10 @@ subroutine specified_calving_velocity(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. !where ((li_mask_is_floating_ice(cellMask) .and. li_mask_is_dynamic_ice(cellMask) .and. & @@ -1847,6 +1807,10 @@ subroutine specified_calving_velocity(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! 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 @@ -1867,11 +1831,8 @@ 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 @@ -2260,11 +2221,10 @@ subroutine von_Mises_calving(domain, err) 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) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") block => block % next enddo ! associated(block) @@ -2535,6 +2495,12 @@ subroutine ismip6_retreat(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") + call remove_small_islands(meshPool, geometryPool) + deallocate(submergedArea) end subroutine ismip6_retreat @@ -3565,7 +3531,10 @@ subroutine damage_calving(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. ! This criteria below only remove too-thin ice at the new calving front, ! meaning just one 'row' of cells per timestep. This could be expanded to continue @@ -3581,7 +3550,10 @@ subroutine damage_calving(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! 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 @@ -3602,12 +3574,9 @@ 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 @@ -3978,7 +3947,10 @@ subroutine li_finalize_damage_after_advection(domain, err) ! make sure masks are up to date. May not be necessary, but safer to do anyway. call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") if (config_preserve_damage) then do iCell = 1, nCells if (damage(iCell) < damageMax(iCell)) then @@ -4266,6 +4238,12 @@ subroutine mask_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") + call remove_small_islands(meshPool, geometryPool) + block => block % next enddo @@ -4440,7 +4418,10 @@ subroutine remove_icebergs(domain) ! do anyway. call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") call mpas_log_write("Iceberg-detection flood-fill: updated masks.") newMaskCountLocal = 0 do iCell = 1, nCellsSolve @@ -4588,7 +4569,10 @@ subroutine li_flood_fill(seedMask, growMask, domain) ! do anyway. call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") !call mpas_log_write("Flood-fill: updated masks.") newMaskCountLocal = sum(seedMask) !call mpas_log_write("Initialized $i cells to local mask", intArgs=(/newMaskCountLocal/)) 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..4e16df1c1c34 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 @@ -287,7 +287,7 @@ function li_core_init(domain, startTimeStamp) result(err) call mpas_pool_get_array(meshPool, 'cellProcID', cellProcID) cellProcID(:) = domain % dminfo % my_proc_id - call landice_init_block(block, domain % dminfo, err_tmp) + call landice_init_block(block, domain, domain % dminfo, err_tmp) err = ior(err, err_tmp) block => block % next @@ -804,6 +804,10 @@ function li_core_initial_solve(domain) result(err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") call li_update_geometry(geometryPool) block => block % next @@ -924,7 +928,7 @@ end function li_core_initial_solve !> This routine initializes blocks for the land ice core. ! !----------------------------------------------------------------------- - subroutine landice_init_block(block, dminfo, err) + subroutine landice_init_block(block, domain, dminfo, err) use mpas_derived_types use mpas_pool_routines @@ -942,14 +946,13 @@ subroutine landice_init_block(block, dminfo, err) ! !----------------------------------------------------------------- type (dm_info), intent(in) :: dminfo !< Input: Domain info - !----------------------------------------------------------------- ! ! input/output variables ! !----------------------------------------------------------------- type (block_type), intent(inout) :: block !< Input/output: Block object - + type (domain_type), intent(inout) :: domain !< Input/output: Domain !----------------------------------------------------------------- ! ! output variables @@ -1039,7 +1042,10 @@ subroutine landice_init_block(block, dminfo, err) ! Mask init identifies initial ice extent call li_calculate_mask_init(geometryPool, err=err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! Initialize thicknessOld variable used to calculate dHdt call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'thicknessOld', thicknessOld) 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 9d74defa2201..b830d41179ce 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 @@ -175,7 +175,10 @@ subroutine li_face_melt_grounded_ice(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") block => block % next enddo ! associated(block) @@ -414,7 +417,10 @@ subroutine li_basal_melt_floating_ice(domain, err) ! calculate masks - so we know where the ice is floating call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! calculate a mask to identify ice that is thick enough to be thermally active do iCell = 1, nCellsSolve if (thickness(iCell) > config_thermal_thickness) then 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 4be457cb50c9..7364966826e9 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,6 +172,10 @@ subroutine li_SGH_init(domain, err) call calc_hydro_mask(domain) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! remove invalid values - not necessary on restart, but shouldn't hurt call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness) waterThickness = max(0.0_RKIND, waterThickness) 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 7fec0b821932..ddbb019551c0 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 @@ -860,7 +860,10 @@ subroutine li_thermal_solver(domain, err) ! calculate masks - so we know where the ice is floating call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! calculate a mask to identify ice that is thick enough to be thermally active do iCell = 1, nCellsSolve if (thickness(iCell) > config_thermal_thickness) then 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 dad87f4d62b1..2249f7f2c141 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 @@ -1206,6 +1206,10 @@ subroutine advection_solver(domain, err) endif call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") 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 acffc35328ae..df9a46bb2c7d 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 @@ -334,6 +334,10 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") 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 29f7a6fab809..9289a480beed 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 @@ -876,7 +876,10 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! code is currently organized, that would be a circular dependency.a call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) - + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') + call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') + call mpas_timer_stop("halo updates") ! 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) ) lowerSurface = config_sea_level - thickness * (config_ice_density / config_ocean_density) 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 19cf9c1ef6ee..7bef9d6f970e 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -532,29 +532,29 @@ 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 - ! 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 + ! Update groundedToFloatingThickness based on new masks + do iCell = 1, nCells + ! find locations that had grounded ice before but now have floating ice + if ( (groundedMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (floatingMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) + thickness(iCell) + ! find locations that had floating ice before but now have grounded ice + else if ( (floatingMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (groundedMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) - thickness(iCell) + endif + enddo + ! Now use flood fill to identify all floating ice connected to the open + ! ocean. Subglacial lakes are any floating ice that are not in the + ! resulting flood fill mask. Commented code below does not yet work + ! because flood fill needs to be moved from li_calving to a shared + ! routine. + ! call li_flood_fill(seedMask=seedMask, growMask=li_mask_is_floating_ice(cellMask), domain) + + !where (li_mask_is_floating(cellMask) and seedMask .eq. 0) + ! groundedMaskForBudget = 1 + ! gloatingMaskForMassBudget = 0 + !end where endif !call mpas_timer_stop('calculate mask cell') From 93309abd984ed159d271958099db036fb99fcb38 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 15 Jul 2022 15:40:55 -0600 Subject: [PATCH 12/33] Revert "Add halo updates after calling li_calculate_mask" This reverts commit ff97c4eb3114dcd0f4fe14f099d8be2a67d0ac43. The changes in that commit did not help close budgets, but they may be helpful at some point. --- .../src/mode_forward/mpas_li_advection.F | 12 --- .../src/mode_forward/mpas_li_bedtopo.F | 8 -- .../src/mode_forward/mpas_li_calving.F | 74 +++---------------- .../src/mode_forward/mpas_li_core.F | 16 ++-- .../src/mode_forward/mpas_li_iceshelf_melt.F | 10 +-- .../mode_forward/mpas_li_subglacial_hydro.F | 4 - .../src/mode_forward/mpas_li_thermal.F | 5 +- .../mpas_li_time_integration_fe_rk.F | 4 - .../src/mode_forward/mpas_li_velocity.F | 4 - .../mode_forward/mpas_li_velocity_external.F | 5 +- .../src/shared/mpas_li_mask.F | 46 ++++++------ 11 files changed, 41 insertions(+), 147 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 ffa72d5c73a6..c5abcbba500c 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 @@ -374,10 +374,6 @@ subroutine li_advection_thickness_tracers(& err = ior(err, err_tmp) layerThicknessEdgeFlux(:,:) = 0.0_RKIND - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") !----------------------------------------------------------------- ! Horizontal transport of thickness and tracers !----------------------------------------------------------------- @@ -581,10 +577,6 @@ subroutine li_advection_thickness_tracers(& thickness = sum(layerThickness, 1) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") !----------------------------------------------------------------- ! Add the surface and basal mass balance to the layer thickness !----------------------------------------------------------------- @@ -671,10 +663,6 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! 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. 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 4f9ec5d6915c..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 @@ -253,10 +253,6 @@ subroutine li_bedtopo_solve(domain, err) call li_update_geometry(geometryPool) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") block => block % next end do @@ -768,10 +764,6 @@ subroutine slmodel_solve(slmTimeStep, domain) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! deallocate memory deallocate(globalArrayThickness) deallocate(gatheredArrayThickness) 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 e2835597f8a0..7164ab211e41 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 @@ -207,10 +207,6 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving) err = ior(err, err_tmp) call li_update_geometry(geometryPool) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! based on the calving timescale, set the fraction of ice that calves if (config_calving_timescale > 0.0_RKIND) then calvingFraction = min(deltat/config_calving_timescale, 1.0_RKIND) @@ -416,10 +412,7 @@ 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) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + block => block % next end do if (present(solveVeloAfterCalving)) then @@ -629,10 +622,6 @@ subroutine li_restore_calving_front(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! initialize restoreThickness = 0.0_RKIND @@ -705,10 +694,6 @@ subroutine li_restore_calving_front(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) call li_update_geometry(geometryPool) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") block => block % next end do @@ -1588,10 +1573,6 @@ subroutine eigencalving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. !where ((li_mask_is_floating_ice(cellMask) .and. li_mask_is_dynamic_ice(cellMask) .and. & @@ -1615,10 +1596,7 @@ subroutine eigencalving(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! 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 @@ -1779,10 +1757,6 @@ subroutine specified_calving_velocity(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. !where ((li_mask_is_floating_ice(cellMask) .and. li_mask_is_dynamic_ice(cellMask) .and. & @@ -1807,10 +1781,6 @@ subroutine specified_calving_velocity(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! 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 @@ -2221,10 +2191,8 @@ subroutine von_Mises_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + call remove_small_islands(meshPool, geometryPool) + block => block % next enddo ! associated(block) @@ -2494,11 +2462,6 @@ subroutine ismip6_retreat(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") call remove_small_islands(meshPool, geometryPool) deallocate(submergedArea) @@ -3531,10 +3494,7 @@ subroutine damage_calving(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! Now also remove thin floating, dynamic ice (based on chosen thickness threshold) after mask is updated. ! This criteria below only remove too-thin ice at the new calving front, ! meaning just one 'row' of cells per timestep. This could be expanded to continue @@ -3550,10 +3510,7 @@ subroutine damage_calving(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! 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 @@ -3947,10 +3904,7 @@ subroutine li_finalize_damage_after_advection(domain, err) ! make sure masks are up to date. May not be necessary, but safer to do anyway. call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + if (config_preserve_damage) then do iCell = 1, nCells if (damage(iCell) < damageMax(iCell)) then @@ -4238,10 +4192,6 @@ subroutine mask_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") call remove_small_islands(meshPool, geometryPool) block => block % next @@ -4418,10 +4368,7 @@ subroutine remove_icebergs(domain) ! do anyway. call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + call mpas_log_write("Iceberg-detection flood-fill: updated masks.") newMaskCountLocal = 0 do iCell = 1, nCellsSolve @@ -4569,10 +4516,7 @@ subroutine li_flood_fill(seedMask, growMask, domain) ! do anyway. call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + !call mpas_log_write("Flood-fill: updated masks.") newMaskCountLocal = sum(seedMask) !call mpas_log_write("Initialized $i cells to local mask", intArgs=(/newMaskCountLocal/)) 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 4e16df1c1c34..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 @@ -287,7 +287,7 @@ function li_core_init(domain, startTimeStamp) result(err) call mpas_pool_get_array(meshPool, 'cellProcID', cellProcID) cellProcID(:) = domain % dminfo % my_proc_id - call landice_init_block(block, domain, domain % dminfo, err_tmp) + call landice_init_block(block, domain % dminfo, err_tmp) err = ior(err, err_tmp) block => block % next @@ -804,10 +804,6 @@ function li_core_initial_solve(domain) result(err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") call li_update_geometry(geometryPool) block => block % next @@ -928,7 +924,7 @@ end function li_core_initial_solve !> This routine initializes blocks for the land ice core. ! !----------------------------------------------------------------------- - subroutine landice_init_block(block, domain, dminfo, err) + subroutine landice_init_block(block, dminfo, err) use mpas_derived_types use mpas_pool_routines @@ -946,13 +942,14 @@ subroutine landice_init_block(block, domain, dminfo, err) ! !----------------------------------------------------------------- type (dm_info), intent(in) :: dminfo !< Input: Domain info + !----------------------------------------------------------------- ! ! input/output variables ! !----------------------------------------------------------------- type (block_type), intent(inout) :: block !< Input/output: Block object - type (domain_type), intent(inout) :: domain !< Input/output: Domain + !----------------------------------------------------------------- ! ! output variables @@ -1042,10 +1039,7 @@ subroutine landice_init_block(block, domain, dminfo, err) ! Mask init identifies initial ice extent call li_calculate_mask_init(geometryPool, err=err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! Initialize thicknessOld variable used to calculate dHdt call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'thicknessOld', thicknessOld) 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 b830d41179ce..9d74defa2201 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 @@ -175,10 +175,7 @@ subroutine li_face_melt_grounded_ice(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + block => block % next enddo ! associated(block) @@ -417,10 +414,7 @@ subroutine li_basal_melt_floating_ice(domain, err) ! calculate masks - so we know where the ice is floating call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! calculate a mask to identify ice that is thick enough to be thermally active do iCell = 1, nCellsSolve if (thickness(iCell) > config_thermal_thickness) then 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 7364966826e9..4be457cb50c9 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,10 +172,6 @@ subroutine li_SGH_init(domain, err) call calc_hydro_mask(domain) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") ! remove invalid values - not necessary on restart, but shouldn't hurt call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness) waterThickness = max(0.0_RKIND, waterThickness) 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 ddbb019551c0..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 @@ -860,10 +860,7 @@ subroutine li_thermal_solver(domain, err) ! calculate masks - so we know where the ice is floating call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! calculate a mask to identify ice that is thick enough to be thermally active do iCell = 1, nCellsSolve if (thickness(iCell) > config_thermal_thickness) then 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 2249f7f2c141..dad87f4d62b1 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 @@ -1206,10 +1206,6 @@ subroutine advection_solver(domain, err) endif call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") 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 df9a46bb2c7d..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 @@ -334,10 +334,6 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") 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 9289a480beed..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 @@ -876,10 +876,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! code is currently organized, that would be a circular dependency.a call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err) - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'groundedMaskForMassBudget') - call mpas_dmpar_field_halo_exch(domain, 'floatingMaskForMassBudget') - call mpas_timer_stop("halo updates") + ! 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) ) lowerSurface = config_sea_level - thickness * (config_ice_density / config_ocean_density) 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 7bef9d6f970e..19cf9c1ef6ee 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -532,29 +532,29 @@ 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 - ! 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 + ! Update groundedToFloatingThickness based on new masks + do iCell = 1, nCells + ! find locations that had grounded ice before but now have floating ice + if ( (groundedMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (floatingMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) + thickness(iCell) + ! find locations that had floating ice before but now have grounded ice + else if ( (floatingMaskForMassBudgetTemporaryField % array(iCell) .eq. 1) .and. & + (groundedMaskForMassBudget(iCell) .eq. 1) ) then + groundedToFloatingThickness(iCell) = groundedToFloatingThickness(iCell) - thickness(iCell) + endif + enddo + ! Now use flood fill to identify all floating ice connected to the open + ! ocean. Subglacial lakes are any floating ice that are not in the + ! resulting flood fill mask. Commented code below does not yet work + ! because flood fill needs to be moved from li_calving to a shared + ! routine. + ! call li_flood_fill(seedMask=seedMask, growMask=li_mask_is_floating_ice(cellMask), domain) + + !where (li_mask_is_floating(cellMask) and seedMask .eq. 0) + ! groundedMaskForBudget = 1 + ! gloatingMaskForMassBudget = 0 + !end where endif !call mpas_timer_stop('calculate mask cell') From 5976ab00e8b83420b9ddf225d1c224af42ac42e9 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 3 Aug 2022 18:09:22 -0600 Subject: [PATCH 13/33] 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 19cf9c1ef6ee..4e9e1e9b9c0a 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -513,8 +513,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 @@ -526,8 +526,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 02ee6a8739bc64e75909e0154519798840015e72 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 5 Oct 2022 11:40:34 -0600 Subject: [PATCH 14/33] 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 7164ab211e41..bfc963a5785c 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 @@ -1039,6 +1039,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 @@ -1053,11 +1054,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 797cf991ae61e5a8889509e9cd45fa4b1c848e5f Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 5 Oct 2022 12:56:48 -0600 Subject: [PATCH 15/33] 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 c5abcbba500c..73709f2bab4f 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 @@ -586,18 +586,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 a5bee82292707bdc6afe6da342be345fc16571e3 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 1 Nov 2022 19:30:05 -0600 Subject: [PATCH 16/33] 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 | 271 +++++++++--------- 1 file changed, 135 insertions(+), 136 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 bfc963a5785c..e46717235d31 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 @@ -515,6 +515,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 @@ -660,6 +662,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. @@ -672,7 +677,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) @@ -680,10 +685,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)) @@ -824,6 +831,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 @@ -891,7 +900,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 @@ -981,14 +991,14 @@ subroutine thickness_calving(domain, calvingFraction, err) ! Calve ice in ocean cells that are not on the protected inactive margin where (oceanMask == 1 .and. inactiveMarginMask == 0 .and. thickness > 0.0_RKIND) - calvingThickness = thickness * calvingFraction + dCalvingThickness = thickness * calvingFraction endwhere calvingThicknessFromThreshold = calvingThickness ! === apply calving === - thickness(:) = thickness(:) - calvingThickness(:) + thickness(:) = thickness(:) - dCalvingThickness(:) - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_small_islands(meshPool, geometryPool) block => block % next @@ -998,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 @@ -1037,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 @@ -1055,22 +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(:) - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_small_islands(meshPool, geometryPool) + deallocate(dCalvingThickness) block => block % next enddo @@ -1101,6 +1116,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 :: thickness real (kind=RKIND), dimension(:), pointer :: bedTopography integer, dimension(:), pointer :: cellMask, & @@ -1144,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 @@ -1299,8 +1316,8 @@ subroutine remove_small_islands(domain, err) call li_flood_fill(seedMask, growMask, domain) do iCell = 1, nCells if (li_mask_is_floating_ice(cellMask(iCell)) .and. seedMask(iCell) == 0) then - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) calvingThicknessFromThreshold(iCell) = calvingThicknessFromThreshold(iCell) + thickness(iCell) + dCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif enddo @@ -1311,6 +1328,9 @@ 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) + + deallocate(dCalvingThickness) end subroutine remove_small_islands @@ -1348,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 @@ -1376,20 +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(:) - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_small_islands(meshPool, geometryPool) + deallocate(dCalvingThickness) + block => block % next enddo @@ -1447,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 @@ -1528,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 @@ -1535,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) @@ -1547,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 & @@ -1569,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) @@ -1590,39 +1619,20 @@ 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 - call update_calving_budget(domain) + call remove_icebergs(domain) call remove_small_islands(meshPool, geometryPool) + deallocate(dCalvingThickness) block => block % next enddo @@ -1674,6 +1684,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 @@ -1727,6 +1738,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 @@ -1745,7 +1758,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, & @@ -1754,8 +1767,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) @@ -1774,40 +1787,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 - call update_calving_budget(domain) call remove_small_islands(meshPool, geometryPool) block => block % next + + deallocate(dCalvingThickness) enddo end subroutine specified_calving_velocity @@ -1868,6 +1863,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 @@ -1980,7 +1976,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 @@ -2116,7 +2114,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) @@ -2140,7 +2138,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 @@ -2151,12 +2149,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 @@ -2164,7 +2169,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 & @@ -2177,25 +2182,17 @@ 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) call remove_small_islands(meshPool, geometryPool) + call remove_icebergs(domain) + deallocate(dCalvingThickness) block => block % next enddo ! associated(block) @@ -2262,6 +2259,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 @@ -2322,6 +2320,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 @@ -2445,7 +2446,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, & @@ -2460,15 +2461,15 @@ 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) call remove_small_islands(meshPool, geometryPool) deallocate(submergedArea) - + deallocate(dCalvingThickness) end subroutine ismip6_retreat !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -3401,6 +3402,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 @@ -3464,6 +3466,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) @@ -3485,15 +3489,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) @@ -3505,38 +3509,21 @@ 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 - call update_calving_budget(domain) + call remove_icebergs(domain) call remove_small_islands(meshPool, geometryPool) + + deallocate(dCalvingThickness) + block => block % next enddo @@ -4011,7 +3998,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 !----------------------------------------------------------------- @@ -4023,12 +4010,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 !----------------------------------------------------------------- @@ -4069,6 +4055,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 @@ -4088,13 +4076,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 @@ -4149,6 +4137,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 @@ -4170,6 +4159,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) @@ -4181,7 +4174,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 @@ -4191,12 +4184,14 @@ 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) call remove_small_islands(meshPool, geometryPool) + deallocate(dCalvingThickness) block => block % next enddo @@ -4311,6 +4306,7 @@ subroutine 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, & @@ -4351,6 +4347,9 @@ subroutine remove_icebergs(domain) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + allocate(dCalvingThickness(nCells+1)) + dCalvingThickness(:) = 0.0_RKIND + seedMask(:) = 0 growMask(:) = 0 @@ -4424,19 +4423,15 @@ subroutine 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 @@ -4454,6 +4449,8 @@ subroutine 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") + + deallocate(dCalvingThickness) end subroutine remove_icebergs !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| @@ -4620,23 +4617,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) @@ -4644,10 +4640,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 b89fb71bd32c44d22970cba721dc17eb20d8a991 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 1 Nov 2022 20:21:56 -0600 Subject: [PATCH 17/33] 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_advection.F | 2 -- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 2 +- 2 files changed, 1 insertion(+), 3 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 73709f2bab4f..cb65d6e98380 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 @@ -194,7 +194,6 @@ subroutine li_advection_thickness_tracers(& basalTracersField ! scratch field containing values of basal tracers type (field1DInteger), pointer :: & - cellMaskTemporaryField, & ! scratch field containing old values of cellMask thermalCellMaskField ! Allocatable arrays need for flux-corrected transport advection @@ -731,7 +730,6 @@ subroutine li_advection_thickness_tracers(& call mpas_deallocate_scratch_field(layerThicknessOldField, .true.) call mpas_deallocate_scratch_field(basalTracersField, .true.) call mpas_deallocate_scratch_field(surfaceTracersField, .true.) - call mpas_deallocate_scratch_field(cellMaskTemporaryField, .true.) call mpas_deallocate_scratch_field(thermalCellMaskField, .true.) ! === error check 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 e46717235d31..fbed6f6b5386 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 @@ -662,7 +662,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 236074d612fdb7779ecb4a458e8067a2e4c57ac8 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 3 Nov 2022 10:21:51 -0600 Subject: [PATCH 18/33] 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 fbed6f6b5386..7f1e91134431 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 @@ -4640,9 +4640,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 e7fe7051dcf9338211e48adc29e8e1fde67de831 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 23 Feb 2023 16:10:59 -0700 Subject: [PATCH 19/33] 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 3518fc09c394..06c588585f45 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1259,6 +1259,7 @@ is the value of that variable from the *previous* time level! /> @@ -1347,7 +1348,6 @@ is the value of that variable from the *previous* time level! - /> @@ -1420,6 +1420,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 7f1e91134431..bea38c320845 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 @@ -4165,7 +4165,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 884049fd760534ccd46f9f70ad2dcb8930ae4889 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 17 Oct 2023 14:16:17 -0600 Subject: [PATCH 20/33] 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 8e85cd8e5833ec3427366616b854d7cb1e131b1b Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 17 Oct 2023 16:26:14 -0600 Subject: [PATCH 21/33] 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 bea38c320845..2c93e65140da 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,6 +663,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 eb6748ff857afa88a56bc8c14cca3c35ecf8df6b Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 Oct 2023 16:00:19 -0600 Subject: [PATCH 22/33] 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 cb65d6e98380..674af3a9d01a 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 @@ -1052,27 +1052,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 6f13e2b9e71c9341987bdcc9cbba38e856aef883 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 18 Oct 2023 17:05:25 -0600 Subject: [PATCH 23/33] 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 6a2d6af7b2f8d7b749b1d3ff6da086182f81c37e Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 19 Oct 2023 20:34:19 -0600 Subject: [PATCH 24/33] 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 674af3a9d01a..ac721a7b5bcf 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 @@ -566,6 +566,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 4b77dba818c694fcb3c5e42420273766c4db8568 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 20 Oct 2023 10:49:40 -0600 Subject: [PATCH 25/33] 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 ac721a7b5bcf..d122bff05dd3 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 @@ -371,38 +371,45 @@ subroutine li_advection_thickness_tracers(& ! update masks call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) + ! 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 @@ -476,7 +483,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 4f756137687a3dffc2f9c8a1efbb1bb81a931c87 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 20 Oct 2023 10:51:00 -0600 Subject: [PATCH 26/33] 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. --- .../mpas_li_time_integration_fe_rk.F | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 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 dad87f4d62b1..5437d4c0de6b 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 @@ -238,15 +238,6 @@ 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 -! === 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") - -! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? *** - ! === Basal melting for floating ice =========== call mpas_timer_start("basal melting for floating ice") call li_basal_melt_floating_ice(domain, err_tmp) @@ -552,6 +543,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") +!TODO: Determine whether grounded melting should in fact be called first +! === 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 ad484aa5112f6cbcaf6b62b091e67ea1f5c21084 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 16 Nov 2023 11:46:14 -0800 Subject: [PATCH 27/33] 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 | 11 ++--------- 2 files changed, 2 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 d122bff05dd3..8582e8580d05 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 5437d4c0de6b..7ed7c717da91 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 @@ -84,7 +84,6 @@ 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 !----------------------------------------------------------------- ! input variables @@ -377,7 +376,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) thickness = sum(layerThickness, 1) ! Calculate masks after updating thickness - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) if (trim(config_thermal_solver) .ne. 'none') then @@ -481,7 +480,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) @@ -493,10 +492,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(:) @@ -1103,7 +1098,6 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & - domain, & err_tmp, & advectTracersIn = .true.) @@ -1123,7 +1117,6 @@ subroutine advection_solver(domain, err) geometryPool, & thermalPool, & scratchPool, & - domain, & err_tmp, & advectTracersIn = .false.) From 797b72af2caae6aabd0680aff14a1e83452e3e6d Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 16 Nov 2023 14:23:32 -0700 Subject: [PATCH 28/33] 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 06c588585f45..d681ab8d8489 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1449,7 +1449,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 611a77a46adb18963f0e61b710237ee7c1c526f2 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 16 Nov 2023 15:10:03 -0700 Subject: [PATCH 29/33] 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 8582e8580d05..41614ce60f4d 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 @@ -238,7 +238,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 @@ -321,6 +321,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) @@ -371,6 +372,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") !----------------------------------------------------------------- @@ -383,7 +385,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) @@ -416,6 +418,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 1bdc08b8e94f83103963f7b79cf78744ca0ab6c3 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 16 Nov 2023 19:16:58 -0700 Subject: [PATCH 30/33] 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 41614ce60f4d..5571caa865d0 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 @@ -381,47 +381,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') @@ -582,6 +541,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 0e6025e47dc670bf8cb4806df55c3f87fab64381 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 17 Apr 2024 14:41:03 -0600 Subject: [PATCH 31/33] 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, 7 insertions(+), 19 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 2c93e65140da..39b67175fd97 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 @@ -331,7 +331,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 @@ -340,16 +340,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 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 @@ -409,7 +409,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) @@ -1000,7 +1000,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) @@ -1629,7 +1626,6 @@ subroutine eigencalving(domain, err) err = ior(err, err_tmp) call remove_icebergs(domain) - call remove_small_islands(meshPool, geometryPool) deallocate(dCalvingThickness) block => block % next @@ -1796,9 +1792,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 remove_small_islands(meshPool, geometryPool) - block => block % next deallocate(dCalvingThickness) @@ -2188,7 +2181,6 @@ subroutine von_Mises_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call remove_small_islands(meshPool, geometryPool) call remove_icebergs(domain) deallocate(dCalvingThickness) @@ -2465,7 +2457,6 @@ subroutine ismip6_retreat(domain, err) ! update mask call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call remove_small_islands(meshPool, geometryPool) deallocate(submergedArea) deallocate(dCalvingThickness) @@ -3517,9 +3508,8 @@ subroutine damage_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call update_calving_budget(domain) + call update_calving_budget(geometryPool, dCalvingThickness) call remove_icebergs(domain) - call remove_small_islands(meshPool, geometryPool) deallocate(dCalvingThickness) @@ -4188,8 +4178,6 @@ subroutine mask_calving(domain, err) call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - call remove_small_islands(meshPool, geometryPool) - deallocate(dCalvingThickness) block => block % next enddo From 0afc6db076169d66b6ab3068b3c90f3ec06982a3 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 10 Jul 2024 13:39:06 -0700 Subject: [PATCH 32/33] 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 39b67175fd97..6eec95ccca7c 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 @@ -351,29 +351,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 a54c8cbb934ea6a2905d5817bf420cf6bf8b02a9 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 10 Jul 2024 13:41:01 -0700 Subject: [PATCH 33/33] 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 | 216 +++++++++--------- 2 files changed, 105 insertions(+), 114 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index d681ab8d8489..91bd2953b3cc 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1254,6 +1254,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 6eec95ccca7c..7ddec51a8d98 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 @@ -492,8 +492,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 @@ -547,6 +547,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) @@ -639,8 +640,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 @@ -655,7 +655,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) @@ -663,9 +663,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 @@ -708,7 +707,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 @@ -809,8 +808,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 @@ -878,8 +877,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 @@ -969,14 +967,14 @@ subroutine thickness_calving(domain, calvingFraction, err) ! Calve ice in ocean cells that are not on the protected inactive margin where (oceanMask == 1 .and. inactiveMarginMask == 0 .and. thickness > 0.0_RKIND) - 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) block => block % next enddo @@ -985,7 +983,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 +1022,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 +1044,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 +1087,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 +1120,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 +1131,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 +1287,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 +1298,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 +1337,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 +1363,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(:) + thickness(:) = thickness(:) - incrementalCalvingThickness(:) - call update_calving_budget(geometryPool, dCalvingThickness) - - deallocate(dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) block => block % next enddo @@ -1446,7 +1439,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 +1497,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 +1522,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 +1530,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 +1542,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 +1561,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,19 +1585,18 @@ 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) call remove_icebergs(domain) - deallocate(dCalvingThickness) block => block % next enddo @@ -1656,7 +1648,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 @@ -1702,6 +1694,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) @@ -1710,8 +1703,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 @@ -1730,17 +1722,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) @@ -1759,19 +1753,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 @@ -1832,7 +1825,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 @@ -1897,6 +1890,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) @@ -1945,8 +1939,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 @@ -2083,7 +2076,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) @@ -2095,7 +2088,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 @@ -2107,7 +2100,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 @@ -2118,15 +2111,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) @@ -2138,7 +2131,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 & @@ -2152,15 +2145,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) @@ -2227,7 +2219,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 @@ -2271,6 +2263,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) @@ -2289,8 +2282,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 @@ -2414,7 +2406,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, & @@ -2425,18 +2417,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 !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -3369,7 +3360,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 @@ -3422,6 +3413,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) @@ -3433,8 +3425,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) @@ -3456,15 +3447,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) @@ -3476,20 +3470,18 @@ 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 update_calving_budget(geometryPool, dCalvingThickness) + call update_calving_budget(geometryPool, incrementalCalvingThickness) call remove_icebergs(domain) - deallocate(dCalvingThickness) - block => block % next enddo @@ -3964,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 !----------------------------------------------------------------- @@ -3980,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 !----------------------------------------------------------------- @@ -4021,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) ) @@ -4042,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 @@ -4103,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 @@ -4121,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) @@ -4140,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 @@ -4150,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 @@ -4270,7 +4261,7 @@ subroutine 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, & @@ -4311,9 +4302,6 @@ subroutine remove_icebergs(domain) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - allocate(dCalvingThickness(nCells+1)) - dCalvingThickness(:) = 0.0_RKIND - seedMask(:) = 0 growMask(:) = 0 @@ -4374,6 +4362,7 @@ subroutine remove_icebergs(domain) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + call mpas_pool_get_array(geometryPool, '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) @@ -4388,14 +4377,14 @@ subroutine 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 @@ -4414,7 +4403,6 @@ subroutine 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 !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| @@ -4581,12 +4569,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 +4592,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