From d1a86ce81a91771ee46b4376e20aa1c4cacb2327 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Fri, 27 Oct 2023 16:00:40 -0400 Subject: [PATCH 1/2] A few bug fixes so that the GL_couple=.true. option works correctly. Setting GL_couple=.true. will determine the grounding based on ocean column thickness rather than the typical the hydrostatic equilibrium condition. This has the advantage of accounting for changes in sea level, tides, etc. However, it has the disadvantage of not working with the same thoroughly-tested sub-element grounding line parameterization used for the hydrostatic condition. Instead, it accounts for sub-element grounding line movement by, during the SSA solution, using a grounding mask averaged over all ocean (sub)steps that completed since the last SSA solve. Unlike the hydrostatic sub-element parameterization, the dependence of the GL_couple=.true. scheme on grid resolution has not yet been determined. Qualitatively similar grounding line retreat/advance behavior is achieved with both approaches for MISOMIP IceOcean1 on a 2km grid, but GL_couple=.true. results in a rougher grounding line position with less retreat. Note that this commit also fixed a bug in applying the hydrostatic grounding line approach without its sub-element parameterization (though the sub-element parameterization should also be used anyway). --- src/ice_shelf/MOM_ice_shelf.F90 | 1 + src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 85 ++++++++++++++---------- 2 files changed, 52 insertions(+), 34 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 84858f17bc..062b7b23a6 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1395,6 +1395,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", & default=.false., do_not_log=CS%GL_regularize) if (CS%GL_regularize) CS%GL_couple = .false. + if (CS%solo_ice_sheet) CS%GL_couple = .false. endif call get_param(param_file, mdl, "SHELF_THERMO", CS%isthermo, & diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 2965f6eac4..6ae3247190 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -402,6 +402,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", & default=.false., do_not_log=CS%GL_regularize) if (CS%GL_regularize) CS%GL_couple = .false. + if (present(solo_ice_sheet_in)) then + if (solo_ice_sheet_in) CS%GL_couple = .false. + endif if (CS%GL_regularize .and. (CS%n_sub_regularize == 0)) call MOM_error (FATAL, & "GROUNDING_LINE_INTERP_SUBGRID_N must be a positive integer if GL regularization is used") call get_param(param_file, mdl, "ICE_SHELF_CFL_FACTOR", CS%CFL_factor, & @@ -759,6 +762,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel) elseif (update_ice_vel) then call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) + CS%GL_couple=.false. endif @@ -922,8 +926,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i real, dimension(SZDIB_(G),SZDJB_(G)) :: Au, Av ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: u_last, v_last ! Previous velocities [L T-1 ~> m s-1] real, dimension(SZDIB_(G),SZDJB_(G)) :: H_node ! Ice shelf thickness at corners [Z ~> m]. - real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! An array indicating where the ice - ! shelf is floating: 0 if floating, 1 if not. + real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! If GL_regularize=true, an array indicating where the ice + ! shelf is floating: 0 if floating, 1 if not real, dimension(SZDIB_(G),SZDJB_(G)) :: Normvec ! Used for convergence character(len=160) :: mesg ! The text of an error message integer :: conv_flag, i, j, k,l, iter @@ -949,18 +953,18 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i ! need to make these conditional on GL interpolation float_cond(:,:) = 0.0 ; H_node(:,:) = 0.0 - CS%ground_frac(:,:) = 0.0 + !CS%ground_frac(:,:) = 0.0 allocate(Phisub(nsub,nsub,2,2,2,2), source=0.0) - do j=G%jsc,G%jec - do i=G%isc,G%iec + if (.not. CS%GL_couple) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then - float_cond(i,j) = 1.0 + if (CS%GL_regularize) float_cond(i,j) = 1.0 CS%ground_frac(i,j) = 1.0 CS%OD_av(i,j) =0.0 endif - enddo - enddo + enddo ; enddo + endif call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) @@ -1010,10 +1014,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") call pass_var(CS%Ee,G%domain) ! This makes sure basal stress is only applied when it is supposed to be - do j=G%jsd,G%jed ; do i=G%isd,G%ied -! CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) - CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) - enddo ; enddo + if (CS%GL_regularize) then + do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) + enddo ; enddo + else + do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) + enddo ; enddo + endif if (CS%nonlin_solve_err_mode == 1) then ! call apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, CS%ice_visc, & @@ -1085,11 +1094,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") call pass_var(CS%Ee,G%domain) ! makes sure basal stress is only applied when it is supposed to be - - do j=G%jsd,G%jed ; do i=G%isd,G%ied -! CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) - CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) - enddo ; enddo + if (CS%GL_regularize) then + do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) + enddo ; enddo + else + do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) + enddo ; enddo + endif if (CS%nonlin_solve_err_mode == 1) then !u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0 @@ -1196,8 +1209,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H intent(in) :: H_node !< The ice shelf thickness at nodal (corner) !! points [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< An array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not. + intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice + !! shelf is floating: 0 if floating, 1 if not real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf @@ -1940,18 +1953,22 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) rhoi_rhow = rho/rhow ! prelim - go through and calculate S - S(:,:) = -CS%bed_elev(:,:) + ISS%h_shelf(:,:) - ! check whether the ice is floating or grounded + if (CS%GL_couple) then + S(:,:) = -CS%bed_elev(:,:) + OD(:,:) + ISS%h_shelf(:,:) + else + S(:,:) = -CS%bed_elev(:,:) + ISS%h_shelf(:,:) + ! check whether the ice is floating or grounded - do j=jsc-G%domain%njhalo,jec+G%domain%njhalo - do i=isc-G%domain%nihalo,iec+G%domain%nihalo - if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then - S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) - else - S(i,j) = ISS%h_shelf(i,j)-CS%bed_elev(i,j) - endif + do j=jsc-G%domain%njhalo,jec+G%domain%njhalo + do i=isc-G%domain%nihalo,iec+G%domain%nihalo + if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then + S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) + else + S(i,j) = ISS%h_shelf(i,j)-CS%bed_elev(i,j) + endif + enddo enddo - enddo + endif call pass_var(S, G%domain) @@ -2214,8 +2231,8 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form !! and units depend on the basal law exponent. real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< An array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not. + intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice + !! shelf is floating: 0 if floating, 1 if not real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points !! relative to sea-level [Z ~> m]. @@ -2385,8 +2402,8 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< An array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not. + intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice + !! shelf is floating: 0 if floating, 1 if not real, dimension(SZDIB_(G),SZDJB_(G)), & intent(in) :: H_node !< The ice shelf thickness at nodal !! (corner) points [Z ~> m]. @@ -2916,7 +2933,7 @@ subroutine update_OD_ffrac(CS, G, US, ocean_mass, find_avg) CS%ground_frac(i,j) = 1.0 - (CS%ground_frac_rt(i,j) * I_counter) CS%OD_av(i,j) = CS%OD_rt(i,j) * I_counter - CS%OD_rt(i,j) = 0.0 ; CS%ground_frac_rt(i,j) = 0.0 + CS%OD_rt(i,j) = 0.0 ; CS%ground_frac_rt(i,j) = 0.0; CS%OD_rt_counter = 0 enddo ; enddo call pass_var(CS%ground_frac, G%domain, complete=.false.) From 21c045716f7c8407ae1ae5b997c528f6463e1fa3 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Tue, 31 Oct 2023 15:31:02 -0400 Subject: [PATCH 2/2] for calculation of ice sheet surface height, eliminated some array syntax and added parentheses to specify order of addition --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 6ae3247190..d44748a728 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -1954,11 +1954,13 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! prelim - go through and calculate S if (CS%GL_couple) then - S(:,:) = -CS%bed_elev(:,:) + OD(:,:) + ISS%h_shelf(:,:) + do j=jsc-G%domain%njhalo,jec+G%domain%njhalo + do i=isc-G%domain%nihalo,iec+G%domain%nihalo + S(i,j) = -CS%bed_elev(i,j) + (OD(i,j) + ISS%h_shelf(i,j)) + enddo + enddo else - S(:,:) = -CS%bed_elev(:,:) + ISS%h_shelf(:,:) ! check whether the ice is floating or grounded - do j=jsc-G%domain%njhalo,jec+G%domain%njhalo do i=isc-G%domain%nihalo,iec+G%domain%nihalo if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then