Skip to content

Some ice shelf updates #1009

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

Open
alex-huth wants to merge 4 commits into
base: dev/gfdl
Choose a base branch
from
Open

Some ice shelf updates #1009

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9fb747b
Optimized the ice-shelf CG scheme by reducing the number of times rep…
alex-huth Jun 9, 2025
2a50e10
Check that frazil is allocated before adding it to ice-shelf water fl…
alex-huth Dec 29, 2025
62fc738
Added melt_mask for ice shelves
alex-huth Sep 25, 2025
b6a9a73
Added melt_mask to ice-shelf restart
alex-huth Dec 4, 2025
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
22 changes: 15 additions & 7 deletions src/ice_shelf/MOM_ice_shelf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -359,7 +359,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
logical :: coupled_GL ! If true, the grounding line position is determined based on
! coupled ice-ocean dynamics.

logical :: add_frazil ! If true, allow frazil formation to modify ice-shelf water flux
real, parameter :: c2_3 = 2.0/3.0 ! Two thirds [nondim]
character(len=320) :: mesg ! The text of an error message
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
Expand Down Expand Up @@ -502,7 +502,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)

do i=is,ie
if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. &
(ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then
(ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo &
.and. ISS%melt_mask(i,j)>0.0) then

if (CS%threeeq) then
! Iteratively determine a self-consistent set of fluxes, with the ocean
Expand Down Expand Up @@ -821,6 +822,11 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
enddo ! i-loop
enddo ! j-loop

if (allocated(sfc_state%frazil)) then
add_frazil = .true.
else
add_frazil = .false.
endif

do j=js,je ; do i=is,ie
! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1]
Expand Down Expand Up @@ -871,7 +877,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
mass_flux(i,j) = ISS%water_flux(i,j) * ISS%area_shelf_h(i,j)

!Add frazil formation
if (ISS%hmask(i,j) == 1 .or. ISS%hmask(i,j) == 2) &
if (add_frazil .and. (ISS%hmask(i,j) == 1 .or. ISS%hmask(i,j) == 2)) &
ISS%water_flux(i,j) = ISS%water_flux(i,j) - sfc_state%frazil(i,j) * I_dt_LHF
fluxes%iceshelf_melt(i,j) = ISS%water_flux(i,j)
enddo ; enddo ! i- and j-loops
Expand Down Expand Up @@ -1920,8 +1926,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,

if (new_sim) then
! new simulation, initialize ice thickness as in the static case
call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, CS%Grid, CS%Grid_in, US, param_file, &
CS%rotate_index, CS%turns)
call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, ISS%melt_mask, CS%Grid, CS%Grid_in, &
US, param_file, CS%rotate_index, CS%turns)

! next make sure mass is consistent with thickness
do j=G%jsd,G%jed ; do i=G%isd,G%ied
Expand Down Expand Up @@ -1955,6 +1961,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,
"Ice shelf area in cell", "m2", conversion=US%L_to_m**2)
call register_restart_field(ISS%h_shelf, "h_shelf", .true., CS%restart_CSp, &
"ice sheet/shelf thickness", "m", conversion=US%Z_to_m)
call register_restart_field(ISS%melt_mask, "melt_mask", .false., CS%restart_CSp, &
"Mask that is >0 where ice-shelf melting is allowed", "none")

