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 d906663476c2..922013ff5d91 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 @@ -184,11 +184,6 @@ subroutine li_advection_thickness_tracers(& surfaceTracersField, & ! scratch field containing values of surface tracers basalTracersField ! scratch field containing values of basal tracers - type (field1DInteger), pointer :: & - groundedMaskForMassBudgetTemporaryField, & ! scratch field containing old values of groundedMaskForMassBudget - floatingMaskForMassBudgetTemporaryField ! scratch field containing old values of floatingMaskForMassBudget - - integer, dimension(:), pointer :: & cellMask, & ! integer bitmask for cells edgeMask, & ! integer bitmask for edges @@ -328,12 +323,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.) - - ! given the old thickness, compute the thickness in each layer call li_calculate_layerThickness(meshPool, thickness, layerThickness) @@ -341,10 +330,6 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, domain, err_tmp) err = ior(err, err_tmp) - ! 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 !----------------------------------------------------------------- @@ -466,8 +451,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 !----------------------------------------------------------------- @@ -540,16 +523,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. @@ -603,8 +576,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(groundedMaskForMassBudgetTemporaryField, .true.) - call mpas_deallocate_scratch_field(floatingMaskForMassBudgetTemporaryField, .true.) ! === error check if (err > 0) then @@ -1934,59 +1905,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 grounded_to_floating - !*********************************************************************** end module li_advection diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index 1607d8ced6dd..227e6e6a5d43 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -283,6 +283,7 @@ subroutine prepare_advection(domain, err) real (kind=RKIND), dimension(:,:), pointer :: normalVelocity real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity + real (kind=RKIND), dimension(:), pointer :: groundedToFloatingThickness real (kind=RKIND), pointer :: calvingCFLdt integer, pointer :: processLimitingTimestep integer, dimension(:), pointer :: edgeMask @@ -374,8 +375,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 f7486a225e15..506f13ea2a81 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