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 b70fc3e4efd5..08ad9c1e6585 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 @@ -416,6 +416,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 @@ -521,6 +523,8 @@ subroutine li_restore_calving_front(domain, err) calvingThickness = 0.0_RKIND restoreThickness = 0.0_RKIND + allocate(dCalvingThickness(nCellsSolve+1)) + ! loop over locally owned cells do iCell = 1, nCells @@ -564,7 +568,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) = thickness(iCell) + dCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND endif ! li_mask_is_initial_ice @@ -573,10 +577,12 @@ subroutine li_restore_calving_front(domain, err) 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)) @@ -715,6 +721,9 @@ subroutine thickness_calving(domain, calvingFraction, err) real (kind=RKIND), dimension(:), pointer :: & calvingThickness ! thickness of ice that calves (computed in this subroutine) + real (kind=RKIND), dimension(:), allocatable :: & + dCalvingThickness + integer, pointer :: nCells integer :: iCell, iCellOnCell, iCellNeighbor real (kind=RKIND) :: flotationThickness ! thickness at which marine-based ice starts to float @@ -780,7 +789,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 @@ -870,13 +880,13 @@ 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 ! === 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 @@ -886,6 +896,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 @@ -923,10 +934,11 @@ subroutine floating_calving(domain, calvingFraction, err) type (block_type), pointer :: block type (mpas_pool_type), pointer :: geometryPool 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(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:), pointer :: thickness integer, dimension(:), pointer :: cellMask integer, dimension(:), pointer :: floatingMaskForMassBudget + integer, pointer :: nCells err = 0 @@ -937,25 +949,26 @@ subroutine floating_calving(domain, calvingFraction, err) ! get pools 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, '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 ! === 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 @@ -983,6 +996,7 @@ subroutine remove_small_islands(meshPool, geometryPool) 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, & @@ -1012,6 +1026,9 @@ subroutine remove_small_islands(meshPool, geometryPool) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) + allocate(dCalvingThickness(nCellsSolve+1)) + dCalvingThickness(:) = 0.0_RKIND + do iCell = 1, nCellsSolve if (li_mask_is_ice(cellMask(iCell))) then ! might as well do for both grounded or floating ! (1 or 2 cell floating masses are icebergs) @@ -1029,14 +1046,8 @@ subroutine remove_small_islands(meshPool, geometryPool) enddo 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) + dCalvingThickness(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 @@ -1053,30 +1064,20 @@ subroutine remove_small_islands(meshPool, geometryPool) if ((nIceNeighbors2 == 1) .and. (nOpenOceanNeighbors2 == nEdgesOnCell(iCell)-1)) then ! <- only neighbor with ice must have been iCell ! kill both cells - calvingThickness(iCell) = calvingThickness(iCell) + thickness(iCell) + dCalvingThickness(iCell) = thickness(iCell) thickness(iCell) = 0.0_RKIND - ! Update calving budgets - if (groundedMaskForMassBudget(iCell) .eq. 1) then - groundedCalvingThickness(iCell) = groundedCalvingThickness(iCell) + thickness(iCell) - elseif (floatingMaskForMassBudget(iCell) .eq. 1) then - floatingCalvingThickness(iCell) = floatingCalvingThickness(iCell) + thickness(iCell) - endif - calvingThickness(neighborWithIce) = calvingThickness(neighborWithIce) + thickness(neighborWithIce) + dCalvingThickness(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 - + call update_calving_budget(geometryPool, dCalvingThickness) + + deallocate(dCalvingThickness) end subroutine remove_small_islands @@ -1113,11 +1114,13 @@ subroutine topographic_calving(domain, calvingFraction, err) type (block_type), pointer :: block 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(:), 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 @@ -1140,19 +1143,23 @@ 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 ! === 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 @@ -1210,6 +1217,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 @@ -1291,6 +1299,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 @@ -1298,7 +1308,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) @@ -1310,7 +1320,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 & @@ -1332,8 +1342,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) @@ -1353,39 +1363,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 @@ -1437,6 +1428,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 @@ -1490,6 +1482,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 @@ -1508,7 +1502,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, & @@ -1517,8 +1511,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) @@ -1537,40 +1531,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 @@ -1629,6 +1605,7 @@ subroutine von_Mises_calving(domain, err) xvelmean, yvelmean, calvingThickness, & floatingVonMisesThresholdStress, & groundedVonMisesThresholdStress + real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness real (kind=RKIND), dimension(:,:), pointer :: flowParamA, & temperature, layerThickness real (kind=RKIND), pointer :: config_default_flowParamA @@ -1731,7 +1708,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 @@ -1776,7 +1755,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) @@ -1800,7 +1779,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 @@ -1811,12 +1790,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 @@ -1824,7 +1810,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 & @@ -1837,25 +1823,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) @@ -1922,6 +1900,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 @@ -1982,6 +1961,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 @@ -2105,7 +2087,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, & @@ -2120,15 +2102,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 !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! @@ -2959,6 +2941,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 @@ -3022,6 +3005,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) @@ -3043,15 +3028,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) @@ -3063,37 +3048,20 @@ 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 @@ -3568,7 +3536,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 !----------------------------------------------------------------- @@ -3580,12 +3548,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 !----------------------------------------------------------------- @@ -3624,6 +3591,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 @@ -3643,11 +3612,11 @@ 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) endif enddo - calvingThickness(iCell) = thickness(iCell) + dCalvingThickness(iCell) = thickness(iCell) endif ! if cell is calving margin enddo ! cell loop @@ -3696,8 +3665,10 @@ subroutine mask_calving(domain, err) type (mpas_pool_type), pointer :: velocityPool real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: calvingThickness + real (kind=RKIND), dimension(:), allocatable :: dCalvingThickness integer, dimension(:), pointer :: calvingMask integer, dimension(:), pointer :: cellMask + integer, pointer :: nCells integer :: err_tmp err = 0 @@ -3716,21 +3687,27 @@ 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 calvingThickness = 0.0_RKIND ! === apply calving === where (li_mask_is_floating_ice(cellMask) .and. (calvingMask >= 1)) - calvingThickness = thickness + dCalvingThickness = thickness thickness = 0.0_RKIND end where + 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 @@ -3844,6 +3821,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, & @@ -3884,6 +3862,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 @@ -3956,18 +3937,14 @@ 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 + 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 @@ -3985,6 +3962,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 !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| @@ -4152,23 +4131,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) @@ -4176,10 +4154,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