if (CS%calve_ice_shelf_bergs) then
call register_restart_field(ISS%calving, "shelf_calving", .true., CS%restart_CSp, &
Expand Down Expand Up @@ -2001,8 +2009,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,

if (new_sim .and. (.not. (CS%override_shelf_movement .and. CS%mass_from_file))) then
! This model is initialized internally or from a file.
call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, CS%Grid, CS%Grid_in, US, param_file,&
CS%rotate_index, CS%turns)
call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, ISS%melt_mask, CS%Grid, CS%Grid_in, &
US, param_file, CS%rotate_index, CS%turns)
! next make sure mass is consistent with thickness
do j=G%jsd,G%jed ; do i=G%isd,G%ied
if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 2) .or. (ISS%hmask(i,j) == 3)) then
Expand Down
49 changes: 25 additions & 24 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1738,9 +1738,11 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H
RHSu, RHSv, & ! Right hand side of the stress balance [R L3 Z T-2 ~> m kg s-2]
Au, Av, & ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2]
Du, Dv, & ! Velocity changes [L T-1 ~> m s-1]
sum_vec, sum_vec_2, sum_vec_3 !, &
!ubd, vbd ! Boundary stress contributions [R L3 Z T-2 ~> kg m s-2]
real :: beta_k, dot_p1, resid0tol2, cg_halo, max_cg_halo
sum_vec
real, dimension(SZDIB_(G),SZDJB_(G),3) :: sum_vec_3d
real :: beta_k, resid0tol2, cg_halo, max_cg_halo
real :: sv3dsum ! sum of sum_vec_3d
real :: sv3dsums(3) ! layer sums of sum_vec_3d
real :: alpha_k ! A scaling factor for iterative corrections [nondim]
real :: resid_scale ! A scaling factor for redimensionalizing the global residuals [m2 L-2 ~> 1]
! [m2 L-2 ~> 1] [R L3 Z T-2 ~> m kg s-2]
Expand All @@ -1763,7 +1765,6 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H
Zu(:,:) = 0 ; Zv(:,:) = 0 ; DIAGu(:,:) = 0 ; DIAGv(:,:) = 0
Ru(:,:) = 0 ; Rv(:,:) = 0 ; Au(:,:) = 0 ; Av(:,:) = 0 ; RHSu(:,:) = 0 ; RHSv(:,:) = 0
Du(:,:) = 0 ; Dv(:,:) = 0
dot_p1 = 0

! Determine the loop limits for sums, bearing in mind that the arrays will be starting at 1.
! Includes the edge of the tile is at the western/southern bdry (if symmetric)
Expand Down Expand Up @@ -1848,23 +1849,24 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H

call pass_vector(Au,Av,G%domain, TO_ALL, BGRID_NE)

sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
sum_vec_3d(:,:,1:2) = 0.0

do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq
if (CS%umask(I,J) == 1) then
sum_vec(I,J) = resid_scale * (Zu(I,J) * Ru(I,J))
sum_vec_2(I,J) = resid_scale * (Du(I,J) * Au(I,J))
sum_vec_3d(I,J,1) = resid_scale * (Zu(I,J) * Ru(I,J))
sum_vec_3d(I,J,2) = resid_scale * (Du(I,J) * Au(I,J))
Ru_old(I,J) = Ru(I,J) ; Zu_old(I,J) = Zu(I,J)
endif
if (CS%vmask(I,J) == 1) then
sum_vec(I,J) = sum_vec(I,J) + resid_scale * (Zv(I,J) * Rv(I,J))
sum_vec_2(I,J) = sum_vec_2(I,J) + resid_scale * (Dv(I,J) * Av(I,J))
sum_vec_3d(I,J,1) = sum_vec_3d(I,J,1) + resid_scale * (Zv(I,J) * Rv(I,J))
sum_vec_3d(I,J,2) = sum_vec_3d(I,J,2) + resid_scale * (Dv(I,J) * Av(I,J))
Rv_old(I,J) = Rv(I,J) ; Zv_old(I,J) = Zv(I,J)
endif
enddo ; enddo

alpha_k = reproducing_sum( sum_vec, Is_sum, Ie_sum, Js_sum, Je_sum ) / &
reproducing_sum( sum_vec_2, Is_sum, Ie_sum, Js_sum, Je_sum )
sv3dsum = reproducing_sum( sum_vec_3d(:,:,1:2), Is_sum, Ie_sum, Js_sum, Je_sum, sums=sv3dsums(1:2) )

alpha_k = sv3dsums(1)/sv3dsums(2)

do J=js,je-1 ; do I=is,ie-1
if (CS%umask(I,J) == 1) then
Expand All @@ -1883,23 +1885,24 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H
! R,u,v,Z valid region moves in by 1

! beta_k = (Z \dot R) / (Zold \dot Rold)
sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0 ; sum_vec_3(:,:) = 0.0
sum_vec_3d(:,:,:) = 0.0; sv3dsums(:)=0.0

