diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 772a21483a..689feeb01b 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -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 @@ -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 @@ -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] @@ -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 @@ -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 @@ -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, & @@ -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 diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 4f8a9f0cf2..a8b53ad306 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -1738,15 +1738,24 @@ 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 ! Sum of squares of residuals in stress calculations [m2 kg2 s-4] + real, dimension(SZDIB_(G),SZDJB_(G),3) :: sum_vec_3d ! Array used for various residuals + ! sum_vec_3d(:,:,1:2) [m s-1] [m kg s-2] + ! sum_vec_3d(:,:,3) [m2 kg2 s-4] + real :: beta_k ! Ratio of residuals used to update search direction [nondim] + real :: resid0tol2 ! Convergence tolerance times the initial residual [m2 kg2 s-4] + real :: sv3dsum ! An unused variable returned when taking global sum of residuals [various] + real :: sv3dsums(3) ! The index-wise global sums of sum_vec_3d + ! sv3dsums(:,:,1:2) [m s-1] [m kg s-2] + ! sv3dsums(:,:,3) [m2 kg2 s-4] 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] + real :: resid_scale ! A scaling factor for redimensionalizing the global residuals + ! [L T-1 ~> m s-1] [R L3 Z T-2 ~> m kg s-2] real :: resid2_scale ! A scaling factor for redimensionalizing the global squared residuals - ! [m2 L-2 ~> 1] [R L3 Z T-2 ~> m kg s-2] + ! [R2 L6 Z2 T-4 ~> m2 kg2 s-4] real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim] + integer :: cg_halo ! Number of halo vertices to include during a CG iteration + integer :: max_cg_halo ! Maximum possible number of halo vertices to include in the CG iterations integer :: iter, i, j, isd, ied, jsd, jed, isc, iec, jsc, jec, is, js, ie, je integer :: Is_sum, Js_sum, Ie_sum, Je_sum ! Loop bounds for global sums or arrays starting at 1. integer :: Isdq, Iedq, Jsdq, Jedq, Iscq, Iecq, Jscq, Jecq, nx_halo, ny_halo @@ -1763,7 +1772,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) @@ -1848,23 +1856,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 @@ -1883,23 +1892,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) @@ -1908,10 +1918,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 diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index a120f8a45e..d3d4ceb0a3 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -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)), & @@ -40,7 +40,9 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, P intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: hmask !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf + !! partly or fully covered by an ice-shelf [nondim] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: melt_mask !< A mask indicating where to allow ice-shelf melting [nondim] 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 @@ -49,9 +51,10 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, P character(len=40) :: mdl = "initialize_ice_thickness" ! This subroutine's name. character(len=200) :: config logical :: rotate = .false. - 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(:,:) :: tmp1_2d ! Temporary array for storing ice shelf input data [Z~>m] + real, allocatable, dimension(:,:) :: tmp2_2d ! Temporary array for storing ice shelf input data [L2~>m2] + real, allocatable, dimension(:,:) :: tmp3_2d ! Temporary array for storing ice shelf input data [nondim] + real, allocatable, dimension(:,:) :: tmp4_2d ! Temporary array for storing ice shelf input data [nondim] call get_param(PF, mdl, "ICE_PROFILE_CONFIG", config, & "This specifies how the initial ice profile is specified. "//& @@ -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 @@ -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]. @@ -94,14 +99,16 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [L2 ~> m2]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: hmask !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf + !! partly or fully covered by an ice-shelf [nondim] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: melt_mask !< A mask indicating where to allow ice-shelf melting [nondim] 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 @@ -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) @@ -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 diff --git a/src/ice_shelf/MOM_ice_shelf_state.F90 b/src/ice_shelf/MOM_ice_shelf_state.F90 index d789c08bd4..6a4dee9a0e 100644 --- a/src/ice_shelf/MOM_ice_shelf_state.F90 +++ b/src/ice_shelf/MOM_ice_shelf_state.F90 @@ -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 [nondim] 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] @@ -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 )