Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fjord calving fix #132

Open
wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

Fjord calving fix #132

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e182483
Adjust calving scheme to handle narrow fjords
matthewhoffman Nov 21, 2024
9962299
Add specified_retreat_velocity option
matthewhoffman Nov 21, 2024
81d8d56
Get velocities from velocityPool not geometryPool
andrewdnolan Nov 22, 2024
66e9570
Add error check to calving edge length scaling
matthewhoffman Nov 25, 2024
90d709f
Support specified_retreat_velocity in adaptive timestepper
matthewhoffman Nov 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions components/mpas-albany-landice/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@
<nml_record name="calving" in_defaults="true">
<nml_option name="config_calving" type="character" default_value="none" units="unitless"
description="Selection of the method for calving ice (as defined below). 'von_Mises_stress' and 'eigencalving' options can be used in combination with damage threshold calving (see descrption of config_damage_calving_method for details)."
possible_values="'none', 'floating', 'topographic_threshold', 'thickness_threshold', 'eigencalving', 'specified_calving_velocity', 'von_Mises_stress', 'damagecalving', 'ismip6_retreat'"
possible_values="'none', 'floating', 'topographic_threshold', 'thickness_threshold', 'eigencalving', 'specified_calving_velocity', 'specified_retreat_velocity', 'von_Mises_stress', 'damagecalving', 'ismip6_retreat'"
/>
<nml_option name="config_apply_calving_mask" type="logical" default_value=".false." units="unitless"
description="Whether the field 'calvingMask' gets used to force calving. This is independent of choice of 'config_calving'; the application of 'calvingMask' can be applied in conjunction with a physical calving law or with no calving law enabled. The calving mask gets applied after a physical calving law, if one is selected."
Expand All @@ -189,11 +189,11 @@
possible_values="any positive real number"
/>
<nml_option name="config_calving_specified_source" type="character" default_value="const" units="none"
description="method of specified calving velocity"
description="method of specified calving velocity (for config_calving='specified_calving_velocity') or specified retreat velocity (for config_calving='specified_retreat_velocity')"
possible_values="'const' (constant calving velocity), 'data' (calving velocity given by data input)"
/>
<nml_option name="config_calving_velocity_const" type="real" default_value="0.0" units="m s^{-1}"
description="constant calving velocity specified in the namelist"
description="constant velocity used by specified_calving_velocity and specified_retreat_velocity options for config_calving. If config_calving='specified_calving_velocity', then config_calving_velocity_const is the calving velocity. If config_calving=specified_retreat_velocity, then config_calving_velocity_const is the retreat velocity."
possible_values="Any positive real value"
/>
<nml_option name="config_data_calving" type="logical" default_value=".false." units="unitless"
Expand Down Expand Up @@ -1321,7 +1321,7 @@ is the value of that variable from the *previous* time level!
description="bookkeeping field for how much volume should be removed from a dynamic edge"
/>
<var name="calvingVelocityData" type="real" dimensions="nCells Time" units="m s^{-1}" time_levs="1"
description="rate of calving front retreat due to calving, represented as a velocity normal to the calving front (in the x-y plane), given by the input netCDF file."
description="If config_calving='specified_calving_velocity', then calvingVelocityData is the calving velocity to be imposed. If config_calving=specified_retreat_velocity, then calvingVelocityData is the retreat velocity to be imposed. In either case, this field is only used if config_calving_specified_source='data'"
/>
<var name="passiveTracer2d" type="real" dimensions="nCells Time" units="none" time_levs="1"
description="passive tracer used to verify advection routines"
Expand Down
109 changes: 80 additions & 29 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -283,7 +283,8 @@ subroutine li_calve_ice(domain, err, solveVeloAfterCalving)
call von_Mises_calving(domain, err_tmp)
err = ior(err, err_tmp)

elseif (trim(config_calving) == 'specified_calving_velocity') then
elseif ((trim(config_calving) == 'specified_calving_velocity') .or. &
(trim(config_calving) == 'specified_retreat_velocity')) then