do J=jscq_sv,jecq ; do i=iscq_sv,iecq
if (CS%umask(I,J) == 1) then
sum_vec(I,J) = resid_scale * (Zu(I,J) * Ru(I,J))
sum_vec_2(I,J) = resid_scale * (Zu_old(I,J) * Ru_old(I,J))
sum_vec_3(I,J) = resid2_scale * Ru(I,J)**2
sum_vec_3d(I,J,1) = resid_scale * (Zu(I,J) * Ru(I,J))
sum_vec_3d(I,J,2) = resid_scale * (Zu_old(I,J) * Ru_old(I,J))
sum_vec_3d(I,J,3) = resid2_scale * Ru(I,J)**2
endif
if (CS%vmask(I,J) == 1) then
sum_vec(I,J) = sum_vec(I,J) + resid_scale * (Zv(I,J) * Rv(I,J))
sum_vec_2(I,J) = sum_vec_2(I,J) + resid_scale * (Zv_old(I,J) * Rv_old(I,J))
sum_vec_3(I,J) = sum_vec_3(I,J) + resid2_scale * Rv(I,J)**2
sum_vec_3d(I,J,1) = sum_vec_3d(I,J,1) + resid_scale * (Zv(I,J) * Rv(I,J))
sum_vec_3d(I,J,2) = sum_vec_3d(I,J,2) + resid_scale * (Zv_old(I,J) * Rv_old(I,J))
sum_vec_3d(I,J,3) = sum_vec_3d(I,J,3) + resid2_scale * Rv(I,J)**2
endif
enddo ; enddo

beta_k = reproducing_sum(sum_vec, Is_sum, Ie_sum, Js_sum, Je_sum ) / &
reproducing_sum(sum_vec_2, Is_sum, Ie_sum, Js_sum, Je_sum )
sv3dsum = reproducing_sum( sum_vec_3d, Is_sum, Ie_sum, Js_sum, Je_sum, sums=sv3dsums )

beta_k = sv3dsums(1)/sv3dsums(2)

do J=js,je-1 ; do I=is,ie-1
if (CS%umask(I,J) == 1) Du(I,J) = Zu(I,J) + beta_k * Du(I,J)
Expand All @@ -1908,10 +1911,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H

! D valid region moves in by 1

dot_p1 = reproducing_sum( sum_vec_3, Is_sum, Ie_sum, Js_sum, Je_sum )

!if sqrt(dot_p1) <= (CS%cg_tolerance * resid0)
if (dot_p1 <= resid0tol2) then
!if sqrt(sv3dsums(3)) <= (CS%cg_tolerance * resid0)
if (sv3dsums(3) <= resid0tol2) then
iters = iter
conv_flag = 1
exit
Expand Down
26 changes: 21 additions & 5 deletions src/ice_shelf/MOM_ice_shelf_initialize.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ module MOM_ice_shelf_initialize
contains

!> Initialize ice shelf thickness
subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, PF, rotate_index, turns)
subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, melt_mask, G, G_in, US, PF, rotate_index, turns)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(ocean_grid_type), intent(in) :: G_in !< The ocean's unrotated grid structure
real, dimension(SZDI_(G),SZDJ_(G)), &
Expand All @@ -41,6 +41,8 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, P
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: hmask !< A mask indicating which tracer points are
!! partly or fully covered by an ice-shelf
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: melt_mask !< A mask indicating where to allow ice-shelf melting
type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters
logical, intent(in), optional :: rotate_index !< If true, this is a rotation test
Expand All @@ -52,6 +54,7 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, P
real, allocatable, dimension(:,:) :: tmp1_2d ! Temporary array for storing ice shelf input data
real, allocatable, dimension(:,:) :: tmp2_2d ! Temporary array for storing ice shelf input data
real, allocatable, dimension(:,:) :: tmp3_2d ! Temporary array for storing ice shelf input data
real, allocatable, dimension(:,:) :: tmp4_2d ! Temporary array for storing ice shelf input data

