diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 1f8d45e88d..c24bd82f73 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -277,24 +277,24 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu') call register_restart_field(CS%v_shelf, "v_shelf", .false., restart_CS, & "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu') - call register_restart_field(CS%u_bdry_val, "u_bdry", .false., restart_CS, & + call register_restart_field(CS%u_bdry_val, "u_bdry_val", .false., restart_CS, & "ice sheet/shelf boundary u-velocity", "m s-1", hor_grid='Bu') - call register_restart_field(CS%v_bdry_val, "v_bdry", .false., restart_CS, & + call register_restart_field(CS%v_bdry_val, "v_bdry_val", .false., restart_CS, & "ice sheet/shelf boundary v-velocity", "m s-1", hor_grid='Bu') - call register_restart_field(CS%u_face_mask_bdry, "u_bdry_mask", .false., restart_CS, & + call register_restart_field(CS%u_face_mask_bdry, "u_face_mask_bdry", .false., restart_CS, & "ice sheet/shelf boundary u-mask", "nondim", hor_grid='Bu') - call register_restart_field(CS%v_face_mask_bdry, "v_bdry_mask", .false., restart_CS, & + call register_restart_field(CS%v_face_mask_bdry, "v_face_mask_bdry", .false., restart_CS, & "ice sheet/shelf boundary v-mask", "nondim", hor_grid='Bu') call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, & "Average open ocean depth in a cell","m") call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & "fractional degree of grounding", "nondim") - call register_restart_field(CS%C_basal_friction, "tau_b_beta", .true., restart_CS, & + call register_restart_field(CS%C_basal_friction, "C_basal_friction", .true., restart_CS, & "basal sliding coefficients", "Pa (m s-1)^n_sliding") - call register_restart_field(CS%AGlen_visc, "A_Glen", .true., restart_CS, & + call register_restart_field(CS%AGlen_visc, "AGlen_visc", .true., restart_CS, & "ice-stiffness parameter", "Pa-3 s-1") - call register_restart_field(CS%h_bdry_val, "h_bdry", .false., restart_CS, & + call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, & "ice thickness at the boundary","m") endif @@ -503,6 +503,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%ground_frac,G%domain) call pass_var(CS%ice_visc,G%domain) call pass_var(CS%basal_traction, G%domain) + call pass_var(CS%AGlen_visc, G%domain) call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) endif @@ -533,30 +534,31 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif ! initialize basal friction coefficients - call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file) - call pass_var(CS%C_basal_friction, G%domain) + if (new_sim) then + call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file) + call pass_var(CS%C_basal_friction, G%domain) - ! initialize ice-stiffness AGlen - call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) - call pass_var(CS%AGlen_visc, G%domain) + ! initialize ice-stiffness AGlen + call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) + call pass_var(CS%AGlen_visc, G%domain) - !initialize boundary conditions - call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & + !initialize boundary conditions + call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & CS%u_bdry_val, CS%v_bdry_val, CS%umask, CS%vmask, CS%h_bdry_val, & CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, US, param_file ) - call pass_var(ISS%hmask, G%domain) - call pass_var(CS%h_bdry_val, G%domain) - call pass_var(CS%thickness_bdry_val, G%domain) - call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) - call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) - - !initialize ice flow velocities from file - call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, & + call pass_var(ISS%hmask, G%domain) + call pass_var(CS%h_bdry_val, G%domain) + call pass_var(CS%thickness_bdry_val, G%domain) + call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) + call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) + + !initialize ice flow velocities from file + call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, & G, US, param_file) - call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) - call pass_var(CS%bed_elev, G%domain,CENTER) - call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) - + call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) + call pass_var(CS%bed_elev, G%domain,CENTER) + call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) + endif ! Register diagnostics. CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesB1, Time, & 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) @@ -564,16 +566,15 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) ! I think that the conversion factors for the next two diagnostics are wrong. - RWH CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesB1, Time, & - 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) + 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, & - 'y-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) + 'y-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa) CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, & 'mask for u-nodes', 'none') CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, & 'mask for v-nodes', 'none') CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, & 'fraction of cell that is grounded', 'none') -! CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & @@ -582,12 +583,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'taub', 'Pa yr m-1', conversion=1e-6*US%Z_to_m) CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) - if (new_sim) then - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") - call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) - endif endif + call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") + call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) end subroutine initialize_ice_shelf_dyn @@ -683,46 +682,46 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding ! - call ice_shelf_advect(CS, ISS, G, time_step, Time) - CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step - if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. - - if (coupled_GL) then - 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(:,:)) - endif + call ice_shelf_advect(CS, ISS, G, time_step, Time) + CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step + if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. + if (coupled_GL) then + 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(:,:)) + endif - if (update_ice_vel) then - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) - endif + + if (update_ice_vel) then + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) + endif ! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) - if (update_ice_vel) then - call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) - if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) - if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) - if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) + if (update_ice_vel) then + call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) + if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) + if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) + if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) ! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) - if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) - if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) - if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) - if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) - if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) - if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) + if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) + if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) + if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag) !! - if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) - if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) - if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag) - if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag) + if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) + if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) + if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag) + if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag) ! if (CS%id_t_mask > 0) call post_data(CS%id_t_mask,CS%tmask,CS%diag) - call disable_averaging(CS%diag) + call disable_averaging(CS%diag) - CS%elapsed_velocity_time = 0.0 - endif + CS%elapsed_velocity_time = 0.0 + endif end subroutine update_ice_shelf @@ -869,13 +868,13 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite CS%ground_frac(:,:) = 0.0 allocate(Phisub(nsub,nsub,2,2,2,2)) ; Phisub(:,:,:,:,:,:) = 0.0 - do j=G%jsc,G%jec - do i=G%isc,G%iec - if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then + do j=G%jsc,G%jec + do i=G%isc,G%iec + if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then float_cond(i,j) = 1.0 CS%ground_frac(i,j) = 1.0 - endif - enddo + endif + enddo enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) @@ -1043,16 +1042,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite call MOM_mesg(mesg, 5) if (err_max <= CS%nonlinear_tolerance * err_init) then - write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init - call MOM_mesg(mesg) - write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" -! call MOM_mesg(mesg, 5) - call MOM_mesg(mesg) exit endif enddo + write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init + call MOM_mesg(mesg) + write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" + call MOM_mesg(mesg) deallocate(Phi) deallocate(Phisub) @@ -1086,6 +1084,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H !! iterations have converged to the specified tolerance integer, intent(out) :: iters !< The number of iterations used in the solver. type(time_type), intent(in) :: Time !< The current model time + character(len=160) :: mesg ! The text of an error message real, dimension(8,4,SZDI_(G),SZDJ_(G)), & intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian !! quadrature points surrounding the cell vertices [L-1 ~> m-1]. @@ -2586,7 +2585,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - & (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) -! CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging +! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) endif @@ -2938,7 +2937,7 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face do j=js,G%jed do i=is,G%ied - if (hmask(i,j) == 1) then + if ((hmask(i,j) == 1) .OR. (hmask(i,j) == 3)) then umask(I,j) = 1. vmask(I,j) = 1. @@ -2947,10 +2946,10 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%u_face_mask_bdry(I-1+k,j))) case (3) - ! vmask(I-1+k,J-1)=0. + vmask(I-1+k,J-1)=3. u_face_mask(I-1+k,j)=3. umask(I-1+k,J)=3. - !vmask(I-1+k,J)=0. + vmask(I-1+k,J)=3. vmask(I-1+k,J)=3. case (2) u_face_mask(I-1+k,j)=2. @@ -2973,9 +2972,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%v_face_mask_bdry(i,J-1+k))) case (3) vmask(I-1,J-1+k)=3. - umask(I-1,J-1+k)=0. + umask(I-1,J-1+k)=3. vmask(I,J-1+k)=3. - umask(I,J-1+k)=0. + umask(I,J-1+k)=3. v_face_mask(i,J-1+k)=3. case (2) v_face_mask(i,J-1+k)=2. diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index f3a5f210fc..73db36596e 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -101,7 +101,7 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U ! 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 ! Variable name in file + character(len=200) :: thickness_varname, area_varname, hmask_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 real :: len_sidestress, mask, udh @@ -125,45 +125,46 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U call get_param(PF, mdl, "ICE_AREA_VARNAME", area_varname, & "The name of the area variable in ICE_THICKNESS_FILE.", & default="area_shelf_h") - + hmask_varname="h_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) call MOM_read_data(filename,trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) - + call MOM_read_data(filename, trim(hmask_varname), hmask, G%Domain) isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec - do j=jsc,jec - do i=isc,iec + if (len_sidestress > 0.) then + do j=jsc,jec + do i=isc,iec ! taper ice shelf in area where there is no sidestress - ! but do not interfere with hmask - if ((G%geoLonCv(i,j) > len_sidestress).and. & - (len_sidestress > 0.)) then - udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j) - if (udh <= 25.0) then - h_shelf(i,j) = 0.0 - area_shelf_h(i,j) = 0.0 - else + if (G%geoLonCv(i,j) > len_sidestress) then + udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j) + if (udh <= 25.0) then + h_shelf(i,j) = 0.0 + area_shelf_h(i,j) = 0.0 + else h_shelf(i,j) = udh + endif endif - endif ! update thickness mask - if (area_shelf_h (i,j) >= G%areaT(i,j)) then - hmask(i,j) = 1. - elseif (area_shelf_h (i,j) == 0.0) then - hmask(i,j) = 0. - elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= G%areaT(i,j))) then - hmask(i,j) = 2. - else - call MOM_error(FATAL,mdl// " AREA IN CELL OUT OF RANGE") - endif + if (area_shelf_h (i,j) >= G%areaT(i,j)) then + hmask(i,j) = 1. + area_shelf_h(i,j)=G%areaT(i,j) + elseif (area_shelf_h (i,j) == 0.0) then + hmask(i,j) = 0. + elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= G%areaT(i,j))) then + hmask(i,j) = 2. + else + call MOM_error(FATAL,mdl// " AREA IN CELL OUT OF RANGE") + endif + enddo enddo - enddo - + endif end subroutine initialize_ice_thickness_from_file !> Initialize ice shelf thickness for a channel configuration