call specified_calving_velocity(domain, err_tmp)
err = ior(err, err_tmp)
Expand Down Expand Up @@ -1378,7 +1379,6 @@ subroutine eigencalving(domain, err)
real (kind=RKIND), dimension(:), pointer :: eigencalvingParameter
real (kind=RKIND), dimension(:), pointer :: calvingVelocity
real (kind=RKIND), dimension(:), pointer :: eMax, eMin
real (kind=RKIND), dimension(:), pointer :: angleEdge
real (kind=RKIND), dimension(:), pointer :: thickness
real (kind=RKIND), dimension(:), pointer :: calvingThickness
real (kind=RKIND), pointer :: calvingCFLdt
Expand Down Expand Up @@ -1429,7 +1429,6 @@ subroutine eigencalving(domain, err)

! get fields
call mpas_pool_get_array(meshPool, 'deltat', deltat)
call mpas_pool_get_array(meshPool, 'angleEdge', angleEdge)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'calvingFrontMask', calvingFrontMask)
call mpas_pool_get_array(geometryPool, 'eigencalvingParameter', eigencalvingParameter)
Expand Down Expand Up @@ -1564,12 +1563,16 @@ end subroutine eigencalving
!
! routine specified_calving_velocity
!
!> \brief use a specified calving velocity given by input data
!> \brief use a specified calving or retreat velocity given by input data
!> \author Matthew Hoffman
!> \date Feb. 2018
!> \details we can specify the calving velocity by i) a constant value
!> given by config_calving_velocity_const in namelist and ii) source data
!> specified by calvingVelocityData
!> This subroutine handles both the specified_calving_velocity and
!> specified_retreat_velocity options. In the latter case, the flow speed
!> is added to the specified retreat rate to obtain a calving rate that will
!> yield the desired retreat rate.
!-----------------------------------------------------------------------
subroutine specified_calving_velocity(domain, err)

Expand Down Expand Up @@ -1597,14 +1600,15 @@ subroutine specified_calving_velocity(domain, err)
type (mpas_pool_type), pointer :: velocityPool
type (mpas_pool_type), pointer :: scratchPool
logical, pointer :: config_print_calving_info
character (len=StrKIND), pointer :: config_calving
real(kind=RKIND), pointer :: config_calving_thickness
real(kind=RKIND), pointer :: config_calving_velocity_const
character (len=StrKIND), pointer :: config_calving_specified_source
real (kind=RKIND), dimension(:), pointer :: calvingVelocity
real (kind=RKIND), dimension(:), pointer :: calvingVelocityData
real (kind=RKIND), dimension(:), pointer :: angleEdge
real (kind=RKIND), dimension(:), pointer :: thickness
real (kind=RKIND), dimension(:), pointer :: calvingThickness
real (kind=RKIND), dimension(:), pointer :: xvelmean, yvelmean
real (kind=RKIND), pointer :: calvingCFLdt
real (kind=RKIND), pointer :: dtCalvingCFLratio
integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell
Expand All @@ -1623,6 +1627,7 @@ subroutine specified_calving_velocity(domain, err)

err = 0

call mpas_pool_get_config(liConfigs, 'config_calving', config_calving)
call mpas_pool_get_config(liConfigs, 'config_print_calving_info', config_print_calving_info)
call mpas_pool_get_config(liConfigs, 'config_calving_thickness', config_calving_thickness)
call mpas_pool_get_config(liConfigs, 'config_calving_velocity_const', config_calving_velocity_const)
Expand All @@ -1643,7 +1648,6 @@ subroutine specified_calving_velocity(domain, err)

! get fields
call mpas_pool_get_array(meshPool, 'deltat', deltat)
call mpas_pool_get_array(meshPool, 'angleEdge', angleEdge)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'calvingFrontMask', calvingFrontMask)
call mpas_pool_get_array(geometryPool, 'calvingVelocity', calvingVelocity)
Expand Down Expand Up @@ -1675,6 +1679,12 @@ subroutine specified_calving_velocity(domain, err)
realArgs=(/minval(calvingVelocity), maxval(calvingVelocity)/))
endif

if (trim(config_calving) == 'specified_retreat_velocity') then
call mpas_pool_get_array(velocityPool, 'xvelmean', xvelmean)
call mpas_pool_get_array(velocityPool, 'yvelmean', yvelmean)
calvingVelocity = calvingVelocity + sqrt(xvelmean**2.0_RKIND + yvelmean**2.0_RKIND)
endif