call get_param(PF, mdl, "ICE_PROFILE_CONFIG", config, &
"This specifies how the initial ice profile is specified. "//&
Expand All @@ -64,20 +67,22 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, P
allocate(tmp1_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed), source=0.0)
allocate(tmp2_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed), source=0.0)
allocate(tmp3_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed), source=0.0)
allocate(tmp4_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed), source=1.0)
select case ( trim(config) )
case ("CHANNEL") ; call initialize_ice_thickness_channel (tmp1_2d, tmp2_2d, tmp3_2d, G_in, US, PF)
case ("FILE") ; call initialize_ice_thickness_from_file (tmp1_2d, tmp2_2d, tmp3_2d, G_in, US, PF)
case ("FILE") ; call initialize_ice_thickness_from_file (tmp1_2d, tmp2_2d, tmp3_2d, tmp4_2d, G_in, US, PF)
case ("USER") ; call USER_init_ice_thickness (tmp1_2d, tmp2_2d, tmp3_2d, G_in, US, PF)
case default ; call MOM_error(FATAL,"MOM_initialize: Unrecognized ice profile setup "//trim(config))
end select
call rotate_array(tmp1_2d,turns, h_shelf)
call rotate_array(tmp2_2d,turns, area_shelf_h)
call rotate_array(tmp3_2d,turns, hmask)
call rotate_array(tmp4_2d,turns, melt_mask)
deallocate(tmp1_2d,tmp2_2d,tmp3_2d)
else
select case ( trim(config) )
case ("CHANNEL") ; call initialize_ice_thickness_channel (h_shelf, area_shelf_h, hmask, G, US, PF)
case ("FILE") ; call initialize_ice_thickness_from_file (h_shelf, area_shelf_h, hmask, G, US, PF)
case ("FILE") ; call initialize_ice_thickness_from_file (h_shelf, area_shelf_h, hmask, melt_mask, G, US, PF)
case ("USER") ; call USER_init_ice_thickness (h_shelf, area_shelf_h, hmask, G, US, PF)
case default ; call MOM_error(FATAL,"MOM_initialize: Unrecognized ice profile setup "//trim(config))
end select
Expand All @@ -86,7 +91,7 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, P
end subroutine initialize_ice_thickness

!> Initialize ice shelf thickness from file
subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, US, PF)
subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, melt_mask, G, US, PF)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
Expand All @@ -95,13 +100,15 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: hmask !< A mask indicating which tracer points are
!! partly or fully covered by an ice-shelf
real, dimension(SZDI_(G),SZDJ_(G)), &
intent(inout) :: melt_mask !< A mask indicating where to allow ice-shelf melting
type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters

! This subroutine reads ice thickness and area from a file and puts it into
! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask
character(len=200) :: filename,thickness_file,inputdir ! Strings for file/path
character(len=200) :: thickness_varname, area_varname, hmask_varname ! Variable name in file
character(len=200) :: thickness_varname, area_varname, hmask_varname, melt_mask_varname ! Variable name in file
character(len=40) :: mdl = "initialize_ice_thickness_from_file" ! This subroutine's name.
integer :: i, j, isc, jsc, iec, jec
logical :: hmask_set
Expand All @@ -127,6 +134,9 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
"The name of the area variable in ICE_THICKNESS_FILE.", &
default="area_shelf_h")
hmask_varname="h_mask"
call get_param(PF, mdl, "MELT_MASK_VARNAME", melt_mask_varname, &
"The name of the melt mask variable in ICE_THICKNESS_FILE.", &
default="melt_mask")
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
" initialize_topography_from_file: Unable to open "//trim(filename))
call MOM_read_data(filename, trim(thickness_varname), h_shelf, G%Domain, scale=US%m_to_Z)
Expand All @@ -139,6 +149,12 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U
"from variable "//trim(hmask_varname)//", which does not exist in "//trim(filename))
hmask_set = .false.
endif
if (field_exists(filename, trim(melt_mask_varname), MOM_domain=G%Domain)) then
call MOM_read_data(filename, trim(melt_mask_varname), melt_mask, G%Domain)
else
melt_mask(:,:)=1.0
endif

isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec

if (.not.hmask_set) then
Expand Down
2 changes: 2 additions & 0 deletions src/ice_shelf/MOM_ice_shelf_state.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module MOM_ice_shelf_state
real, pointer, dimension(:,:) :: &
mass_shelf => NULL(), & !< The mass per unit area of the ice shelf or sheet [R Z ~> kg m-2].
area_shelf_h => NULL(), & !< The area per cell covered by the ice shelf [L2 ~> m2].
melt_mask => NULL(), & !< Mask is > 0 where melting is allowed
h_shelf => NULL(), & !< the thickness of the shelf [Z ~> m], redundant with mass but may
!! make the code more readable
dhdt_shelf => NULL(), & !< the change in thickness of the shelf over time [Z T-1 ~> m s-1]
Expand Down Expand Up @@ -74,6 +75,7 @@ subroutine ice_shelf_state_init(ISS, G)

allocate(ISS%mass_shelf(isd:ied,jsd:jed), source=0.0 )
allocate(ISS%area_shelf_h(isd:ied,jsd:jed), source=0.0 )
allocate(ISS%melt_mask(isd:ied,jsd:jed), source=1.0 )
allocate(ISS%h_shelf(isd:ied,jsd:jed), source=0.0 )
allocate(ISS%dhdt_shelf(isd:ied,jsd:jed), source=0.0 )
allocate(ISS%hmask(isd:ied,jsd:jed), source=-2.0 )
Expand Down