! Convert calvingVelocity to calvingThickness
call li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool, calvingThickness, calvingVelocity, &
applyToGrounded=.true., applyToFloating=.true., applyToGroundingLine=.false., &
Expand Down Expand Up @@ -2495,7 +2505,7 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
real (kind=RKIND), pointer :: config_calving_error_threshold
logical, pointer :: config_distribute_unablatedVolumeDynCell
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge
real (kind=RKIND), dimension(:), pointer :: angleEdge
real (kind=RKIND), dimension(:), pointer :: xCell, yCell, xEdge, yEdge
real (kind=RKIND), dimension(:), pointer :: calvingVelocity, faceMeltSpeed
real (kind=RKIND), dimension(:), pointer :: areaCell
real (kind=RKIND), dimension(:,:), pointer :: uReconstructX
Expand Down Expand Up @@ -2536,7 +2546,10 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
call mpas_pool_get_array(meshPool, 'angleEdge', angleEdge)
call mpas_pool_get_array(meshPool, 'xCell', xCell)
call mpas_pool_get_array(meshPool, 'yCell', yCell)
call mpas_pool_get_array(meshPool, 'xEdge', xEdge)
call mpas_pool_get_array(meshPool, 'yEdge', yEdge)
call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
Expand Down Expand Up @@ -2850,10 +2863,11 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
do iNeighbor = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(iNeighbor, iCell)
jCell = cellsOnCell(iNeighbor, iCell)
if (li_mask_is_margin(edgeMask(iEdge)) .and. & !< edge is a margin
.not. li_mask_is_ice(cellMask(jCell)) .and. bedTopography(jCell) < config_sea_level & ! ensure margin is w/ocn
) then
edgeLengthScaling = scale_edge_length(angleEdge(iEdge), uvelForAblation(iCell), vvelForAblation(iCell))
if (li_mask_is_margin(edgeMask(iEdge)) .or. &
(li_mask_is_ice(cellMask(jCell)) .and. bedTopography(jCell) > config_sea_level)) then
call scale_edge_length(xCell(iCell), yCell(iCell), xEdge(iEdge), yEdge(iEdge), &
uvelForAblation(iCell), vvelForAblation(iCell), edgeLengthScaling, err_tmp)
err = ior(err, err_tmp)
requiredAblationVolumeNonDynEdge(iEdge) = ablationVelocity(iCell) * &
edgeLengthScaling * dvEdge(iEdge) * thicknessForAblation(iCell) * deltat
requiredAblationVolumeNonDynCell(iCell) = requiredAblationVolumeNonDynCell(iCell) + &
Expand Down Expand Up @@ -2913,9 +2927,12 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
.and. li_mask_is_margin(cellMask(iCell)) ) then ! that is at the edge of the ice
do iNeighbor = 1, nEdgesOnCell(iCell)
jCell = cellsOnCell(iNeighbor, iCell)
if ((.not. li_mask_is_ice(cellMask(jCell))) .and. (bedTopography(jCell) < config_sea_level)) then
if ((.not. li_mask_is_ice(cellMask(jCell))) .and. (bedTopography(jCell) < config_sea_level) .or. &
(bedTopography(jCell) > config_sea_level) ) then
iEdge = edgesOnCell(iNeighbor, iCell)
edgeLengthScaling = scale_edge_length(angleEdge(iEdge), uvelForAblation(iCell), vvelForAblation(iCell))
call scale_edge_length(xCell(iCell), yCell(iCell), xEdge(iEdge), yEdge(iEdge), &
uvelForAblation(iCell), vvelForAblation(iCell), edgeLengthScaling, err_tmp)
err = ior(err, err_tmp)
requiredAblationVolumeDynEdge(iEdge) = ablationVelocity(iCell) * &
edgeLengthScaling * dvEdge(iEdge) * thicknessForAblation(iCell) * deltat
endif
Expand All @@ -2939,7 +2956,9 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
do iNeighbor = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(iNeighbor, iCell)
if (li_mask_is_dynamic_ice(edgeMask(iEdge))) then
edgeLengthScaling = scale_edge_length(angleEdge(iEdge), uvelForAblation(iCell), vvelForAblation(iCell))
call scale_edge_length(xCell(iCell), yCell(iCell), xEdge(iEdge), yEdge(iEdge), &
uvelForAblation(iCell), vvelForAblation(iCell), edgeLengthScaling, err_tmp)
err = ior(err, err_tmp)
ablateLengthEdge = edgeLengthScaling * dvEdge(iEdge)
ablateLengthCell = ablateLengthCell + ablateLengthEdge
endif
Expand All @@ -2951,7 +2970,9 @@ subroutine li_apply_front_ablation_velocity(meshPool, geometryPool, velocityPool
iEdge = edgesOnCell(iNeighbor, iCell)
jCell = cellsOnCell(iNeighbor, iCell)
if (li_mask_is_dynamic_ice(edgeMask(iEdge))) then
edgeLengthScaling = scale_edge_length(angleEdge(iEdge), uvelForAblation(iCell), vvelForAblation(iCell))
call scale_edge_length(xCell(iCell), yCell(iCell), xEdge(iEdge), yEdge(iEdge), &
uvelForAblation(iCell), vvelForAblation(iCell), edgeLengthScaling, err_tmp)
err = ior(err, err_tmp)
ablateLengthEdge = edgeLengthScaling * dvEdge(iEdge)
if (requiredAblationVolumeDynEdge(iEdge) > 0.0_RKIND) then
call mpas_log_write("Unexpectedly found a dynamic edge that already has calving assigned to it." // &
Expand Down Expand Up @@ -3178,19 +3199,49 @@ end subroutine li_apply_front_ablation_velocity

! Helper function for subroutine li_apply_front_ablation_velocity
! Calculates the amount to scale an edge length based on the orientation of the edge with the surface velocity
function scale_edge_length(angleEdgeHere, u, v)
real(kind=RKIND), intent(in) :: angleEdgeHere
real(kind=RKIND), intent(in) :: u
real(kind=RKIND), intent(in) :: v
real(kind=RKIND) :: scale_edge_length

real(kind=RKIND) :: mag

mag = sqrt(u**2 + v**2)
if (mag == 0.0_RKIND) mag = 1.0_RKIND
scale_edge_length = abs(u/mag * cos(angleEdgeHere) + v/mag * sin(angleEdgeHere)) ! dot product of unit vectors
!scale_edge_length = abs(cos(angleEdgeHere - atan2(v, u)))
end function scale_edge_length
subroutine scale_edge_length(xc, yc, xe, ye, vx, vy, scale_factor, err)
real(kind=RKIND), intent(in) :: xc, yc ! x,y coords of cell center
real(kind=RKIND), intent(in) :: xe, ye ! x,y coords of edge neighboring cell
real(kind=RKIND), intent(in) :: vx, vy ! x,y components of velocity vector
real(kind=RKIND), intent(out) :: scale_factor
integer, intent(out) :: err

real(kind=RKIND) :: ux, uy, umag, vmag

err = 0
scale_factor = 0.0_RKIND

! Calculate vector pointing from cell center to edge
ux = xe - xc
uy = ye - yc
umag = sqrt(ux**2 + uy**2)

vmag = sqrt(vx**2 + vy**2)

! check for 0-magnitude vectors
if (umag == 0.0_RKIND) then
call mpas_log_write("In scale_edge_length, cell to edge vector has magnitude zero", MPAS_LOG_ERR)
err = 1
return
endif
if (vmag == 0.0_RKIND) then
call mpas_log_write("In scale_edge_length, velocity vector has magnitude zero", MPAS_LOG_ERR)
err = 1
return
endif

! dot product of unit vectors
scale_factor = (ux / umag) * (vx / vmag) + (uy / umag) * (vy / vmag)
! set scale to 0 where vectors are not pointing the same way
! (i.e. where vector to edge points away from flow direction)
scale_factor = max(0.0, scale_factor)

if (scale_factor > 1.0_RKIND) then
call mpas_log_write("In scale_edge_length, scale_factor is greater than 1.", MPAS_LOG_ERR)
err = 1
endif

end subroutine scale_edge_length


!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Expand Down
1 change: 1 addition & 0 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1354,6 +1354,7 @@ subroutine set_timestep(allowableAdvecDt, allowableDiffDt, calvingCFLdt, faceMel
! Only include calving CFL if requested and we are using a calving option that calculates it
if (config_adaptive_timestep_include_calving .and. ( &
trim(config_calving) == 'specified_calving_velocity' .or. &
trim(config_calving) == 'specified_retreat_velocity' .or. &
trim(config_calving) == 'eigencalving' .or. &
trim(config_calving) == 'von_Mises_stress' .or. &
trim(config_calving) == 'ismip6_retreat' .or. &
Expand Down