diff --git a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 index 42444878ee..f125490c1d 100644 --- a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 @@ -32,7 +32,7 @@ program Shelf_main use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_domains, only : MOM_domains_init, clone_MOM_domain, pass_var use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid - use MOM_error_handler, only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe + use MOM_error_handler, only : MOM_set_verbosity, MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_file_parser, only : close_param_file @@ -149,7 +149,7 @@ program Shelf_main character(len=9) :: month character(len=16) :: calendar = 'noleap' integer :: calendar_type=-1 - + integer :: verbosity integer :: unit, io_status, ierr logical :: symmetric @@ -196,6 +196,9 @@ program Shelf_main ! Also calls the subroutine that opens run-time parameter files. call Get_MOM_Input(param_file, dirs) + call get_param(param_file, mod_name, "VERBOSITY", verbosity, default=5) + call MOM_set_verbosity(verbosity) + ! Read ocean_solo restart, which can override settings from the namelist. if (file_exists(trim(dirs%restart_input_dir)//'ice_solo.res')) then call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ice_solo.res', action=READONLY_FILE) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index eb6dd7f4d9..716d213403 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -14,6 +14,7 @@ module MOM_ice_shelf_dynamics use MOM_IS_diag_mediator, only : diag_ctrl, time_type, enable_averages, disable_averaging use MOM_domains, only : MOM_domains_init, clone_MOM_domain use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, AGRID, CORNER, CENTER +use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : MOM_grid_init, ocean_grid_type @@ -108,13 +109,6 @@ module MOM_ice_shelf_dynamics !! quadrature point for Newton iterations [T-1 ~> s-1] real, pointer, dimension(:,:,:) :: newton_str_sh => NULL() !< Engineering shear strain-rate uy+vx at each !! viscosity quadrature point for Newton iterations [T-1 ~> s-1] - real, pointer, dimension(:,:) :: newton_umid => NULL() !< Cell-averaged zonal velocity u at the current outer - !! iterate, for Newton basal drag correction [L T-1 ~> m s-1] - real, pointer, dimension(:,:) :: newton_vmid => NULL() !< Cell-averaged meridional velocity v at the current - !! outer iterate, for Newton basal drag correction [L T-1 ~> m s-1] - real, pointer, dimension(:,:) :: newton_drag_coef => NULL() !< Newton basal drag correction coefficient: - !! 2 * d(basal_trac)/d(|u|^2) * area = d(tau_b_i)/d(u_j) - basal_trac*delta_ij - !! expressed as the u_i*u_j tensor coefficient [R Z T ~> kg m-2 s] real, pointer, dimension(:,:) :: AGlen_visc => NULL() !< Ice-stiffness parameter in Glen's law ice viscosity, !! often in [Pa-3 s-1] if n_Glen is 3. real, pointer, dimension(:,:) :: u_bdry_val => NULL() !< The zonal ice velocity at inflowing boundaries @@ -129,12 +123,15 @@ module MOM_ice_shelf_dynamics !! the same as G%bathyT+Z_ref, when below sea-level. !! Sign convention: positive below sea-level, negative above. - real, pointer, dimension(:,:) :: basal_traction => NULL() !< The area-integrated taub_beta field - !! (m2 Pa s m-1, or kg s-1) related to the nonlinear part - !! of "linearized" basal stress (Pa) [R Z L2 T-1 ~> kg s-1] - !! The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011 real, pointer, dimension(:,:) :: C_basal_friction => NULL()!< Coefficient in sliding law tau_b = C u^(n_basal_fric), !! units of [R L Z T-2 (s m-1)^(n_basal_fric) ~> Pa (s m-1)^(n_basal_fric)] + real, pointer, dimension(:,:) :: coef_prefactor => NULL() !< Pre-computed area*C_basal_friction*L_T_to_m_s for + !! basal friction quadrature evaluation [R L2 Z T-1 ~> kg s-1]. + real, pointer, dimension(:,:) :: fB_elem => NULL() !< Pre-computed element-level Coulomb fB parameter + !! [(T L-1)^CF_PostPeak]; 0 for Weertman. + !! Updated each outer iteration by calc_shelf_basal_prefactors. + real :: alpha_coulomb = 1.0 !< Coulomb prefactor (CF_PostPeak-1)^(CF_PostPeak-1)/CF_PostPeak^CF_PostPeak [nondim] + real :: coulomb_pp_n !< CF_PostPeak/n_basal_fric [nondim] real, pointer, dimension(:,:) :: OD_rt => NULL() !< A running total for calculating OD_av [Z ~> m]. real, pointer, dimension(:,:) :: ground_frac_rt => NULL() !< A running total for calculating ground_frac. real, pointer, dimension(:,:) :: OD_av => NULL() !< The time average open ocean depth [Z ~> m]. @@ -214,24 +211,42 @@ module MOM_ice_shelf_dynamics !! circulation or thermodynamics. It is used to estimate the !! gravitational driving force at the shelf front (until we think of !! a better way to do it, but any difference will be negligible). + real :: rhoi_rhow !< The density of ice divided by a typical water density [nondim] + real :: rhow_rhoi !< A typical water density divided by the density of ice [nondim] real :: thresh_float_col_depth !< The water column depth over which the shelf if considered to be floating logical :: moving_shelf_front !< Specify whether to advance shelf front (and calve). logical :: calve_to_mask !< If true, calve off the ice shelf when it passes the edge of a mask. real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m]. real :: T_shelf_missing !< An ice shelf temperature to use where there is no ice shelf [C ~> degC] - real :: cg_tolerance !< The tolerance in the CG solver, relative to initial residual, that + real :: cg_tolerance !< For Picard iterations, the tolerance in the CG solver, relative to initial residual, that !! determines when to stop the conjugate gradient iterations [nondim]. - real :: cg_tol_newton !< Working CG tolerance for the current inner solve [nondim]. + real :: cg_newton_tolerance !< For inexact Newton iterations, the initial tolerance in the CG solver, relative to + !! initial residual, that determines when to stop the CG iterations [nondim]. + real :: cg_tol_current !< Working CG tolerance for the current inner solve [nondim]. real :: nonlinear_tolerance !< The fractional nonlinear tolerance, relative to the initial error, !! that sets when to stop the iterative velocity solver [nondim] real :: newton_after_tolerance !< The fractional nonlinear tolerance, relative to the initial error, at - !! which to switch from Picard to Newton iterations in the velocity solver [nondim] + !! which to switch from Picard to Newton iterations in the velocity solver + !! If set to <= 0, no Picard [nondim] + type(group_pass_type) :: pass_visc_and_newton !< Handle for Newton-and-viscosity-related group passes + type(group_pass_type) :: pass_newton !< Handle for Newton-related group passes logical :: newton_adapt_cg_tol !< Use an adaptive CG tolerance during Newton iterations + real :: ew_gamma !< Gamma in Eisenstat-Walker adaptive Newton tolerance [nondim]. + real :: ew_alpha !< Alpha in Eisenstat-Walker adaptive Newton tolerance [nondim]. + integer :: ew_safety !< Safeguard Eisenstat-Walker using: + !!(0) no safeguard, (1) EW choice 2 threshold or (2) PETSc option 3 (Chacon 2008) + real :: ew_1_thres !< Threshold for Eisenstat-Walker version 1 [nondim] + real :: ew_eta_max !< Maximum allowed Eisenstat-Walker eta [nondim] integer :: cg_max_iterations !< The maximum number of iterations that can be used in the CG solver - integer :: nonlin_solve_err_mode !< 1: exit vel solve based on nonlin residual + integer :: nonlin_solve_err_mode !< 1: exit based on nonlin residual | F | / | F_0 | where | | is infty-norm !! 2: exit based on "fixed point" metric (|u - u_last| / |u| < tol) where | | is infty-norm - !! 3: exit based on change of norm - + !! 3: exit based on change of solution norm 2*abs(|u|-|u_last|)/(|u|+|u_last|) where | | is L2-norm + !! 4: exit based on nonlin residual | F | / | F_0 | where | | is L2-norm + !! 5: exit based on relative residual | F | / | tau | where | | is L2-norm + logical :: ssa_add_rel_resid !< Nonlinear error in velocity solve will also depend on the + !! L2 residual norm relative to RHS norm + real :: rr_nonlinear_tolerance !< If ssa_add_rel_resid, the additional nonlin tolerance in the iterative + !! velocity solve used for the relative residual [nondim] ! for write_ice_shelf_energy type(time_type) :: energysavedays !< The interval between writing the energies !! and other integral quantities of the run. @@ -378,7 +393,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) call get_param(param_file, mdl, "NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS", CS%visc_qps, & "Number of ice viscosity quadrature points. Either 1 (cell-centered) for 4", & - units="none", default=1) + units="none", default=4) if (CS%visc_qps/=1 .and. CS%visc_qps/=4) call MOM_error (FATAL, & "NUMBER OF ICE_VISCOSITY_QUADRATURE_POINTS must be 1 or 4") @@ -401,13 +416,11 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%newton_str_ux(isd:ied,jsd:jed,CS%visc_qps), source=0.0) allocate(CS%newton_str_vy(isd:ied,jsd:jed,CS%visc_qps), source=0.0) allocate(CS%newton_str_sh(isd:ied,jsd:jed,CS%visc_qps), source=0.0) - allocate(CS%newton_umid(isd:ied,jsd:jed), source=0.0) - allocate(CS%newton_vmid(isd:ied,jsd:jed), source=0.0) - allocate(CS%newton_drag_coef(isd:ied,jsd:jed), source=0.0) allocate(CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25) ! [Pa-3 s-1] - allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0) ! [R Z L2 T-1 ~> kg s-1] allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10*US%Pa_to_RLZ_T2) ! Units of [R L Z T-2 (s m-1)^n_sliding ~> Pa (s m-1)^n_sliding] + allocate(CS%coef_prefactor(isd:ied,jsd:jed), source=0.0) + allocate(CS%fB_elem(isd:ied,jsd:jed), source=0.0) allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0) allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0) allocate(CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0) @@ -421,6 +434,16 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%v_face_mask_bdry(IsdB:iedB,JsdB:JedB), source=-2.0) allocate(CS%h_bdry_val(isd:ied,jsd:jed), source=0.0) + ! Create group pass handles + call create_group_pass(CS%pass_visc_and_newton, CS%ice_visc, G%domain) + call create_group_pass(CS%pass_visc_and_newton, CS%newton_str_sh, G%domain) + call create_group_pass(CS%pass_visc_and_newton, CS%newton_visc_factor, G%domain) + call create_group_pass(CS%pass_visc_and_newton, CS%newton_str_ux, CS%newton_str_vy, G%domain, TO_ALL, AGRID) + + call create_group_pass(CS%pass_newton, CS%newton_str_sh, G%domain) + call create_group_pass(CS%pass_newton, CS%newton_visc_factor, G%domain) + call create_group_pass(CS%pass_newton, CS%newton_str_ux, CS%newton_str_vy, G%domain, TO_ALL, AGRID) + ! additional restarts for ice shelf state call register_restart_field(CS%u_shelf, "u_shelf", .false., restart_CS, & "ice sheet/shelf u-velocity", & @@ -594,20 +617,55 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "CF_Max", CS%CF_Max, & "Coulomb friction maximum coefficient", & units="none", default=0.5, fail_if_missing=.false.) + ! Pre-compute Coulomb prefactor alpha = (q-1)^(q-1)/q^q for q=CF_PostPeak [nondim]. + ! Default is 1.0; only update when Coulomb is active and q /= 1. + ! Also store CS%coulomb_pp_n = CF_PostPeak/n_basal_fric [nondim] + if (CS%CoulombFriction) then + if (CS%CF_PostPeak /= 1.0) then + CS%alpha_coulomb = (CS%CF_PostPeak-1.0)**(CS%CF_PostPeak-1.0) / CS%CF_PostPeak**CS%CF_PostPeak + endif + CS%coulomb_pp_n = CS%CF_PostPeak/CS%n_basal_fric + endif call get_param(param_file, mdl, "DENSITY_ICE", CS%density_ice, & "A typical density of ice.", units="kg m-3", default=917.0, scale=US%kg_m3_to_R) + + ! Precompute commonly-used density ratios + CS%rhoi_rhow=CS%density_ice / CS%density_ocean_avg + CS%rhow_rhoi=CS%density_ocean_avg / CS%density_ice + call get_param(param_file, mdl, "CONJUGATE_GRADIENT_TOLERANCE", CS%cg_tolerance, & - "tolerance in CG solver, relative to initial residual", units="nondim", default=1.e-6) - CS%cg_tol_newton = CS%cg_tolerance ! Will be tightened adaptively during Newton iterations + "For Picard iterations, the tolerance in CG solver, relative to initial residual", & + units="nondim", default=1.e-6) + call get_param(param_file, mdl, "NEWTON_CONJUGATE_GRADIENT_TOLERANCE", CS%cg_newton_tolerance, & + "For inexact Newton iterations, the initial tolerance in CG solver, relative to initial residual", & + units="nondim", default=CS%cg_tolerance) + CS%cg_tol_current = CS%cg_tolerance ! Can be tightened adaptively during inexact Newton iterations call get_param(param_file, mdl, "ICE_NONLINEAR_TOLERANCE", CS%nonlinear_tolerance, & "nonlin tolerance in iterative velocity solve", units="nondim", default=1.e-6) call get_param(param_file, mdl, "NEWTON_AFTER_TOLERANCE", CS%newton_after_tolerance, & "Switch from Picard to Newton iterations in the nonlinear ice velocity solve when "//& - "the fractional nonlinear residual falls below this tolerance.",& + "the fractional nonlinear residual falls below this tolerance. If <=0, no Picard.",& units="none", default=CS%nonlinear_tolerance) call get_param(param_file, mdl, "NEWTON_ADAPT_CG_TOL", CS%newton_adapt_cg_tol, & "Use an adaptive CG tolerance during Newton iterations.", default=.true.) + call get_param(param_file, mdl, "NEWTON_EW_GAMMA", CS%ew_gamma, & + "Gamma in Eisenstat-Walker adaptive Newton tolerance", units="nondim", default=0.9, & + do_not_log=(.not. CS%newton_adapt_cg_tol)) + call get_param(param_file, mdl, "NEWTON_EW_ALPHA", CS%ew_alpha, & + "Alpha in Eisenstat-Walker adaptive Newton tolerance", units="nondim", default=2.0, & + do_not_log=(.not. CS%newton_adapt_cg_tol)) + call get_param(param_file, mdl, "NEWTON_EW_SAFETY", CS%ew_safety, & + "Safeguard Eisenstat-Walker using (0) no safeguard, (1) EW choice 2 threshold "//& + "or (2) PETSc option 3 (Chacon 2008)", default=2, do_not_log=(.not. CS%newton_adapt_cg_tol)) + call get_param(param_file, mdl, "NEWTON_EW_1_THRESHOLD", CS%ew_1_thres, & + "Eisenstat-Walker version 1 threshold", & + units="nondim", default=0.1, do_not_log=(.not. CS%newton_adapt_cg_tol)) + call get_param(param_file, mdl, "NEWTON_EW_ETA_MAX", CS%ew_eta_max, & + "Maximum allowed Eisenstat-Walker eta (between 0 and 1)", & + units="nondim", default=0.9, do_not_log=(.not. CS%newton_adapt_cg_tol)) + if (CS%ew_eta_max<=0 .or. CS%ew_eta_max>= 1) & + call MOM_error(FATAL, "NEWTON_EW_ETA_MAX must be between 0 and 1.") call get_param(param_file, mdl, "ICE_SHELF_INNER_SOLVER", inner_solver_str, & "Choice of inner linear solver for the ice-shelf SSA velocity system. "//& "Valid choices are CG (default), CR, and MINRES.", & @@ -633,8 +691,19 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ units="m", default=1.e-3, scale=US%m_to_Z) call get_param(param_file, mdl, "NONLIN_SOLVE_ERR_MODE", CS%nonlin_solve_err_mode, & "Choose whether nonlin error in vel solve is based on nonlinear "//& - "residual (1), relative change since last iteration (2), or change in norm (3)", default=3) - + "Linf norm residual (1), Linf norm relative change since last iteration (2), "//& + "change in solution L2 norm (3), L2 norm residual (4), L2 backward norm (5)", default=3) + if (CS%nonlin_solve_err_mode /= 5) then + call get_param(param_file, mdl, "SSA_ADD_REL_RESID", CS%ssa_add_rel_resid, & + "Nonlinear error in vel solve will also depend on "// & + "L2 residual norm relative to RHS norm.", default=.false.) + else + CS%ssa_add_rel_resid = .false. !Avoids redundantly calculating err_mode 5 twice + endif + call get_param(param_file, mdl, "ICE_RR_NONLINEAR_TOLERANCE", CS%rr_nonlinear_tolerance, & + "if ssa_add_rel_resid, the additional nonlin tolerance "//& + "in the iterative velocity solve for the residual norm relative to RHS norm", & + units="nondim", default=1.e-4) call get_param(param_file, mdl, "SHELF_MOVING_FRONT", CS%moving_shelf_front, & "Specify whether to advance shelf front (and calve).", & default=.false.) @@ -766,7 +835,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%OD_av,G%domain, complete=.false.) call pass_var(CS%ground_frac, G%domain, complete=.false.) - call pass_var(CS%basal_traction, G%domain, complete=.false.) call pass_var(CS%AGlen_visc, G%domain, complete=.false.) call pass_var(CS%bed_elev, G%domain, complete=.false.) call pass_var(CS%C_basal_friction, G%domain, complete=.false.) @@ -956,8 +1024,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ !IS_dynamics_post_data is called before update_ice_shelf if (CS%id_taudx_shelf>0 .or. CS%id_taudy_shelf>0) & call calc_shelf_driving_stress(CS, ISS, G, US, CS%taudx_shelf, CS%taudy_shelf, CS%OD_av) - if (CS%id_taub>0) & - call calc_shelf_taub(CS, ISS, G, US, CS%u_shelf, CS%v_shelf) if (CS%id_visc_shelf>0) & call calc_shelf_visc(CS, ISS, G, US, CS%u_shelf, CS%v_shelf) endif @@ -978,17 +1044,15 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) type(time_type), intent(in) :: Time !< The current model time integer :: i, j, iters, isd, ied, jsd, jed - real :: rhoi_rhow real :: OD ! Depth of open water below the ice shelf [Z ~> m] type(time_type) :: dummy_time ! - rhoi_rhow = CS%density_ice / CS%density_ocean_avg dummy_time = set_time(0,0) isd=G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed do j=jsd,jed do i=isd,ied - OD = CS%bed_elev(i,j) - rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) + OD = CS%bed_elev(i,j) - CS%rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -1099,13 +1163,10 @@ subroutine volume_above_floatation(CS, G, ISS, vaf, hemisphere) real, dimension(SZI_(G),SZJ_(G)) :: vaf_cell !< cell-wise volume above floatation [Z L2 ~> m3] integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: is,ie,js,je,i,j - real :: rhoi_rhow, rhow_rhoi if (CS%GL_couple) & call MOM_error(FATAL, "MOM_ice_shelf_dyn, volume above floatation calculation assumes GL_couple=.FALSE..") - rhoi_rhow = CS%density_ice / CS%density_ocean_avg - rhow_rhoi = CS%density_ocean_avg / CS%density_ice is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec if (present(hemisphere)) then @@ -1135,7 +1196,7 @@ subroutine volume_above_floatation(CS, G, ISS, vaf, hemisphere) vaf_cell(i,j) = ISS%h_shelf(i,j) * ISS%area_shelf_h(i,j) else !grounded if vaf_cell(i,j) > 0 - vaf_cell(i,j) = max(ISS%h_shelf(i,j) - rhow_rhoi * CS%bed_elev(i,j), 0.0) * ISS%area_shelf_h(i,j) + vaf_cell(i,j) = max(ISS%h_shelf(i,j) - CS%rhow_rhoi * CS%bed_elev(i,j), 0.0) * ISS%area_shelf_h(i,j) endif endif enddo ; enddo @@ -1218,9 +1279,7 @@ subroutine IS_dynamics_post_data(time_step, Time, CS, ISS, G) call post_data(CS%id_visc_shelf, ice_visc, CS%diag) endif if (CS%id_taub > 0) then - do j=G%jsc,G%jec ; do i=G%isc,G%iec - basal_tr(i,j) = CS%basal_traction(i,j)*G%IareaT(i,j) - enddo ; enddo + call calc_shelf_taub(CS, ISS, G, basal_tr) call post_data(CS%id_taub, basal_tr, CS%diag) endif if (CS%id_u_mask > 0) call post_data(CS%id_u_mask, CS%umask, CS%diag) @@ -1530,16 +1589,22 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! If GL_regularize=true, indicates cells containing ! the grounding line (float_cond=1) or not (float_cond=0) real, dimension(SZDIB_(G),SZDJB_(G)) :: Normvec ! Velocities used for convergence [L2 T-2 ~> m2 s-2] + logical :: converged ! Indicates nonlinear convergence + logical :: calc_Au_for_convergence ! Used for convergence criteria than need a CG_action character(len=160) :: mesg ! The text of an error message integer :: conv_flag, i, j, k,l, iter, nodefloat integer :: Isdq, Iedq, Jsdq, Jedq, isd, ied, jsd, jed integer :: Iscq, Iecq, Jscq, Jecq, isc, iec, jsc, jec real :: err_max, err_tempu, err_tempv, err_init ! Errors in [R L3 Z T-2 ~> kg m s-2] or [L T-1 ~> m s-1] - real :: ew_prev_err ! Previous outer residual for Eisenstat-Walker CG tolerance (same units as err_max) + real :: norm_tau, err_rr ! Errors in [R L3 Z T-2 ~> kg m s-2] for relative residual + real :: ew_resid = 0.0 ! L2 norm of stress residual ||A(u)u - tau|| for Eisenstat-Walker [kg m s-2] + real :: ew_prev_resid = 0.0 ! Previous ew_resid; 0.0 flags first Newton call [kg m s-2] + real :: ew_eta = 0.0 ! Current EW inner tolerance [nondim] + real :: ew_eta_prev = 0.0 ! Previous EW inner tolerance for Chacon 2008 sharp-decrease safeguard [nondim] + real :: ew_stol ! Temporary safeguard tolerance [nondim] real :: max_vel ! The maximum velocity magnitude [L T-1 ~> m s-1] real :: tempu, tempv ! Temporary variables with velocity magnitudes [L T-1 ~> m s-1] real :: Norm, PrevNorm ! Velocities used to assess convergence [L T-1 ~> m s-1] - real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim] integer :: Is_sum, Js_sum, Ie_sum, Je_sum ! Loop bounds for global sums or arrays starting at 1. integer :: Iscq_sv, Jscq_sv ! Starting loop bound for sum_vec @@ -1547,7 +1612,22 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i Iscq = G%IscB ; Iecq = G%IecB ; Jscq = G%JscB ; Jecq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec - rhoi_rhow = CS%density_ice / CS%density_ocean_avg + + ! 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) + if (CS%nonlin_solve_err_mode >= 3 .or. CS%ssa_add_rel_resid) then + if ((isc+G%idg_offset==G%isg) .and. (.not. CS%reentrant_x)) then + Is_sum = Iscq + (1-Isdq) ; Iscq_sv = Iscq + else + Is_sum = isc + (1-Isdq) ; Iscq_sv = isc + endif + if ((jsc+G%jdg_offset==G%jsg) .and. (.not. CS%reentrant_y)) then + Js_sum = Jscq + (1-Jsdq) ; Jscq_sv = Jscq + else + Js_sum = jsc + (1-Jsdq) ; Jscq_sv = jsc + endif + Ie_sum = Iecq + (1-Isdq) ; Je_sum = Jecq + (1-Jsdq) + endif taudx(:,:) = 0.0 ; taudy(:,:) = 0.0 Au(:,:) = 0.0 ; Av(:,:) = 0.0 @@ -1558,22 +1638,25 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (.not. CS%GL_couple) then do j=G%jsc,G%jec ; do i=G%isc,G%iec - if (rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) > 0) then + if (CS%rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) > 0) then CS%ground_frac(i,j) = 1.0 CS%OD_av(i,j) =0.0 endif enddo ; enddo endif + ! Warning: This turns off Picard entirely and may not converge. + if (CS%newton_after_tolerance<=0.0) CS%doing_newton=.true. + + ! Calculate RHS 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) + ! This is to determine which cells contain the grounding line, the criterion being that the cell ! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by ! assuming topography is cellwise constant and H is bilinear in a cell; floating where ! rho_i/rho_w * H_node - D is negative - ! need to make this conditional on GL interp - if (CS%GL_regularize) then call interpolate_H_to_B(G, ISS%h_shelf, ISS%hmask, H_node, CS%min_h_shelf) @@ -1583,7 +1666,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i do l=0,1 ; do k=0,1 if ((ISS%hmask(i,j) == 1 .or. ISS%hmask(i,j)==3) .and. & - (rhoi_rhow * H_node(i-1+k,j-1+l) - CS%bed_elev(i,j) <= 0)) then + (CS%rhoi_rhow * H_node(i-1+k,j-1+l) - CS%bed_elev(i,j) <= 0)) then nodefloat = nodefloat + 1 endif enddo ; enddo @@ -1594,35 +1677,30 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i enddo ; enddo call pass_var(CS%float_cond, G%Domain, complete=.false.) - call pass_var(CS%ground_frac, G%domain, complete=.false.) + call pass_var(CS%ground_frac, G%domain, complete=.true.) endif - call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%basal_traction, G%domain, complete=.true.) + ! Calculate basal drag constants and initial velocity + call calc_shelf_basal_prefactors(CS, ISS, G, US) call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%ice_visc, G%domain) - - ! This makes sure basal stress is only applied when it is supposed to be - if (CS%GL_regularize) then - do j=G%jsd,G%jed ; do i=G%isd,G%ied - if (CS%ground_frac(i,j)/=1.0) CS%basal_traction(i,j) = 0.0 - enddo ; enddo + if (CS%doing_newton) then + ! halo pass for ice_visc, newton_str_sh, newton_visc_factor, newton_str_x + call do_group_pass(CS%pass_visc_and_newton, G%domain) 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 + call pass_var(CS%ice_visc, G%domain, complete=.true.) endif - if (CS%nonlin_solve_err_mode == 1) then - + ! Calculate err_init, the denominator for some convergence criteria + if (CS%nonlin_solve_err_mode == 1 .or. CS%nonlin_solve_err_mode == 4) then Au(:,:) = 0.0 ; Av(:,:) = 0.0 - call CG_action(CS, Au, Av, u_shlf, v_shlf, CS%Phi, CS%Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & - CS%ice_visc, CS%float_cond, CS%bed_elev, CS%basal_traction, & - G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow, use_newton_in=.false.) - call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) + CS%ice_visc, CS%float_cond, CS%bed_elev, u_shlf, v_shlf, & + G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, CS%rhoi_rhow, use_newton_in=.false.) + call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) ! TODO: is this needed? + endif + if (CS%nonlin_solve_err_mode == 1) then err_init = 0 ; err_tempu = 0 ; err_tempv = 0 do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB if (CS%umask(I,J) == 1) then @@ -1634,39 +1712,57 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (err_tempv >= err_init) err_init = err_tempv endif enddo ; enddo - call max_across_PEs(err_init) + elseif (CS%nonlin_solve_err_mode == 3) then Normvec(:,:) = 0.0 + do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq + if (CS%umask(I,J) == 1) Normvec(I,J) = (u_shlf(I,J)**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (v_shlf(I,J)**2) + enddo ; enddo + Norm = sqrt( reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, unscale=US%L_T_to_m_s**2 ) ) - ! 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) - if ((isc+G%idg_offset==G%isg) .and. (.not. CS%reentrant_x)) then - Is_sum = Iscq + (1-Isdq) ; Iscq_sv = Iscq - else - Is_sum = isc + (1-Isdq) ; Iscq_sv = isc - endif - if ((jsc+G%jdg_offset==G%jsg) .and. (.not. CS%reentrant_y)) then - Js_sum = Jscq + (1-Jsdq) ; Jscq_sv = Jscq - else - Js_sum = jsc + (1-Jsdq) ; Jscq_sv = jsc - endif - Ie_sum = Iecq + (1-Isdq) ; Je_sum = Jecq + (1-Jsdq) + elseif (CS%nonlin_solve_err_mode == 4) then + Normvec(:,:) = 0.0 + do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq + if (CS%umask(I,J) == 1) Normvec(I,J) = ((Au(I,J) - taudx(I,J))**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + ((Av(I,J) - taudy(I,J))**2) + enddo ; enddo + err_init = sqrt(reproducing_sum(Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, & + unscale=((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2)) + endif + if (CS%nonlin_solve_err_mode == 5 .or. CS%ssa_add_rel_resid) then + Normvec(:,:) = 0.0 do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq - if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + u_shlf(I,J)**2 - if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + v_shlf(I,J)**2 + if (CS%umask(I,J) == 1) Normvec(I,J) = (taudx(I,J)**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (taudy(I,J)**2) enddo ; enddo - Norm = sqrt( reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, unscale=US%L_T_to_m_s**2 ) ) + if (CS%nonlin_solve_err_mode == 5) then + err_init = sqrt(reproducing_sum(Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, & + unscale=((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2)) + else + norm_tau = sqrt(reproducing_sum(Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, & + unscale=((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2)) + endif endif u_last(:,:) = u_shlf(:,:) ; v_last(:,:) = v_shlf(:,:) - CS%cg_tol_newton = CS%cg_tolerance + if (CS%doing_newton) then + CS%cg_tol_current = CS%cg_newton_tolerance + else + CS%cg_tol_current = CS%cg_tolerance + endif + ew_prev_resid = 0.0 + converged = .false. + calc_Au_for_convergence = (CS%nonlin_solve_err_mode == 1 .or. CS%nonlin_solve_err_mode == 4 .or. & + CS%nonlin_solve_err_mode == 5 .or. CS%ssa_add_rel_resid) !! begin loop do iter=1,50 + ! The linear solve call ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H_node, CS%float_cond, & ISS%hmask, conv_flag, iters, time, CS%Phi, CS%Phisub) @@ -1678,53 +1774,58 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations" call MOM_mesg(mesg, 5) - call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%basal_traction, G%domain, complete=.true.) + ! Update viscosity call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%ice_visc, G%domain, complete=.false.) - call pass_var(CS%newton_str_sh, G%domain, complete=.false.) - call pass_var(CS%newton_visc_factor, G%domain, complete=.true.) - call pass_var(CS%newton_drag_coef, G%domain) - call pass_vector(CS%newton_str_ux, CS%newton_str_vy, G%domain, TO_ALL, AGRID) - call pass_vector(CS%newton_umid, CS%newton_vmid, G%domain, TO_ALL, AGRID) - - ! makes sure basal stress is only applied when it is supposed to be - if (CS%GL_regularize) then - do j=G%jsd,G%jed ; do i=G%isd,G%ied - if (CS%ground_frac(i,j)/=1.0) CS%basal_traction(i,j) = 0.0 - enddo ; enddo + + if (CS%doing_newton) then + ! halo pass for ice_visc, newton_str_sh, newton_visc_factor, newton_str_x + call do_group_pass(CS%pass_visc_and_newton, G%domain) 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 + call pass_var(CS%ice_visc, G%domain, complete=.true.) endif - if (CS%nonlin_solve_err_mode == 1) then - + ! Calculate convergence norms + if (calc_Au_for_convergence) then Au(:,:) = 0 ; Av(:,:) = 0 + call CG_action(CS, Au, Av, u_shlf, v_shlf, CS%Phi, CS%Phisub, CS%umask, CS%vmask, ISS%hmask, & + H_node, CS%ice_visc, CS%float_cond, CS%bed_elev, u_shlf, v_shlf, & + G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, CS%rhoi_rhow, use_newton_in=.false.) - call CG_action(CS, Au, Av, u_shlf, v_shlf, CS%Phi, CS%Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & - CS%ice_visc, CS%float_cond, CS%bed_elev, CS%basal_traction, & - G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow, use_newton_in=.false.) + if (CS%nonlin_solve_err_mode == 1) then + err_max = 0 - call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) + do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB + if (CS%umask(I,J) == 1) then + err_tempu = ABS(Au(I,J) - taudx(I,J)) + if (err_tempu >= err_max) err_max = err_tempu + endif + if (CS%vmask(I,J) == 1) then + err_tempv = ABS(Av(I,J) - taudy(I,J)) + if (err_tempv >= err_max) err_max = err_tempv + endif + enddo ; enddo - err_max = 0 + call max_across_PEs(err_max) + endif - do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB - if (CS%umask(I,J) == 1) then - err_tempu = ABS(Au(I,J) - taudx(I,J)) - if (err_tempu >= err_max) err_max = err_tempu - endif - if (CS%vmask(I,J) == 1) then - err_tempv = ABS(Av(I,J) - taudy(I,J)) - if (err_tempv >= err_max) err_max = err_tempv + if (CS%nonlin_solve_err_mode == 4 .or. CS%nonlin_solve_err_mode == 5 .or. CS%ssa_add_rel_resid) then + Normvec(:,:) = 0.0 + do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq + if (CS%umask(I,J) == 1) Normvec(I,J) = ((Au(I,J) - taudx(I,J))**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + ((Av(I,J) - taudy(I,J))**2) + enddo ; enddo + if (CS%nonlin_solve_err_mode == 4 .or. CS%nonlin_solve_err_mode == 5) then + err_max = sqrt(reproducing_sum(Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, & + unscale=((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2)) + if (CS%ssa_add_rel_resid) err_rr = err_max + elseif (CS%ssa_add_rel_resid) then + err_rr = sqrt(reproducing_sum(Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, & + unscale=((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2)) endif - enddo ; enddo - - call max_across_PEs(err_max) + endif + endif - elseif (CS%nonlin_solve_err_mode == 2) then + if (CS%nonlin_solve_err_mode == 2) then err_max=0. ; max_vel = 0 ; tempu = 0 ; tempv = 0 ; err_tempu = 0 do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB @@ -1753,43 +1854,121 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i elseif (CS%nonlin_solve_err_mode == 3) then PrevNorm = Norm ; Norm = 0.0 ; Normvec=0.0 do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq - if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + u_shlf(I,J)**2 - if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + v_shlf(I,J)**2 + if (CS%umask(I,J) == 1) Normvec(I,J) = (u_shlf(I,J)**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (v_shlf(I,J)**2) enddo ; enddo Norm = sqrt( reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, unscale=US%L_T_to_m_s**2 ) ) err_max = 2.*abs(Norm-PrevNorm) ; err_init = Norm+PrevNorm endif - write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init - call MOM_mesg(mesg, 5) + !Test convergence + if (err_max <= CS%nonlinear_tolerance * err_init) then + if (CS%ssa_add_rel_resid) then + if (err_rr <= CS%rr_nonlinear_tolerance * norm_tau) converged = .true. + else + converged = .true. + endif + endif - if (err_max <= CS%newton_after_tolerance * err_init .and. .not. CS%doing_newton) then - CS%doing_newton = .true. - ew_prev_err = err_max ! seed Eisenstat-Walker with residual at the Newton switch point - write(mesg,*) "ice_shelf_solve_outer: switching to Newton iterations at iter = ", iter + if (converged) then + exit + else + write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init call MOM_mesg(mesg, 5) - endif - ! Eisenstat-Walker Choice II (Eisenstat & Walker 1994): η_k = γ*(||F_k||/||F_{k-1}||)^α - ! with γ=0.9, α=2. Uses the ratio of consecutive outer residuals so that the CG - ! tolerance scales linearly with the current error (enabling quadratic outer convergence) - ! without over-tightening at later Newton steps. The first Newton step uses the standard - ! cg_tolerance (ratio = 1 on entry). - if (CS%doing_newton .and. CS%newton_adapt_cg_tol) then - CS%cg_tol_newton = min(CS%cg_tolerance, 0.9 * (err_max / ew_prev_err)**2) - ew_prev_err = err_max - endif + if (CS%ssa_add_rel_resid) then + write(mesg,*) "ice_shelf_solve_outer: nonlinear relative stress residual = ", err_rr/norm_tau + call MOM_mesg(mesg, 5) + endif - if (err_max <= CS%nonlinear_tolerance * err_init) then - exit - endif + ! Activate Newton + if (err_max <= CS%newton_after_tolerance * err_init .and. .not. CS%doing_newton) then + CS%doing_newton = .true. + write(mesg,*) "ice_shelf_solve_outer: switching to Newton iterations at iter = ", iter + call MOM_mesg(mesg, 7) + ! halo pass for newton_str_sh, newton_visc_factor, newton_str_x + call do_group_pass(CS%pass_newton, G%domain) + CS%cg_tol_current = CS%cg_newton_tolerance + endif + + ! Inexact Newton: Adapt inner solver tolerance to prevent oversolving + ! Based on Eisenstat-Walker Choice II (Eisenstat & Walker 1994): η_k = γ*(||F_k||/||F_{k-1}||)^α + ! with γ=0.9, α=2 as default. Uses the L2 norm of the nonlinear stress residual ||Au - tau||_2, + ! consistent with the inner solver's convergence check (sv3dsums(3)). + ! The first Newton step uses the standard cg_tolerance. + if (CS%doing_newton .and. CS%newton_adapt_cg_tol) then + !calculate residual needed for EW; some convergence criteria already did this + if (CS%nonlin_solve_err_mode >= 4) then + ew_resid=err_max + elseif (CS%ssa_add_rel_resid) then + ew_resid=err_rr + else + if (.not. calc_Au_for_convergence) then + Au(:,:) = 0 ; Av(:,:) = 0 + call CG_action(CS, Au, Av, u_shlf, v_shlf, CS%Phi, CS%Phisub, CS%umask, CS%vmask, ISS%hmask, & + H_node, CS%ice_visc, CS%float_cond, CS%bed_elev, u_shlf, v_shlf, & + G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, CS%rhoi_rhow, use_newton_in=.false.) + endif + Normvec(:,:) = 0.0 + do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq + if (CS%umask(I,J) == 1) Normvec(I,J) = ((Au(I,J) - taudx(I,J))**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + ((Av(I,J) - taudy(I,J))**2) + enddo ; enddo + ew_resid = sqrt(reproducing_sum(Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, & + unscale=((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2)) + endif + if (ew_prev_resid == 0.0) then + ! First Newton iteration: seed residuals; use initial newton cg_tolerance this step + ew_prev_resid = ew_resid + CS%cg_tol_current = CS%cg_newton_tolerance + ew_eta_prev = CS%cg_tol_current + else + ! Safeguarding and oversolving adjustments: + ! Eisenstat-Walker Choice II safeguard base formula + ew_eta = CS%ew_gamma * (ew_resid / ew_prev_resid)**CS%ew_alpha + ew_stol = CS%ew_gamma * ew_eta_prev**CS%ew_alpha + !Safeguards to sharp decrease/oversolving: + if (CS%ew_safety==1) then + ! Eisenstat-Walker Choice II safeguard: + write(mesg,*) "ice_shelf_solve_outer: ew_stol = ", ew_stol + call MOM_mesg(mesg, 8) + if (ew_stol > CS%ew_1_thres) ew_eta = max(ew_eta, ew_stol) + elseif (CS%ew_safety==2) then + ! PETSc choice 3 safeguard (e,g, Chacon 2008, J. Phys: Conf. Ser. 125 012041): + ! Avoid steep decreases in ew_eta + ew_eta = min(CS%cg_newton_tolerance, max(ew_eta, ew_stol)) + ! Avoid oversolving in last Newton iters: + ! The original is technically only applicable for nonlin_solve_err_mode=4: + ! ew_stol = CS%ew_gamma * ew_resid_first * CS%nonlinear_tolerance / ew_resid + ! Here, adapt for all nonlin_solve_err_modes: + ew_stol = CS%ew_gamma * err_init * CS%nonlinear_tolerance / err_max + if (CS%ssa_add_rel_resid) then + ew_stol = min(ew_stol, CS%ew_gamma * norm_tau * CS%rr_nonlinear_tolerance / err_rr) + endif + ew_eta = min(CS%cg_newton_tolerance, max(ew_eta, ew_stol)) + write(mesg,*) "ice_shelf_solve_outer: ew_stol = ", ew_stol + call MOM_mesg(mesg, 8) + endif + ew_eta = min(ew_eta,CS%ew_eta_max) + CS%cg_tol_current = ew_eta + ew_eta_prev = ew_eta + ew_prev_resid = ew_resid + write(mesg,*) "ice_shelf_solve_outer: New inner tolerance = ", CS%cg_tol_current + call MOM_mesg(mesg, 8) + endif + endif + endif enddo CS%doing_newton = .false. - CS%cg_tol_newton = CS%cg_tolerance + CS%cg_tol_current = CS%cg_tolerance write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init call MOM_mesg(mesg) + if (CS%ssa_add_rel_resid) then + write(mesg,*) "ice_shelf_solve_outer: nonlinear relative residual = ", err_rr/norm_tau + call MOM_mesg(mesg, 5) + endif write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" call MOM_mesg(mesg) @@ -1836,7 +2015,6 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H Au, Av, & ! Matrix-vector product A*x [R L3 Z T-2 ~> kg m s-2] DIAGu, DIAGv, & ! Diagonals [R L2 Z T-1 ~> kg s-1] IDIAGu, IDIAGv ! Reciprocal diagonals [R-1 L-2 Z-1 T ~> kg-1 s] - real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim] real :: resid_scale ! A scaling factor for redimensionalizing the global residuals ! [T3 kg m2 R-1 Z-1 L-4 s-3 ~> 1] integer :: Is_sum, Js_sum, Ie_sum, Je_sum ! Loop bounds for global sums or arrays starting at 1. @@ -1849,8 +2027,6 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H Iscq = G%IscB ; Iecq = G%IecB ; Jscq = G%JscB ; Jecq = G%JecB isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec - rhoi_rhow = CS%density_ice / CS%density_ocean_avg - ! Initialize shared arrays Au(:,:) = 0 ; Av(:,:) = 0 ; DIAGu(:,:) = 0 ; DIAGv(:,:) = 0 @@ -1871,13 +2047,13 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H RHSu(:,:) = taudx(:,:) ; RHSv(:,:) = taudy(:,:) call pass_vector(RHSu, RHSv, G%domain, TO_ALL, BGRID_NE, complete=.false.) - call matrix_diagonal(CS, G, US, float_cond, H_node, CS%ice_visc, CS%basal_traction, & - hmask, rhoi_rhow, Phi, Phisub, DIAGu, DIAGv) + call matrix_diagonal(CS, G, US, float_cond, H_node, CS%ice_visc, u_shlf, v_shlf, & + hmask, CS%rhoi_rhow, Phi, Phisub, DIAGu, DIAGv) call pass_vector(DIAGu, DIAGv, G%domain, TO_ALL, BGRID_NE, complete=.false.) call CG_action(CS, Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & - G, US, isc-1, iec+1, jsc-1, jec+1, rhoi_rhow, use_newton_in=.false.) + H_node, CS%ice_visc, float_cond, CS%bed_elev, u_shlf, v_shlf, & + G, US, isc-1, iec+1, jsc-1, jec+1, CS%rhoi_rhow, use_newton_in=.false.) call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE, complete=.true.) ! Precompute reciprocal diagonal @@ -1894,17 +2070,17 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H case (INNER_CG) call ice_shelf_solve_inner_CG(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, Av, & IDIAGu, IDIAGv, H_node, float_cond, hmask, & - rhoi_rhow, resid_scale, Phi, Phisub, conv_flag, iters, & + CS%rhoi_rhow, resid_scale, Phi, Phisub, conv_flag, iters, & Is_sum, Js_sum, Ie_sum, Je_sum, Iscq_sv, Jscq_sv) case (INNER_MINRES) call ice_shelf_solve_inner_MINRES(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, Av, & IDIAGu, IDIAGv, H_node, float_cond, hmask, & - rhoi_rhow, resid_scale, Phi, Phisub, conv_flag, iters, & + CS%rhoi_rhow, resid_scale, Phi, Phisub, conv_flag, iters, & Is_sum, Js_sum, Ie_sum, Je_sum, Iscq_sv, Jscq_sv) case (INNER_CR) call ice_shelf_solve_inner_CR(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, Av, & IDIAGu, IDIAGv, H_node, float_cond, hmask, & - rhoi_rhow, resid_scale, Phi, Phisub, conv_flag, iters, & + CS%rhoi_rhow, resid_scale, Phi, Phisub, conv_flag, iters, & Is_sum, Js_sum, Ie_sum, Je_sum, Iscq_sv, Jscq_sv) end select @@ -1977,6 +2153,10 @@ subroutine ice_shelf_solve_inner_CG(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A integer, intent(in) :: Iscq_sv !< Starting i-index for sum_vec arrays integer, intent(in) :: Jscq_sv !< Starting j-index for sum_vec arrays + real, dimension(SZDIB_(G),SZDJB_(G)) :: u_curr !< Frozen current iterate u^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] + real, dimension(SZDIB_(G),SZDJB_(G)) :: v_curr !< Frozen current iterate v^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] real, dimension(SZDIB_(G),SZDJB_(G)) :: & Ru, Rv, & ! Residuals [R L3 Z T-2 ~> m kg s-2] Zu, Zv, & ! Preconditioned residuals [L T-1 ~> m s-1] @@ -2013,6 +2193,9 @@ subroutine ice_shelf_solve_inner_CG(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A Ru(:,:) = (RHSu(:,:) - Au(:,:)) ; Rv(:,:) = (RHSv(:,:) - Av(:,:)) + ! current velocities used in CG_action for basal drag + u_curr(:,:) = u_shlf(:,:) ; v_curr(:,:) = v_shlf(:,:) + do J=Jsdq,Jedq ; do I=Isdq,Iedq if (CS%umask(I,J) == 1) Zu(I,J) = Ru(I,J) * IDIAGu(I,J) if (CS%vmask(I,J) == 1) Zv(I,J) = Rv(I,J) * IDIAGv(I,J) @@ -2037,7 +2220,7 @@ subroutine ice_shelf_solve_inner_CG(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A rho_old = sv3dsums(1) !resid0 = sqrt(sv3dsums(2)) - resid0tol2 = CS%cg_tol_newton**2 * sv3dsums(2) + resid0tol2 = CS%cg_tol_current**2 * sv3dsums(2) if (G%symmetric) then max_cg_halo=min(nx_halo,ny_halo) @@ -2070,7 +2253,7 @@ subroutine ice_shelf_solve_inner_CG(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A Au(:,:) = 0 ; Av(:,:) = 0 call CG_action(CS, Au, Av, Du, Dv, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & + H_node, CS%ice_visc, float_cond, CS%bed_elev, u_curr, v_curr, & G, US, is, ie, js, je, rhoi_rhow) sum_vec(:,:) = 0.0 @@ -2201,6 +2384,10 @@ subroutine ice_shelf_solve_inner_MINRES(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, A integer, intent(in) :: Iscq_sv !< Starting i-index for sum_vec arrays integer, intent(in) :: Jscq_sv !< Starting j-index for sum_vec arrays + real, dimension(SZDIB_(G),SZDJB_(G)) :: u_curr !< Frozen current iterate u^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] + real, dimension(SZDIB_(G),SZDJB_(G)) :: v_curr !< Frozen current iterate v^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] real, dimension(SZDIB_(G),SZDJB_(G)) :: & V_old_u, V_old_v, V_curr_u, V_curr_v, V_new_u, V_new_v, & ! Lanczos basis vectors [R L3 Z T-2 ~> m kg s-2] Z_curr_u, Z_curr_v, Z_new_u, Z_new_v, & ! Preconditioned Lanczos vectors [L T-1 ~> m s-1] @@ -2240,6 +2427,9 @@ subroutine ice_shelf_solve_inner_MINRES(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, A ! Initial Residual V_curr_u(:,:) = (RHSu(:,:) - Au(:,:)) ; V_curr_v(:,:) = (RHSv(:,:) - Av(:,:)) + ! current velocities used in CG_action for basal drag + u_curr(:,:) = u_shlf(:,:) ; v_curr(:,:) = v_shlf(:,:) + do J=Jscq,Jecq ; do I=Iscq,Iecq if (CS%umask(I,J) == 1) Z_curr_u(I,J) = V_curr_u(I,J) * IDIAGu(I,J) if (CS%vmask(I,J) == 1) Z_curr_v(I,J) = V_curr_v(I,J) * IDIAGv(I,J) @@ -2278,7 +2468,7 @@ subroutine ice_shelf_solve_inner_MINRES(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, A call pass_vector(Z_curr_u, Z_curr_v, G%domain, TO_ALL, BGRID_NE) eta = beta1 - resid0tol = CS%cg_tol_newton * beta1 + resid0tol = CS%cg_tol_current * beta1 conv_flag = 0 c0 = 1.0 ; s0 = 0.0 ; c1 = 1.0 ; s1 = 0.0 @@ -2294,7 +2484,7 @@ subroutine ice_shelf_solve_inner_MINRES(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, A ! --- STEP 1: Matrix Vector Product --- Qu(:,:) = 0 ; Qv(:,:) = 0 call CG_action(CS, Qu, Qv, Z_curr_u, Z_curr_v, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & + H_node, CS%ice_visc, float_cond, CS%bed_elev, u_curr, v_curr, & G, US, isc-1, iec+1, jsc-1, jec+1, rhoi_rhow) ! --- STEP 2: alpha = q dot z_curr --- sum_vec_3d(:,:) = 0.0 @@ -2437,6 +2627,10 @@ subroutine ice_shelf_solve_inner_CR(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A integer, intent(in) :: Iscq_sv !< Starting i-index for sum_vec arrays integer, intent(in) :: Jscq_sv !< Starting j-index for sum_vec arrays + real, dimension(SZDIB_(G),SZDJB_(G)) :: u_curr !< Frozen current iterate u^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] + real, dimension(SZDIB_(G),SZDJB_(G)) :: v_curr !< Frozen current iterate v^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] real, dimension(SZDIB_(G),SZDJB_(G)) :: & Ru, Rv, & ! Residuals (r) [R L3 Z T-2 ~> m kg s-2] Zu, Zv, & ! Preconditioned residuals (z = M^-1 r) [L T-1 ~> m s-1] @@ -2474,6 +2668,9 @@ subroutine ice_shelf_solve_inner_CR(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A ! r_0 = b - A*x_0 Ru(:,:) = (RHSu(:,:) - Au(:,:)) ; Rv(:,:) = (RHSv(:,:) - Av(:,:)) + ! current velocities used in CG_action for basal drag + u_curr(:,:) = u_shlf(:,:) ; v_curr(:,:) = v_shlf(:,:) + ! z_0 = M^-1 r_0 do J=Jsdq,Jedq ; do I=Isdq,Iedq if (CS%umask(I,J) == 1) Zu(I,J) = Ru(I,J) * IDIAGu(I,J) @@ -2486,7 +2683,7 @@ subroutine ice_shelf_solve_inner_CR(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A ! Compute A * z_0 Au(:,:) = 0 ; Av(:,:) = 0 call CG_action(CS, Au, Av, Zu, Zv, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & + H_node, CS%ice_visc, float_cond, CS%bed_elev, u_curr, v_curr, & G, US, isc-1, iec+1, jsc-1, jec+1, rhoi_rhow) call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) @@ -2510,7 +2707,7 @@ subroutine ice_shelf_solve_inner_CR(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A r_norm_sq = sv3dsums(1) z_w_sum = sv3dsums(2) - resid0tol2 = CS%cg_tol_newton**2 * r_norm_sq + resid0tol2 = CS%cg_tol_current**2 * r_norm_sq conv_flag = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -2564,7 +2761,7 @@ subroutine ice_shelf_solve_inner_CR(CS, G, US, u_shlf, v_shlf, RHSu, RHSv, Au, A ! --- STEP 3: w_{k+1} = A z_{k+1} --- Au(:,:) = 0 ; Av(:,:) = 0 call CG_action(CS, Au, Av, Zu, Zv, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & + H_node, CS%ice_visc, float_cond, CS%bed_elev, u_curr, v_curr, & G, US, isc-1, iec+1, jsc-1, jec+1, rhoi_rhow) call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) @@ -3062,7 +3259,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) real, dimension(SIZE(OD,1),SIZE(OD,2)) :: S ! surface elevation [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G)) :: sx_e, sy_e !element contributions to driving stress - real :: rho, rhow, rhoi_rhow ! Ice and ocean densities [R ~> kg m-3] + real :: rho, rhow ! Ice and ocean densities [R ~> kg m-3] real :: sx, sy ! Ice shelf top slopes at tracer points [Z L-1 ~> nondim] real :: neumann_val ! [R Z L2 T-2 ~> kg s-2] real :: grav ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] @@ -3087,7 +3284,6 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) rho = CS%density_ice rhow = CS%density_ocean_avg grav = CS%g_Earth - rhoi_rhow = rho/rhow ! prelim - go through and calculate S if (CS%GL_couple) then @@ -3097,8 +3293,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) else ! check whether the ice is floating or grounded do j=jsc-2,jec+2 ; do i=isc-2,iec+2 - if (rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) <= 0) then - S(i,j) = (1 - rhoi_rhow)*max(ISS%h_shelf(i,j),CS%min_h_shelf) + if (CS%rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) <= 0) then + S(i,j) = (1 - CS%rhoi_rhow)*max(ISS%h_shelf(i,j),CS%min_h_shelf) else S(i,j) = max(ISS%h_shelf(i,j),CS%min_h_shelf)-CS%bed_elev(i,j) endif @@ -3191,7 +3387,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (CS%ground_frac(i,j) == 1) then neumann_val = ((.5 * grav) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2 - rhow * CS%bed_elev(i,j)**2)) else - neumann_val = (.5 * grav) * ((1-rho/rhow) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2)) + neumann_val = (.5 * grav) * ((1-CS%rhoi_rhow) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2)) endif if ((CS%u_face_mask_bdry(I-1,j) == 2) .OR. & ((ISS%hmask(i-1,j) == 0 .OR. ISS%hmask(i-1,j) == 2) .AND. (CS%reentrant_x .OR. (i+i_off /= gisc)))) then @@ -3243,7 +3439,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) end subroutine calc_shelf_driving_stress subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmask, H_node, & - ice_visc, float_cond, bathyT, basal_trac, G, US, is, ie, js, je, dens_ratio, use_newton_in) + ice_visc, float_cond, bathyT, u_curr, v_curr, G, US, is, ie, js, je, dens_ratio, use_newton_in) type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. @@ -3277,14 +3473,17 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's !! flow law [R L4 Z T-1 ~> kg m2 s-1]. real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not + intent(in) :: float_cond !< If GL_regularize=true, indicates cells containing + !! the grounding line (float_cond=1) or not (float_cond=0) real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points !! relative to sea-level [Z ~> m]. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear - !! part of the "linearized" basal stress [R Z L2 T-1 ~> kg s-1]. + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(in) :: u_curr !< Frozen current iterate u^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(in) :: v_curr !< Frozen current iterate v^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] real, intent(in) :: dens_ratio !< The density of ice divided by the density !! of seawater, nondimensional @@ -3315,16 +3514,25 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, ! Phi_k is equal to 1 at vertex k, and 0 at vertex l /= k, and bilinear real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1] - real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1] - real :: strx_n, stry_n, strsh_n, dstrain_n, inner_dot_n ! Newton correction variables [T-1 ~> s-1], [T-2 ~> s-2] + real :: uq, vq ! Interpolated direction-vector δu at quadrature point [L T-1 ~> m s-1] + real :: strx_n, stry_n, strsh_n, dstrain_n ! Newton viscosity correction variables [T-1 ~> s-1], [T-2 ~> s-2] + real :: u_curr_qp, v_curr_qp ! Current iterate u^k at quadrature point [L T-1 ~> m s-1] + real :: unorm2_qp ! Regularized squared speed of u^k at quadrature point [L2 T-2 ~> m2 s-2] + real :: basal_coef_qp ! Picard basal friction coefficient at quadrature point [R L2 Z T-1 ~> kg s-1] + real :: drag_newt_qp ! Newton basal drag coefficient at quadrature point [R Z T-1 ~> kg m-2 s-1] + real :: inner_dot_qp ! u^k_qp · δu_qp inner product for Newton basal drag [L2 T-2 ~> m2 s-2] + real :: coef_prefactor_e ! Pre-computed area * C_basal_friction * L_T_to_m_s [R L2 Z T-1 ~> kg s-1] + real :: eps_vel2_e ! Velocity regularization squared for current element [L2 T-2 ~> m2 s-2] + real :: min_trac_e ! min_basal_traction * areaT for current element [R L2 Z T-1 ~> kg s-1] + real :: fB_e ! Pre-computed Coulomb fB for element; 0 for Weertman [(T L-1)^CF_PostPeak] real :: jac_wt ! Per-quadrature-point metric correction |J_q|/areaT [nondim] integer :: iq, jq, iphi, jphi, i, j, ilq, jlq, Itgt, Jtgt, qp, qpv logical :: visc_qp4 logical :: use_newton ! Whether to apply Newton tangent stiffness corrections logical :: do_newton_visc ! Whether to apply viscosity-related Newton tangent stiffness corrections real, dimension(2) :: xquad ! Nondimensional quadrature ratios [nondim] - real, dimension(2,2) :: Ucell, Vcell, Usub, Vsub ! Velocities at the nodal points around the cell [L T-1 ~> m s-1] - real, dimension(2,2) :: Hcell ! Ice shelf thickness at notal (corner) points [Z ~> m] + real, dimension(2,2) :: Usub, Vsub ! Subgrid nodal contributions to basal traction [R L3 Z T-2 ~> kg m s-2] + real, dimension(2,2) :: Hcell ! Ice shelf thickness at nodal (corner) points [Z ~> m] real, dimension(2,2,4) :: uret_qp, vret_qp ! Temporary arrays in [R Z L3 T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G),4) :: uret_b, vret_b ! Temporary arrays in [R Z L3 T-2 ~> kg m s-2] @@ -3348,6 +3556,13 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, uret_qp(:,:,:) = 0.0 ; vret_qp(:,:,:) = 0.0 + ! Pre-computed element-level basal friction quantities (updated each outer Newton iteration + ! by calc_shelf_basal_prefactors; avoids O(N_cg) recomputation of expensive prefactors). + coef_prefactor_e = CS%coef_prefactor(i,j) + eps_vel2_e = CS%eps_glen_min**2 * ((G%dxT(i,j)**2) + (G%dyT(i,j)**2)) + min_trac_e = CS%min_basal_traction * G%areaT(i,j) + fB_e = CS%fB_elem(i,j) ! 0 for Weertman; non-zero for Coulomb + do iq=1,2 ; do jq=1,2 qp = 2*(jq-1)+iq !current quad point @@ -3390,12 +3605,32 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, stry_n = CS%newton_str_vy(i,j,qpv) strsh_n = CS%newton_str_sh(i,j,qpv) dstrain_n = (((2.*strx_n + stry_n)*ux) + ((2.*stry_n + strx_n)*vy)) + & - strsh_n * (uy + vx) * 0.5 + (strsh_n * (uy + vx) * 0.5) endif - ! Newton correction for basal drag: compute inner_dot_n once per quadrature point - if (use_newton) then - inner_dot_n = (CS%newton_umid(i,j)*uq) + (CS%newton_vmid(i,j)*vq) + ! Basal friction and Newton Jacobian evaluated at this quadrature point (fully grounded cells only). + ! Evaluating at quadrature points rather than cell-averaged ensures the Newton correction is the + ! exact Jacobian of the Picard residual, enabling quadratic convergence for all friction exponents. + if (float_cond(i,j) == 0 .and. CS%ground_frac(i,j)>0) then + u_curr_qp = ((u_curr(I-1,J-1) * (xquad(3-iq) * xquad(3-jq))) + & + (u_curr(I,J) * (xquad(iq) * xquad(jq)))) + & + ((u_curr(I,J-1) * (xquad(iq) * xquad(3-jq))) + & + (u_curr(I-1,J) * (xquad(3-iq) * xquad(jq)))) + v_curr_qp = ((v_curr(I-1,J-1) * (xquad(3-iq) * xquad(3-jq))) + & + (v_curr(I,J) * (xquad(iq) * xquad(jq)))) + & + ((v_curr(I,J-1) * (xquad(iq) * xquad(3-jq))) + & + (v_curr(I-1,J) * (xquad(3-iq) * xquad(jq)))) + unorm2_qp = ((u_curr_qp**2) + (v_curr_qp**2)) + eps_vel2_e + call compute_basal_coef(unorm2_qp, coef_prefactor_e, min_trac_e, fB_e, & + CS%n_basal_fric, CS%CoulombFriction, CS%CF_PostPeak, US%L_T_to_m_s, use_newton, & + basal_coef_qp, drag_newt_qp) + ! Apply ground fraction scaling (replaces external scaling of basal_traction) + basal_coef_qp = basal_coef_qp * CS%ground_frac(i,j) + if (use_newton) then + drag_newt_qp = drag_newt_qp * CS%ground_frac(i,j) + ! Inner product u^k_qp . delta_u_qp for the Newton correction. + inner_dot_qp = (u_curr_qp * uq) + (v_curr_qp * vq) + endif endif ! Ratio |J_q|/areaT corrects the uniform-area weight baked into ice_visc for @@ -3410,7 +3645,7 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, (((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & ((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) - ! Newton tangent stiffness correction: add (dη/dε_e^2) * (g·δε) * (g·φ_m) term + ! Newton viscosity tangent stiffness: (dη/dε_e^2) * (g·δε) * (g·φ_m). if (do_newton_visc) then if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_n * & @@ -3422,19 +3657,21 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, ((2.*stry_n + strx_n) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) endif - if (float_cond(i,j) == 0) then + if (float_cond(i,j) == 0 .and. CS%ground_frac(i,j)>0) then ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 - if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & - (jac_wt * (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq))) - if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & - (jac_wt * (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq))) - ! Newton basal drag tangent stiffness: (m-1)*basal_trac/|u|^2 * u_i * (u . delta_u) contribution + ! Picard basal drag: C*|u^k|^(m-1) * δu evaluated at quadrature point, weighted by φ_m + if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & + (jac_wt * (basal_coef_qp * uq) * (xquad(ilq) * xquad(jlq))) + if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & + (jac_wt * (basal_coef_qp * vq) * (xquad(ilq) * xquad(jlq))) + ! Newton basal drag: pointwise Jacobian of the Picard residual. + ! Tangent stiffness = basal_coef_qp*I + drag_newt_qp * u^k_qp ⊗ u^k_qp if (use_newton) then if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & - jac_wt * CS%newton_drag_coef(i,j) * CS%newton_umid(i,j) * inner_dot_n * (xquad(ilq) * xquad(jlq)) + jac_wt * drag_newt_qp * u_curr_qp * inner_dot_qp * (xquad(ilq) * xquad(jlq)) if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & - jac_wt * CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j) * inner_dot_n * (xquad(ilq) * xquad(jlq)) + jac_wt * drag_newt_qp * v_curr_qp * inner_dot_qp * (xquad(ilq) * xquad(jlq)) endif endif enddo ; enddo @@ -3457,46 +3694,22 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, vret_b(I ,J ,1) = 0.25*((vret_qp(2,2,1)+vret_qp(2,2,4))+(vret_qp(2,2,2)+vret_qp(2,2,3))) if (float_cond(i,j) == 1) then - Ucell(:,:) = u_shlf(I-1:I,J-1:J) ; Vcell(:,:) = v_shlf(I-1:I,J-1:J) + ! Subgrid grounding-line: evaluate basal friction at each grounded sub-quadrature point. + ! Picard and Newton Jacobian are both computed inside CG_action_subgrid_basal. Hcell(:,:) = H_node(I-1:I,J-1:J) - - call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, & - bathyT(i,j), dens_ratio, Usub, Vsub, & + call CG_action_subgrid_basal(CS, G, US, Phisub, Hcell, & + u_curr(I-1:I,J-1:J), v_curr(I-1:I,J-1:J), & + u_shlf(I-1:I,J-1:J), v_shlf(I-1:I,J-1:J), & + bathyT(i,j), dens_ratio, i, j, fB_e, use_newton, Usub, Vsub, & G%dxCv(i,j-1), G%dxCv(i,j), G%dyCu(i-1,j), G%dyCu(i,j), G%IareaT(i,j)) - - if (umask(I-1,J-1) == 1) uret_b(I-1,J-1,4) = uret_b(I-1,J-1,4) + (Usub(1,1) * basal_trac(i,j)) - if (umask(I-1,J ) == 1) uret_b(I-1,J ,2) = uret_b(I-1,J ,2) + (Usub(1,2) * basal_trac(i,j)) - if (umask(I ,J-1) == 1) uret_b(I ,J-1,3) = uret_b(I ,J-1,3) + (Usub(2,1) * basal_trac(i,j)) - if (umask(I ,J ) == 1) uret_b(I ,J ,1) = uret_b(I ,J ,1) + (Usub(2,2) * basal_trac(i,j)) - - if (vmask(I-1,J-1) == 1) vret_b(I-1,J-1,4) = vret_b(I-1,J-1,4) + (Vsub(1,1) * basal_trac(i,j)) - if (vmask(I-1,J ) == 1) vret_b(I-1,J ,2) = vret_b(I-1,J ,2) + (Vsub(1,2) * basal_trac(i,j)) - if (vmask(I ,J-1) == 1) vret_b(I ,J-1,3) = vret_b(I ,J-1,3) + (Vsub(2,1) * basal_trac(i,j)) - if (vmask(I ,J ) == 1) vret_b(I ,J ,1) = vret_b(I ,J ,1) + (Vsub(2,2) * basal_trac(i,j)) - - ! Newton basal drag correction for subgrid grounding line cells. - ! inner_dot_sub(m,n) = sum over grounded sub-QPs of (u^k . delta_u) * phi_{m,n} * weight - ! = newton_umid * Usub(m,n) + newton_vmid * Vsub(m,n) - ! Correction to u-node (m,n): newton_drag_coef * newton_umid * inner_dot_sub(m,n) - ! Correction to v-node (m,n): newton_drag_coef * newton_vmid * inner_dot_sub(m,n) - if (use_newton) then - if (umask(I-1,J-1)==1) uret_b(I-1,J-1,4) = uret_b(I-1,J-1,4) + CS%newton_drag_coef(i,j) * & - CS%newton_umid(i,j) * ((CS%newton_umid(i,j)*Usub(1,1)) + (CS%newton_vmid(i,j)*Vsub(1,1))) - if (umask(I-1,J )==1) uret_b(I-1,J ,2) = uret_b(I-1,J ,2) + CS%newton_drag_coef(i,j) * & - CS%newton_umid(i,j) * ((CS%newton_umid(i,j)*Usub(1,2)) + (CS%newton_vmid(i,j)*Vsub(1,2))) - if (umask(I ,J-1)==1) uret_b(I ,J-1,3) = uret_b(I ,J-1,3) + CS%newton_drag_coef(i,j) * & - CS%newton_umid(i,j) * ((CS%newton_umid(i,j)*Usub(2,1)) + (CS%newton_vmid(i,j)*Vsub(2,1))) - if (umask(I ,J )==1) uret_b(I ,J ,1) = uret_b(I ,J ,1) + CS%newton_drag_coef(i,j) * & - CS%newton_umid(i,j) * ((CS%newton_umid(i,j)*Usub(2,2)) + (CS%newton_vmid(i,j)*Vsub(2,2))) - if (vmask(I-1,J-1)==1) vret_b(I-1,J-1,4) = vret_b(I-1,J-1,4) + CS%newton_drag_coef(i,j) * & - CS%newton_vmid(i,j) * ((CS%newton_umid(i,j)*Usub(1,1)) + (CS%newton_vmid(i,j)*Vsub(1,1))) - if (vmask(I-1,J )==1) vret_b(I-1,J ,2) = vret_b(I-1,J ,2) + CS%newton_drag_coef(i,j) * & - CS%newton_vmid(i,j) * ((CS%newton_umid(i,j)*Usub(1,2)) + (CS%newton_vmid(i,j)*Vsub(1,2))) - if (vmask(I ,J-1)==1) vret_b(I ,J-1,3) = vret_b(I ,J-1,3) + CS%newton_drag_coef(i,j) * & - CS%newton_vmid(i,j) * ((CS%newton_umid(i,j)*Usub(2,1)) + (CS%newton_vmid(i,j)*Vsub(2,1))) - if (vmask(I ,J )==1) vret_b(I ,J ,1) = vret_b(I ,J ,1) + CS%newton_drag_coef(i,j) * & - CS%newton_vmid(i,j) * ((CS%newton_umid(i,j)*Usub(2,2)) + (CS%newton_vmid(i,j)*Vsub(2,2))) - endif + if (umask(I-1,J-1) == 1) uret_b(I-1,J-1,4) = uret_b(I-1,J-1,4) + Usub(1,1) + if (umask(I-1,J ) == 1) uret_b(I-1,J ,2) = uret_b(I-1,J ,2) + Usub(1,2) + if (umask(I ,J-1) == 1) uret_b(I ,J-1,3) = uret_b(I ,J-1,3) + Usub(2,1) + if (umask(I ,J ) == 1) uret_b(I ,J ,1) = uret_b(I ,J ,1) + Usub(2,2) + if (vmask(I-1,J-1) == 1) vret_b(I-1,J-1,4) = vret_b(I-1,J-1,4) + Vsub(1,1) + if (vmask(I-1,J ) == 1) vret_b(I-1,J ,2) = vret_b(I-1,J ,2) + Vsub(1,2) + if (vmask(I ,J-1) == 1) vret_b(I ,J-1,3) = vret_b(I ,J-1,3) + Vsub(2,1) + if (vmask(I ,J ) == 1) vret_b(I ,J ,1) = vret_b(I ,J ,1) + Vsub(2,2) endif endif ; enddo ; enddo @@ -3507,85 +3720,188 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, end subroutine CG_action -subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, Vcontr, & +!> Compute subgrid grounding-line basal traction nodal contributions for a CG action. +!! Evaluates basal friction (Picard and Newton Jacobian) at each grounded sub-quadrature point. +!! The sub-qp flotation test accounts for partial grounding; no external ground_frac scaling needed. +subroutine CG_action_subgrid_basal(CS, G, US, Phisub, H, U_curr, V_curr, U_delta, V_delta, & + bathyT, dens_ratio, i_elem, j_elem, fB_e, use_newton, Ucontr, Vcontr, & dxCv_S, dxCv_N, dyCu_W, dyCu_E, IareaT) - real, dimension(:,:,:,:,:,:), & - intent(in) :: Phisub !< Quadrature structure weights at subgridscale - !! locations for finite element calculations [nondim] - real, dimension(2,2), intent(in) :: H !< The ice shelf thickness at nodal (corner) points [Z ~> m]. - real, dimension(2,2), intent(in) :: U !< The zonal ice shelf velocity at vertices [L T-1 ~> m s-1] - real, dimension(2,2), intent(in) :: V !< The meridional ice shelf velocity at vertices [L T-1 ~> m s-1] - real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points - !! relative to sea-level [Z ~> m]. - real, intent(in) :: dens_ratio !< The density of ice divided by the density - !! of seawater [nondim] - real, dimension(2,2), intent(out) :: Ucontr !< The areal average of u-velocities where the ice shelf - !! is grounded, or 0 where it is floating [L T-1 ~> m s-1]. - real, dimension(2,2), intent(out) :: Vcontr !< The areal average of v-velocities where the ice shelf - !! is grounded, or 0 where it is floating [L T-1 ~> m s-1]. - real, intent(in) :: dxCv_S !< The cell width at the southern (v-point) edge [L ~> m] - real, intent(in) :: dxCv_N !< The cell width at the northern (v-point) edge [L ~> m] - real, intent(in) :: dyCu_W !< The cell height at the western (u-point) edge [L ~> m] - real, intent(in) :: dyCu_E !< The cell height at the eastern (u-point) edge [L ~> m] - real, intent(in) :: IareaT !< The inverse of the cell area at the tracer point [L-2 ~> m-2] - - real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: Ucontr_sub, Vcontr_sub ! The contributions to Ucontr and Vcontr - !! at each sub-cell - real, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: uloc_arr !The local sub-cell u-velocity [L T-1 ~> m s-1] - real, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: vloc_arr !The local sub-cell v-velocity [L T-1 ~> m s-1] - real, dimension(2,2) :: Ucontr_q, Vcontr_q !Contributions to a node from each quadrature point in a sub-grid cell - real :: jac_sub_wt ! Per-sub-cell-QP metric correction |J_sub|/areaT [nondim] - real :: a, d ! Interpolated cell-edge spacings at the sub-cell QP [L ~> m] - real :: subarea ! The fractional sub-cell area in reference space [nondim] - real :: hloc ! The local sub-cell ice thickness [Z ~> m] + type(ice_shelf_dyn_CS), intent(in) :: CS !< Ice shelf control structure + type(ocean_grid_type), intent(in) :: G !< The grid structure + type(unit_scale_type), intent(in) :: US !< Unit conversion factors + real, dimension(:,:,:,:,:,:), intent(in) :: Phisub !< Sub-grid quadrature weights [nondim] + real, dimension(2,2), intent(in) :: H !< Ice thickness at element corners [Z ~> m] + real, dimension(2,2), intent(in) :: U_curr !< Frozen u^k at element corners [L T-1 ~> m s-1] + real, dimension(2,2), intent(in) :: V_curr !< Frozen v^k at element corners [L T-1 ~> m s-1] + real, dimension(2,2), intent(in) :: U_delta !< Search direction δu at element corners [L T-1 ~> m s-1] + real, dimension(2,2), intent(in) :: V_delta !< Search direction δv at element corners [L T-1 ~> m s-1] + real, intent(in) :: bathyT !< Ocean bathymetry depth at tracer point [Z ~> m] + real, intent(in) :: dens_ratio !< Ice density / water density [nondim] + integer, intent(in) :: i_elem !< Tracer-grid i-index of the element + integer, intent(in) :: j_elem !< Tracer-grid j-index of the element + real, intent(in) :: fB_e !< Element Coulomb parameter fB; 0 for Weertman [(T L-1)^CF_PostPeak] + logical, intent(in) :: use_newton !< If true, include Newton basal drag correction + real, dimension(2,2), intent(out) :: Ucontr !< Nodal u-contributions with friction applied [R L3 Z T-2 ~> kg m s-2] + real, dimension(2,2), intent(out) :: Vcontr !< Nodal v-contributions with friction applied [R L3 Z T-2 ~> kg m s-2] + real, intent(in) :: dxCv_S !< The cell width at the southern (v-point) edge [L ~> m] + real, intent(in) :: dxCv_N !< The cell width at the northern (v-point) edge [L ~> m] + real, intent(in) :: dyCu_W !< The cell height at the western (u-point) edge [L ~> m] + real, intent(in) :: dyCu_E !< The cell height at the eastern (u-point) edge [L ~> m] + real, intent(in) :: IareaT !< The inverse of the cell area at the tracer point [L-2 ~> m-2] + + real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: Ucontr_sub, Vcontr_sub ! The contributions to Ucontr and Vcontr + !! at each sub-cell + real, dimension(2,2,2,2) :: U_qp_nd, V_qp_nd ! Per-qp nodal contributions (qx,qy,m,n) + ! accumulated then pair-summed for rotation invariance + real :: hloc ! Local sub-cell ice thickness [Z ~> m] + real :: u_curr_loc ! Frozen u^k interpolated to sub-qp [L T-1 ~> m s-1] + real :: v_curr_loc ! Frozen v^k interpolated to sub-qp [L T-1 ~> m s-1] + real :: u_delta_loc ! Search direction δu interpolated to sub-qp [L T-1 ~> m s-1] + real :: v_delta_loc ! Search direction δv interpolated to sub-qp [L T-1 ~> m s-1] + real :: unorm2_loc ! Regularized |u^k|^2 at sub-qp [L2 T-2 ~> m2 s-2] + real :: basal_coef_loc ! Picard friction coefficient at sub-qp [R L2 Z T-1 ~> kg s-1] + real :: drag_newt_loc ! Newton drag coefficient at sub-qp [R Z T ~> kg m-2 s] + real :: inner_dot_loc ! u^k · δu inner product at sub-qp [L2 T-2 ~> m2 s-2] + real :: phi_mn ! Basis function value at sub-qp [nondim] + real :: contrib ! Quadrature weight contribution [nondim] + real :: coef_prefactor ! Pre-computed area * C_basal_friction * L_T_to_m_s [R L2 Z T-1 ~> kg s-1] + real :: min_trac_area ! Minimum area-integrated traction floor [R L2 Z T-1 ~> kg s-1] + real :: eps_vel2 ! Velocity regularization squared [L2 T-2 ~> m2 s-2] + real :: jac_sub_wt ! Per-sub-cell-QP metric correction |J_sub|/areaT [nondim] + real :: a, d ! Interpolated cell-edge spacings at the sub-cell QP [L ~> m] + real :: subarea ! Fractional sub-cell area [nondim] integer :: nsub, i, j, qx, qy, m, n - nsub = size(Phisub,3) + nsub = size(Phisub, 3) subarea = 1.0 / real(nsub)**2 - uloc_arr(:,:,:,:) = 0.0 ; vloc_arr(:,:,:,:)=0.0 + coef_prefactor = CS%coef_prefactor(i_elem,j_elem) + min_trac_area = CS%min_basal_traction * G%areaT(i_elem,j_elem) + eps_vel2 = CS%eps_glen_min**2 * ((G%dxT(i_elem,j_elem)**2) + (G%dyT(i_elem,j_elem)**2)) - do j=1,nsub ; do i=1,nsub ; do qy=1,2 ; do qx=1,2 - hloc = ((Phisub(qx,qy,i,j,1,1)*H(1,1)) + (Phisub(qx,qy,i,j,2,2)*H(2,2))) + & - ((Phisub(qx,qy,i,j,1,2)*H(1,2)) + (Phisub(qx,qy,i,j,2,1)*H(2,1))) - if (dens_ratio * hloc - bathyT > 0) then - uloc_arr(qx,qy,i,j) = (((Phisub(qx,qy,i,j,1,1) * U(1,1)) + (Phisub(qx,qy,i,j,2,2) * U(2,2))) + & - ((Phisub(qx,qy,i,j,1,2) * U(1,2)) + (Phisub(qx,qy,i,j,2,1) * U(2,1)))) - vloc_arr(qx,qy,i,j) = (((Phisub(qx,qy,i,j,1,1) * V(1,1)) + (Phisub(qx,qy,i,j,2,2) * V(2,2))) + & - ((Phisub(qx,qy,i,j,1,2) * V(1,2)) + (Phisub(qx,qy,i,j,2,1) * V(2,1)))) - endif - enddo ; enddo ; enddo ; enddo + Ucontr_sub(:,:,:,:) = 0.0 ; Vcontr_sub(:,:,:,:) = 0.0 - do n=1,2 ; do m=1,2 ; do j=1,nsub ; do i=1,nsub + do j=1,nsub ; do i=1,nsub + U_qp_nd(:,:,:,:) = 0.0 ; V_qp_nd(:,:,:,:) = 0.0 do qy=1,2 ; do qx=1,2 - ! Interpolate cell-edge metrics to the sub-cell QP using the bilinear shape function values - ! from bilinear_shape_functions_subgrid. Marginal sums of Phisub give the interpolation - ! weights: sum over k=1 nodes gives (1-y); k=2 gives y; l=1 gives (1-x); l=2 gives x. - ! This is analogous to jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) in the regular routines. - a = dxCv_S * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,2,1)) + & ! (1-y) * dxCv_S - dxCv_N * (Phisub(qx,qy,i,j,1,2) + Phisub(qx,qy,i,j,2,2)) ! + y * dxCv_N - d = dyCu_W * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,1,2)) + & ! (1-x) * dyCu_W - dyCu_E * (Phisub(qx,qy,i,j,2,1) + Phisub(qx,qy,i,j,2,2)) ! + x * dyCu_E - jac_sub_wt = 0.25 * subarea * (a * d) * IareaT - - !calculate quadrature point contributions for the sub-cell, to each node - Ucontr_q(qx,qy) = jac_sub_wt * Phisub(qx,qy,i,j,m,n) * uloc_arr(qx,qy,i,j) - Vcontr_q(qx,qy) = jac_sub_wt * Phisub(qx,qy,i,j,m,n) * vloc_arr(qx,qy,i,j) + hloc = ((Phisub(qx,qy,i,j,1,1)*H(1,1)) + (Phisub(qx,qy,i,j,2,2)*H(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*H(1,2)) + (Phisub(qx,qy,i,j,2,1)*H(2,1))) + if (dens_ratio * hloc - bathyT > 0) then ! grounded sub-qp + u_curr_loc = (((Phisub(qx,qy,i,j,1,1)*U_curr(1,1)) + (Phisub(qx,qy,i,j,2,2)*U_curr(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*U_curr(1,2)) + (Phisub(qx,qy,i,j,2,1)*U_curr(2,1)))) + v_curr_loc = (((Phisub(qx,qy,i,j,1,1)*V_curr(1,1)) + (Phisub(qx,qy,i,j,2,2)*V_curr(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*V_curr(1,2)) + (Phisub(qx,qy,i,j,2,1)*V_curr(2,1)))) + u_delta_loc = (((Phisub(qx,qy,i,j,1,1)*U_delta(1,1)) + (Phisub(qx,qy,i,j,2,2)*U_delta(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*U_delta(1,2)) + (Phisub(qx,qy,i,j,2,1)*U_delta(2,1)))) + v_delta_loc = (((Phisub(qx,qy,i,j,1,1)*V_delta(1,1)) + (Phisub(qx,qy,i,j,2,2)*V_delta(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*V_delta(1,2)) + (Phisub(qx,qy,i,j,2,1)*V_delta(2,1)))) + + unorm2_loc = ((u_curr_loc**2) + (v_curr_loc**2)) + eps_vel2 + call compute_basal_coef(unorm2_loc, coef_prefactor, min_trac_area, fB_e, & + CS%n_basal_fric, CS%CoulombFriction, CS%CF_PostPeak, US%L_T_to_m_s, use_newton, & + basal_coef_loc, drag_newt_loc) + inner_dot_loc = (u_curr_loc * u_delta_loc) + (v_curr_loc * v_delta_loc) + + ! Interpolate cell-edge metrics to the sub-cell QP using the bilinear shape function values + ! from bilinear_shape_functions_subgrid. Marginal sums of Phisub give the interpolation + ! weights: sum over k=1 nodes gives (1-y); k=2 gives y; l=1 gives (1-x); l=2 gives x. + ! This is analogous to jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) in the regular routines. + a = (dxCv_S * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,2,1))) + & ! (1-y) * dxCv_S + (dxCv_N * (Phisub(qx,qy,i,j,1,2) + Phisub(qx,qy,i,j,2,2))) ! + y * dxCv_N + d = (dyCu_W * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,1,2))) + & ! (1-x) * dyCu_W + (dyCu_E * (Phisub(qx,qy,i,j,2,1) + Phisub(qx,qy,i,j,2,2))) ! + x * dyCu_E + jac_sub_wt = 0.25 * subarea * (a * d) * IareaT + + do n=1,2 ; do m=1,2 + phi_mn = Phisub(qx,qy,i,j,m,n) + contrib = jac_sub_wt * phi_mn + ! Picard: friction matrix applied to search direction δu + U_qp_nd(qx,qy,m,n) = contrib * (basal_coef_loc * u_delta_loc) + V_qp_nd(qx,qy,m,n) = contrib * (basal_coef_loc * v_delta_loc) + ! Newton: Jacobian d(tau_b_i)/d(u_j) = basal_coef*I + drag_newt*u^k_i*u^k_j + if (use_newton) then + U_qp_nd(qx,qy,m,n) = U_qp_nd(qx,qy,m,n) + (contrib * (drag_newt_loc * u_curr_loc * inner_dot_loc)) + V_qp_nd(qx,qy,m,n) = V_qp_nd(qx,qy,m,n) + (contrib * (drag_newt_loc * v_curr_loc * inner_dot_loc)) + endif + enddo ; enddo + endif enddo ; enddo - !calculate sub-cell contribution to each node by summing up quadrature point contributions from the sub-cell - Ucontr_sub(i,j,m,n) = (Ucontr_q(1,1) + Ucontr_q(2,2)) + (Ucontr_q(1,2)+Ucontr_q(2,1)) - Vcontr_sub(i,j,m,n) = (Vcontr_q(1,1) + Vcontr_q(2,2)) + (Vcontr_q(1,2)+Vcontr_q(2,1)) - enddo ; enddo ; enddo ; enddo + do n=1,2 ; do m=1,2 + Ucontr_sub(i,j,m,n) = (U_qp_nd(1,1,m,n) + U_qp_nd(2,2,m,n)) + & + (U_qp_nd(1,2,m,n) + U_qp_nd(2,1,m,n)) + Vcontr_sub(i,j,m,n) = (V_qp_nd(1,1,m,n) + V_qp_nd(2,2,m,n)) + & + (V_qp_nd(1,2,m,n) + V_qp_nd(2,1,m,n)) + enddo ; enddo + enddo ; enddo - !sum up the sub-cell contributions to each node do n=1,2 ; do m=1,2 - call sum_square_matrix(Ucontr(m,n),Ucontr_sub(:,:,m,n),nsub) - call sum_square_matrix(Vcontr(m,n),Vcontr_sub(:,:,m,n),nsub) + call sum_square_matrix(Ucontr(m,n), Ucontr_sub(:,:,m,n), nsub) + call sum_square_matrix(Vcontr(m,n), Vcontr_sub(:,:,m,n), nsub) enddo ; enddo end subroutine CG_action_subgrid_basal +!> Compute the Picard basal friction coefficient and Newton drag coefficient at a +!! single quadrature point. Encapsulates the 3-path dispatch (linear Weertman / nonlinear +!! Weertman / Coulomb) so that CG_action, matrix_diagonal, and their subgrid equivalents +!! remain readable. The ground_frac scaling is NOT applied here; callers do it after the call. +subroutine compute_basal_coef(unorm2_qp, coef_prefactor, min_trac_area, fB_e, & + n_basal_fric, CoulombFriction, CF_PostPeak, L_T_to_m_s, use_newton, & + basal_coef, drag_newt) + real, intent(in) :: unorm2_qp !< Regularized |u^k|^2 > 0 at quadrature point [L2 T-2 ~> m2 s-2] + real, intent(in) :: coef_prefactor !< Pre-computed area * C_basal_friction * L_T_to_m_s [R L2 Z T-1 ~> kg s-1] + real, intent(in) :: min_trac_area !< Pre-computed min_basal_traction * areaT floor [R L2 Z T-1 ~> kg s-1] + real, intent(in) :: fB_e !< Element-level Coulomb fB; 0 for Weertman [(T L-1)^CF_PostPeak] + real, intent(in) :: n_basal_fric !< Friction sliding exponent m [nondim] + logical, intent(in) :: CoulombFriction !< True if using Coulomb friction + real, intent(in) :: CF_PostPeak !< Coulomb post-peak exponent q [nondim] + real, intent(in) :: L_T_to_m_s !< Unit conversion factor from internal [L T-1] to [m s-1] + logical, intent(in) :: use_newton !< If true, evaluate drag_newt; otherwise set to 0 + real, intent(out) :: basal_coef !< Picard friction coefficient at quadrature point [R L2 Z T-1 ~> kg s-1] + real, intent(out) :: drag_newt !< Newton drag coefficient [R Z T ~> kg m-2 s]; 0 without Newton + + real :: unorm ! |u^k| at quadrature point in physical units [m s-1] + real :: raw_coef ! Pre-floor friction coefficient [R L2 Z T-1 ~> kg s-1] + real :: fBuq ! fB_e * |u^k|^q [nondim] + + if (n_basal_fric == 1.0 .and. .not. CoulombFriction) then + ! Linear Weertman: coef is independent of |u|; sqrt and Newton correction not needed + basal_coef = max(coef_prefactor, min_trac_area) + drag_newt = 0.0 + elseif (CoulombFriction) then + ! Schoof/Gagliardini Coulomb friction + unorm = L_T_to_m_s * sqrt(unorm2_qp) + fBuq = fB_e * unorm**CF_PostPeak + raw_coef = coef_prefactor * (unorm**(n_basal_fric-1.0)) / (1.0 + fBuq)**n_basal_fric + if (raw_coef < min_trac_area) then + basal_coef = min_trac_area ; drag_newt = 0.0 + else + basal_coef = raw_coef + if (use_newton) then + drag_newt = (1.0/unorm2_qp) * raw_coef * & + ((n_basal_fric-1.0) - n_basal_fric * CF_PostPeak * fBuq / (1.0 + fBuq)) + else + drag_newt = 0.0 + endif + endif + else + ! Nonlinear Weertman (m > 1) + unorm = L_T_to_m_s * sqrt(unorm2_qp) + raw_coef = coef_prefactor * (unorm**(n_basal_fric-1.0)) + if (raw_coef < min_trac_area) then + basal_coef = min_trac_area ; drag_newt = 0.0 + else + basal_coef = raw_coef + if (use_newton) then + drag_newt = (n_basal_fric-1.0) / unorm2_qp * raw_coef + else + drag_newt = 0.0 + endif + endif + endif + +end subroutine compute_basal_coef !! Returns the sum of the elements in a square matrix. This sum is bitwise identical even if the matrices are rotated. subroutine sum_square_matrix(sum_out, mat_in, n) @@ -3627,8 +3943,8 @@ subroutine sum_square_matrix(sum_out, mat_in, n) end subroutine sum_square_matrix !> returns the diagonal entries of the matrix for a Jacobi preconditioning -subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, hmask, dens_ratio, & - Phi, Phisub, u_diagonal, v_diagonal) +subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, u_curr, v_curr, & + hmask, dens_ratio, Phi, Phisub, u_diagonal, v_diagonal) type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. @@ -3642,9 +3958,12 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, real, dimension(SZDI_(G),SZDJ_(G),CS%visc_qps), & intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's !! flow law [R L4 Z T-1 ~> kg m2 s-1]. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear - !! part of the "linearized" basal stress [R L3 T-1 ~> kg s-1]. + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(in) :: u_curr !< Frozen current iterate u^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(in) :: v_curr !< Frozen current iterate v^k, used to evaluate basal friction + !! at quadrature points [L T-1 ~> m s-1] real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf @@ -3666,13 +3985,20 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, ! returns the diagonal entries of the matrix for a Jacobi preconditioning real :: ux, uy, vx, vy ! Interpolated weight gradients [L-1 ~> m-1] - real :: uq, vq real :: jac_wt ! Per-quadrature-point metric correction |J_q|/areaT [nondim] real :: strx_n, stry_n, strsh_n ! Newton viscosity strain rates [T-1 ~> s-1] real :: dstrain_diag_u, dstrain_diag_v ! Newton viscosity diagonal correction factors [T-1 L-1 ~> s-1 m-1] real :: phi_m_sq ! Squared basis function value at quadrature point [nondim] + real :: u_curr_qp, v_curr_qp ! Current iterate u^k at quadrature point [L T-1 ~> m s-1] + real :: unorm2_qp ! Regularized squared speed of u^k at quadrature point [L2 T-2 ~> m2 s-2] + real :: basal_coef_qp ! Picard basal friction coefficient at quadrature point [R L2 Z T-1 ~> kg s-1] + real :: drag_newt_qp ! Newton basal drag coefficient at quadrature point [R Z T-1 ~> kg m-2 s-1] + real :: coef_prefactor_e ! Pre-computed area * C_basal_friction * L_T_to_m_s [R L2 Z T-1 ~> kg s-1] + real :: eps_vel2_e ! Velocity regularization squared for current element [L2 T-2 ~> m2 s-2] + real :: min_trac_e ! min_basal_traction * areaT for current element [R L2 Z T-1 ~> kg s-1] + real :: fB_e ! Pre-computed Coulomb fB for element; 0 for Weertman [(T L-1)^CF_PostPeak] real, dimension(2) :: xquad - real, dimension(2,2) :: Hcell, sub_ground + real, dimension(2,2) :: Hcell, u_diag_sub, v_diag_sub ! Subgrid diagonal contributions [R L2 Z T-1 ~> kg s-1] real, dimension(2,2,4) :: u_diag_qp, v_diag_qp real, dimension(SZDIB_(G),SZDJB_(G),4) :: u_diag_b, v_diag_b logical :: do_newton_visc ! Whether to apply viscosity-related Newton tangent stiffness corrections @@ -3702,6 +4028,12 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, u_diag_qp(:,:,:) = 0.0 ; v_diag_qp(:,:,:) = 0.0 + ! Pre-computed element-level basal friction quantities (updated each outer iteration). + coef_prefactor_e = CS%coef_prefactor(i,j) + eps_vel2_e = CS%eps_glen_min**2 * ((G%dxT(i,j)**2) + (G%dyT(i,j)**2)) + min_trac_e = CS%min_basal_traction * G%areaT(i,j) + fB_e = CS%fB_elem(i,j) ! 0 for Weertman; non-zero for Coulomb + do iq=1,2 ; do jq=1,2 qp = 2*(jq-1)+iq !current quad point @@ -3718,6 +4050,24 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, strsh_n = CS%newton_str_sh(i,j,qpv) endif + ! Basal friction coefficients at this quadrature point (fully grounded cells only) + if (float_cond(i,j) == 0 .and. CS%ground_frac(i,j)>0) then + u_curr_qp = ((u_curr(I-1,J-1) * (xquad(3-iq) * xquad(3-jq))) + & + (u_curr(I,J) * (xquad(iq) * xquad(jq)))) + & + ((u_curr(I,J-1) * (xquad(iq) * xquad(3-jq))) + & + (u_curr(I-1,J) * (xquad(3-iq) * xquad(jq)))) + v_curr_qp = ((v_curr(I-1,J-1) * (xquad(3-iq) * xquad(3-jq))) + & + (v_curr(I,J) * (xquad(iq) * xquad(jq)))) + & + ((v_curr(I,J-1) * (xquad(iq) * xquad(3-jq))) + & + (v_curr(I-1,J) * (xquad(3-iq) * xquad(jq)))) + unorm2_qp = ((u_curr_qp**2) + (v_curr_qp**2)) + eps_vel2_e + call compute_basal_coef(unorm2_qp, coef_prefactor_e, min_trac_e, fB_e, & + CS%n_basal_fric, CS%CoulombFriction, CS%CF_PostPeak, US%L_T_to_m_s, .true., & + basal_coef_qp, drag_newt_qp) + basal_coef_qp = basal_coef_qp * CS%ground_frac(i,j) + drag_newt_qp = drag_newt_qp * CS%ground_frac(i,j) + endif + do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi ilq = 1 ; if (iq == iphi) ilq = 2 @@ -3744,15 +4094,11 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_diag_u**2 endif - if (float_cond(i,j) == 0) then - uq = xquad(ilq) * xquad(jlq) - u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + & - jac_wt * (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)) - ! Newton basal drag diagonal correction: newton_drag_coef * (umid_i)^2 * phi_m^2 - if (CS%doing_newton) then - u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + & - jac_wt * CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * phi_m_sq - endif + if (float_cond(i,j) == 0 .and. CS%ground_frac(i,j)>0) then + ! Picard diagonal: basal_coef_qp * phi_m^2; Newton diagonal adds drag_newt_qp * u^k_qp^2 * phi_m^2. + u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + jac_wt * basal_coef_qp * phi_m_sq + if (CS%doing_newton) & + u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + jac_wt * drag_newt_qp * u_curr_qp**2 * phi_m_sq endif endif @@ -3775,15 +4121,10 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_diag_v**2 endif - if (float_cond(i,j) == 0) then - vq = xquad(ilq) * xquad(jlq) - v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + & - jac_wt * (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)) - ! Newton basal drag diagonal correction: newton_drag_coef * (vmid_i)^2 * phi_m^2 - if (CS%doing_newton) then - v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + & - jac_wt * CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * phi_m_sq - endif + if (float_cond(i,j) == 0 .and. CS%ground_frac(i,j)>0) then + v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + jac_wt * basal_coef_qp * phi_m_sq + if (CS%doing_newton) & + v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + jac_wt * drag_newt_qp * v_curr_qp**2 * phi_m_sq endif endif enddo ; enddo @@ -3806,41 +4147,23 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, v_diag_b(I ,J ,1) = 0.25*((v_diag_qp(2,2,1)+v_diag_qp(2,2,4))+(v_diag_qp(2,2,2)+v_diag_qp(2,2,3))) if (float_cond(i,j) == 1) then - Hcell(:,:) = H_node(i-1:i,j-1:j) - call CG_diagonal_subgrid_basal(Phisub, Hcell, CS%bed_elev(i,j), dens_ratio, sub_ground, & - G%dxCv(i,j-1), G%dxCv(i,j), G%dyCu(i-1,j), G%dyCu(i,j), G%IareaT(i,j)) - - if (CS%umask(I-1,J-1) == 1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + (sub_ground(1,1) * basal_trac(i,j)) - if (CS%umask(I-1,J ) == 1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + (sub_ground(1,2) * basal_trac(i,j)) - if (CS%umask(I ,J-1) == 1) u_diag_b(I ,J-1,3) = u_diag_b(I ,J-1,3) + (sub_ground(2,1) * basal_trac(i,j)) - if (CS%umask(I ,J ) == 1) u_diag_b(I ,J ,1) = u_diag_b(I ,J ,1) + (sub_ground(2,2) * basal_trac(i,j)) - - if (CS%vmask(I-1,J-1) == 1) v_diag_b(I-1,J-1,4) = v_diag_b(I-1,J-1,4) + (sub_ground(1,1) * basal_trac(i,j)) - if (CS%vmask(I-1,J ) == 1) v_diag_b(I-1,J ,2) = v_diag_b(I-1,J ,2) + (sub_ground(1,2) * basal_trac(i,j)) - if (CS%vmask(I ,J-1) == 1) v_diag_b(I ,J-1,3) = v_diag_b(I ,J-1,3) + (sub_ground(2,1) * basal_trac(i,j)) - if (CS%vmask(I ,J ) == 1) v_diag_b(I ,J ,1) = v_diag_b(I ,J ,1) + (sub_ground(2,2) * basal_trac(i,j)) - - ! Newton basal drag diagonal correction for subgrid grounding line cells. - ! sub_ground(m,n) = sum over grounded sub-QPs of phi_{m,n}^2 * weight, computed by - ! CG_diagonal_subgrid_basal. Newton diagonal = newton_drag_coef * umid^2 * sub_ground (for u-block). - if (CS%doing_newton) then - if (CS%umask(I-1,J-1)==1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + & - CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * sub_ground(1,1) - if (CS%umask(I-1,J )==1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + & - CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * sub_ground(1,2) - if (CS%umask(I ,J-1)==1) u_diag_b(I ,J-1,3) = u_diag_b(I ,J-1,3) + & - CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * sub_ground(2,1) - if (CS%umask(I ,J )==1) u_diag_b(I ,J ,1) = u_diag_b(I ,J ,1) + & - CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * sub_ground(2,2) - if (CS%vmask(I-1,J-1)==1) v_diag_b(I-1,J-1,4) = v_diag_b(I-1,J-1,4) + & - CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * sub_ground(1,1) - if (CS%vmask(I-1,J )==1) v_diag_b(I-1,J ,2) = v_diag_b(I-1,J ,2) + & - CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * sub_ground(1,2) - if (CS%vmask(I ,J-1)==1) v_diag_b(I ,J-1,3) = v_diag_b(I ,J-1,3) + & - CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * sub_ground(2,1) - if (CS%vmask(I ,J )==1) v_diag_b(I ,J ,1) = v_diag_b(I ,J ,1) + & - CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * sub_ground(2,2) - endif + ! Subgrid grounding-line: evaluate basal friction diagonal at each grounded sub-quadrature point. + ! Returns separate u_diag_sub and v_diag_sub (differ in Newton term: u^2 vs v^2). + ! The sub-qp flotation test handles grounding fraction; no external ground_frac scaling needed. + Hcell(:,:) = H_node(I-1:I,J-1:J) + call CG_diagonal_subgrid_basal(CS, G, US, Phisub, Hcell, & + u_curr(I-1:I,J-1:J), v_curr(I-1:I,J-1:J), & + CS%bed_elev(i,j), dens_ratio, i, j, fB_e, u_diag_sub, v_diag_sub, & + G%dxCv(i,j-1), G%dxCv(i,j), G%dyCu(i-1,j), G%dyCu(i,j), G%IareaT(i,j)) + + if (CS%umask(I-1,J-1)==1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + u_diag_sub(1,1) + if (CS%umask(I-1,J )==1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + u_diag_sub(1,2) + if (CS%umask(I ,J-1)==1) u_diag_b(I ,J-1,3) = u_diag_b(I ,J-1,3) + u_diag_sub(2,1) + if (CS%umask(I ,J )==1) u_diag_b(I ,J ,1) = u_diag_b(I ,J ,1) + u_diag_sub(2,2) + if (CS%vmask(I-1,J-1)==1) v_diag_b(I-1,J-1,4) = v_diag_b(I-1,J-1,4) + v_diag_sub(1,1) + if (CS%vmask(I-1,J )==1) v_diag_b(I-1,J ,2) = v_diag_b(I-1,J ,2) + v_diag_sub(1,2) + if (CS%vmask(I ,J-1)==1) v_diag_b(I ,J-1,3) = v_diag_b(I ,J-1,3) + v_diag_sub(2,1) + if (CS%vmask(I ,J )==1) v_diag_b(I ,J ,1) = v_diag_b(I ,J ,1) + v_diag_sub(2,2) endif endif ; enddo ; enddo @@ -3851,66 +4174,113 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, end subroutine matrix_diagonal -subroutine CG_diagonal_subgrid_basal(Phisub, H_node, bathyT, dens_ratio, f_grnd, & +!> Compute subgrid grounding-line basal traction contributions for the preconditioner diagonal. +!! Evaluates friction at each grounded sub-quadrature point. Returns separate u and v diagonals +!! because the Newton term uses u^2 for the u-block and v^2 for the v-block. +!! The sub-qp flotation test handles partial grounding; no external ground_frac scaling needed. +subroutine CG_diagonal_subgrid_basal(CS, G, US, Phisub, H_node, U_curr, V_curr, & + bathyT, dens_ratio, i_elem, j_elem, fB_e, u_diag, v_diag, & dxCv_S, dxCv_N, dyCu_W, dyCu_E, IareaT) - real, dimension(:,:,:,:,:,:), & - intent(in) :: Phisub !< Quadrature structure weights at subgridscale - !! locations for finite element calculations [nondim] - real, dimension(2,2), intent(in) :: H_node !< The ice shelf thickness at nodal (corner) - !! points [Z ~> m]. - real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m]. - real, intent(in) :: dens_ratio !< The density of ice divided by the density - !! of seawater [nondim] - real, dimension(2,2), intent(out) :: f_grnd !< The weighted fraction of the sub-cell where the ice shelf - !! is grounded [nondim] - real, intent(in) :: dxCv_S !< The cell width at the southern (v-point) edge [L ~> m] - real, intent(in) :: dxCv_N !< The cell width at the northern (v-point) edge [L ~> m] - real, intent(in) :: dyCu_W !< The cell height at the western (u-point) edge [L ~> m] - real, intent(in) :: dyCu_E !< The cell height at the eastern (u-point) edge [L ~> m] - real, intent(in) :: IareaT !< The inverse of the cell area at the tracer point [L-2 ~> m-2] - - real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: f_grnd_sub ! The contributions to nodal f_grnd - !! from each sub-cell - integer, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: grnd_stat !0 at floating quad points, 1 at grounded - real, dimension(2,2) :: f_grnd_q !Contributions to a node from each quadrature point in a sub-grid cell - real :: jac_sub_wt ! Per-sub-cell-QP metric correction |J_sub|/areaT [nondim] - real :: a, d ! Interpolated cell-edge spacings at the sub-cell QP [L ~> m] - real :: subarea ! The fractional sub-cell area in reference space [nondim] - real :: hloc ! The local sub-region thickness [Z ~> m] + type(ice_shelf_dyn_CS), intent(in) :: CS !< Ice shelf control structure + type(ocean_grid_type), intent(in) :: G !< The grid structure + type(unit_scale_type), intent(in) :: US !< Unit conversion factors + real, dimension(:,:,:,:,:,:), intent(in) :: Phisub !< Sub-grid quadrature weights [nondim] + real, dimension(2,2), intent(in) :: H_node !< Ice thickness at element corners [Z ~> m] + real, dimension(2,2), intent(in) :: U_curr !< Frozen u^k at element corners [L T-1 ~> m s-1] + real, dimension(2,2), intent(in) :: V_curr !< Frozen v^k at element corners [L T-1 ~> m s-1] + real, intent(in) :: bathyT !< Ocean bathymetry depth at tracer point [Z ~> m] + real, intent(in) :: dens_ratio !< Ice density / water density [nondim] + integer, intent(in) :: i_elem !< Tracer-grid i-index of the element + integer, intent(in) :: j_elem !< Tracer-grid j-index of the element + real, intent(in) :: fB_e !< Element Coulomb parameter fB; 0 for Weertman [(T L-1)^CF_PostPeak] + real, dimension(2,2), intent(out) :: u_diag !< Nodal u-diagonal entries [R L2 Z T-1 ~> kg s-1] + real, dimension(2,2), intent(out) :: v_diag !< Nodal v-diagonal entries [R L2 Z T-1 ~> kg s-1] + real, intent(in) :: dxCv_S !< The cell width at the southern (v-point) edge [L ~> m] + real, intent(in) :: dxCv_N !< The cell width at the northern (v-point) edge [L ~> m] + real, intent(in) :: dyCu_W !< The cell height at the western (u-point) edge [L ~> m] + real, intent(in) :: dyCu_E !< The cell height at the eastern (u-point) edge [L ~> m] + real, intent(in) :: IareaT !< The inverse of the cell area at the tracer point [L-2 ~> m-2] + + real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: u_diag_sub, v_diag_sub + real, dimension(2,2,2,2) :: u_diag_qp_nd, v_diag_qp_nd ! Per-qp nodal diagonal entries (qx,qy,m,n), + ! pair-summed for rotation invariance + real :: hloc ! Local sub-cell ice thickness [Z ~> m] + real :: u_curr_loc ! Frozen u^k interpolated to sub-qp [L T-1 ~> m s-1] + real :: v_curr_loc ! Frozen v^k interpolated to sub-qp [L T-1 ~> m s-1] + real :: unorm2_loc ! Regularized |u^k|^2 at sub-qp [L2 T-2 ~> m2 s-2] + real :: basal_coef_loc ! Picard friction coefficient at sub-qp [R L2 Z T-1 ~> kg s-1] + real :: drag_newt_loc ! Newton drag coefficient at sub-qp [R Z T ~> kg m-2 s] + real :: phi_mn_sq ! Squared basis function value at sub-qp [nondim] + real :: contrib ! Quadrature weight contribution [nondim] + real :: coef_prefactor ! Pre-computed area * C_basal_friction * L_T_to_m_s [R L2 Z T-1 ~> kg s-1] + real :: min_trac_area ! Minimum area-integrated traction floor [R L2 Z T-1 ~> kg s-1] + real :: eps_vel2 ! Velocity regularization squared [L2 T-2 ~> m2 s-2] + real :: jac_sub_wt ! Per-sub-cell-QP metric correction |J_sub|/areaT [nondim] + real :: a, d ! Interpolated cell-edge spacings at the sub-cell QP [L ~> m] + real :: subarea ! Fractional sub-cell area [nondim] integer :: nsub, i, j, qx, qy, m, n - nsub = size(Phisub,3) + nsub = size(Phisub, 3) subarea = 1.0 / real(nsub)**2 - grnd_stat(:,:,:,:)=0 - - do j=1,nsub ; do i=1,nsub ; do qy=1,2 ; do qx=1,2 - hloc = ((Phisub(qx,qy,i,j,1,1)*H_node(1,1)) + (Phisub(qx,qy,i,j,2,2)*H_node(2,2))) + & - ((Phisub(qx,qy,i,j,1,2)*H_node(1,2)) + (Phisub(qx,qy,i,j,2,1)*H_node(2,1))) - if (dens_ratio * hloc - bathyT > 0) grnd_stat(qx,qy,i,j) = 1 - enddo ; enddo ; enddo ; enddo - - do n=1,2 ; do m=1,2 ; do j=1,nsub ; do i=1,nsub - do qy=1,2 ; do qx = 1,2 - ! Interpolate cell-edge metrics to the sub-cell QP using the bilinear shape function values - ! from bilinear_shape_functions_subgrid. Marginal sums of Phisub give the interpolation - ! weights: sum over k=1 nodes gives (1-y); k=2 gives y; l=1 gives (1-x); l=2 gives x. - ! This is analogous to jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) in the regular routines. - a = dxCv_S * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,2,1)) + & ! (1-y) * dxCv_S - dxCv_N * (Phisub(qx,qy,i,j,1,2) + Phisub(qx,qy,i,j,2,2)) ! + y * dxCv_N - d = dyCu_W * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,1,2)) + & ! (1-x) * dyCu_W - dyCu_E * (Phisub(qx,qy,i,j,2,1) + Phisub(qx,qy,i,j,2,2)) ! + x * dyCu_E - jac_sub_wt = 0.25 * subarea * (a * d) * IareaT - - f_grnd_q(qx,qy) = jac_sub_wt * grnd_stat(qx,qy,i,j) * Phisub(qx,qy,i,j,m,n)**2 + coef_prefactor = CS%coef_prefactor(i_elem,j_elem) + min_trac_area = CS%min_basal_traction * G%areaT(i_elem,j_elem) + eps_vel2 = CS%eps_glen_min**2 * ((G%dxT(i_elem,j_elem)**2) + (G%dyT(i_elem,j_elem)**2)) + + u_diag_sub(:,:,:,:) = 0.0 ; v_diag_sub(:,:,:,:) = 0.0 + + do j=1,nsub ; do i=1,nsub + ! Zero the 4-qp per-node buffer so ungrounded qp contribute exactly 0. + u_diag_qp_nd(:,:,:,:) = 0.0 ; v_diag_qp_nd(:,:,:,:) = 0.0 + do qy=1,2 ; do qx=1,2 + hloc = ((Phisub(qx,qy,i,j,1,1)*H_node(1,1)) + (Phisub(qx,qy,i,j,2,2)*H_node(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*H_node(1,2)) + (Phisub(qx,qy,i,j,2,1)*H_node(2,1))) + if (dens_ratio * hloc - bathyT > 0) then ! grounded sub-qp + u_curr_loc = (((Phisub(qx,qy,i,j,1,1)*U_curr(1,1)) + (Phisub(qx,qy,i,j,2,2)*U_curr(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*U_curr(1,2)) + (Phisub(qx,qy,i,j,2,1)*U_curr(2,1)))) + v_curr_loc = (((Phisub(qx,qy,i,j,1,1)*V_curr(1,1)) + (Phisub(qx,qy,i,j,2,2)*V_curr(2,2))) + & + ((Phisub(qx,qy,i,j,1,2)*V_curr(1,2)) + (Phisub(qx,qy,i,j,2,1)*V_curr(2,1)))) + + unorm2_loc = ((u_curr_loc**2) + (v_curr_loc**2)) + eps_vel2 + call compute_basal_coef(unorm2_loc, coef_prefactor, min_trac_area, fB_e, & + CS%n_basal_fric, CS%CoulombFriction, CS%CF_PostPeak, US%L_T_to_m_s, .true., & + basal_coef_loc, drag_newt_loc) + ! Interpolate cell-edge metrics to the sub-cell QP using the bilinear shape function values + ! from bilinear_shape_functions_subgrid. Marginal sums of Phisub give the interpolation + ! weights: sum over k=1 nodes gives (1-y); k=2 gives y; l=1 gives (1-x); l=2 gives x. + ! This is analogous to jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) in the regular routines. + a = (dxCv_S * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,2,1))) + & ! (1-y) * dxCv_S + (dxCv_N * (Phisub(qx,qy,i,j,1,2) + Phisub(qx,qy,i,j,2,2))) ! + y * dxCv_N + d = (dyCu_W * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,1,2))) + & ! (1-x) * dyCu_W + (dyCu_E * (Phisub(qx,qy,i,j,2,1) + Phisub(qx,qy,i,j,2,2))) ! + x * dyCu_E + jac_sub_wt = 0.25 * subarea * (a * d) * IareaT + + do n=1,2 ; do m=1,2 + phi_mn_sq = Phisub(qx,qy,i,j,m,n)**2 + contrib = jac_sub_wt * phi_mn_sq + ! Picard diagonal + Newton diagonal (u_curr^2 for u-block, v_curr^2 for v-block) + if (CS%doing_newton) then + u_diag_qp_nd(qx,qy,m,n) = contrib * (basal_coef_loc + drag_newt_loc * u_curr_loc**2) + v_diag_qp_nd(qx,qy,m,n) = contrib * (basal_coef_loc + drag_newt_loc * v_curr_loc**2) + else + u_diag_qp_nd(qx,qy,m,n) = contrib * basal_coef_loc + v_diag_qp_nd(qx,qy,m,n) = contrib * basal_coef_loc + endif + enddo ; enddo + endif + enddo ; enddo + + do n=1,2 ; do m=1,2 + u_diag_sub(i,j,m,n) = (u_diag_qp_nd(1,1,m,n) + u_diag_qp_nd(2,2,m,n)) + & + (u_diag_qp_nd(1,2,m,n) + u_diag_qp_nd(2,1,m,n)) + v_diag_sub(i,j,m,n) = (v_diag_qp_nd(1,1,m,n) + v_diag_qp_nd(2,2,m,n)) + & + (v_diag_qp_nd(1,2,m,n) + v_diag_qp_nd(2,1,m,n)) enddo ; enddo - !calculate sub-cell contribution to each node by summing up quadrature point contributions from the sub-cell - f_grnd_sub(i,j,m,n) = (f_grnd_q(1,1) + f_grnd_q(2,2)) + (f_grnd_q(1,2)+f_grnd_q(2,1)) - enddo ; enddo ; enddo ; enddo + enddo ; enddo - !sum up the sub-cell contributions to each node do n=1,2 ; do m=1,2 - call sum_square_matrix(f_grnd(m,n),f_grnd_sub(:,:,m,n),nsub) + call sum_square_matrix(u_diag(m,n), u_diag_sub(:,:,m,n), nsub) + call sum_square_matrix(v_diag(m,n), v_diag_sub(:,:,m,n), nsub) enddo ; enddo end subroutine CG_diagonal_subgrid_basal @@ -4072,6 +4442,8 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) real :: Visc_coef, n_g real :: ux, uy, vx, vy real :: eps_min ! Velocity shears [T-1 ~> s-1] + real :: In_g ! inverse of Glen's exponent [nondim] + real :: eps_e2_exp ! (1.-n_g)/(2.*n_g) [nondim] logical :: model_qp1, model_qp4 isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec @@ -4093,6 +4465,8 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) endif n_g = CS%n_glen ; eps_min = CS%eps_glen_min + In_g=1./n_g + eps_e2_exp=(1.-n_g)/(2.*n_g) do j=jsc,jec ; do i=isc,iec @@ -4111,7 +4485,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) elseif (model_qp1) then ! calculate viscosity at 1 cell-centered quadrature point per cell - Visc_coef = (CS%AGlen_visc(i,j))**(-1./n_g) + Visc_coef = (CS%AGlen_visc(i,j))**(-In_g) ! Units of Aglen_visc [Pa-(n_g) s-1] ux = ((u_shlf(I-1,J-1) * CS%PhiC(1,i,j)) + & @@ -4136,24 +4510,23 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) CS%ice_visc(i,j,1) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & max(0.5 * Visc_coef * & - (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)) * & + (US%s_to_T**2 * (((ux**2) + (vy**2)) + ((ux*vy) + 0.25*((uy+vx)**2)) + eps_min**2))**(eps_e2_exp) * & (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) ! Rescale after the fractional power law. ! Store Newton tangent stiffness data: strain rates and coefficient for Newton iterations. - ! The Newton correction coefficient is (1/n-1)/2 * ice_visc / eps_e2, + ! The Newton correction coefficient is (1/n-1) * ice_visc / eps_e2, ! where eps_e2 = ux^2 + vy^2 + ux*vy + (uy+vx)^2/4 + eps_min^2 [T-2]. ! It is zero where ice_visc is limited by min_ice_visc (viscosity is not smooth there). CS%newton_str_ux(i,j,1) = ux ; CS%newton_str_vy(i,j,1) = vy CS%newton_str_sh(i,j,1) = uy + vx CS%newton_visc_factor(i,j,1) = 0.0 if (CS%ice_visc(i,j,1) > CS%min_ice_visc * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf))) then - CS%newton_visc_factor(i,j,1) = (0.5*(1./n_g - 1.) / & + CS%newton_visc_factor(i,j,1) = ((In_g - 1.) / & (((ux**2) + (vy**2)) + ((ux*vy) + 0.25*((uy+vx)**2)) + eps_min**2)) * & CS%ice_visc(i,j,1) endif elseif (model_qp4) then !calculate viscosity at 4 quadrature points per cell - - Visc_coef = (CS%AGlen_visc(i,j))**(-1./n_g) + Visc_coef = (CS%AGlen_visc(i,j))**(-In_g) do iq=1,2 ; do jq=1,2 @@ -4179,15 +4552,15 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) CS%ice_visc(i,j,2*(jq-1)+iq) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & max(0.5 * Visc_coef * & - (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)) * & + (US%s_to_T**2*(((ux**2) + (vy**2)) + ((ux*vy) + 0.25*((uy+vx)**2)) + eps_min**2))**(eps_e2_exp) * & (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) ! Rescale after the fractional power law. ! Store Newton tangent stiffness data at each quadrature point. CS%newton_str_ux(i,j,2*(jq-1)+iq) = ux ; CS%newton_str_vy(i,j,2*(jq-1)+iq) = vy - CS%newton_str_sh(i,j,2*(jq-1)+iq) = uy + vx + CS%newton_str_sh(i,j,2*(jq-1)+iq) = (uy + vx) CS%newton_visc_factor(i,j,2*(jq-1)+iq) = 0.0 if (CS%ice_visc(i,j,2*(jq-1)+iq) > & CS%min_ice_visc * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf))) then - CS%newton_visc_factor(i,j,2*(jq-1)+iq) = (0.5*(1./n_g - 1.) / & + CS%newton_visc_factor(i,j,2*(jq-1)+iq) = ((In_g - 1.) / & (((ux**2) + (vy**2)) + ((ux*vy) + 0.25*((uy+vx)**2)) + eps_min**2)) * & CS%ice_visc(i,j,2*(jq-1)+iq) endif @@ -4198,103 +4571,95 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) end subroutine calc_shelf_visc +!> Pre-compute element-level basal friction prefactors for quadrature-point evaluation. +subroutine calc_shelf_basal_prefactors(CS, ISS, G, US) + type(ice_shelf_dyn_CS), intent(inout) :: CS !< Ice shelf dynamics control structure + type(ice_shelf_state), intent(in) :: ISS !< Ice shelf state (hmask, h_shelf) + type(ocean_grid_type), intent(in) :: G !< The grid structure + type(unit_scale_type), intent(in) :: US !< Unit conversion factors -!> Update basal shear -subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) - type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure - type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe - !! the ice-shelf state - 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(G%IsdB:G%IedB,G%JsdB:G%JedB), & - intent(inout) :: u_shlf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. - real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & - intent(inout) :: v_shlf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + integer :: i, j, isd, ied, jsd, jed + real :: Hf ! Floatation thickness [Z ~> m] + real :: fN ! Effective pressure for Coulomb friction [R Z L T-2 ~> Pa] -! also this subroutine updates the nonlinear part of the basal traction + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed -! this may be subject to change later... to make it "hybrid" + do j = jsd, jed ; do i = isd, ied + CS%coef_prefactor(i,j) = G%areaT(i,j) * CS%C_basal_friction(i,j) * US%L_T_to_m_s + if (CS%CoulombFriction .and. (ISS%hmask(i,j) == 1 .or. ISS%hmask(i,j) == 3)) then + Hf = max(CS%rhow_rhoi * CS%bed_elev(i,j), 0.0) + fN = max((US%L_to_Z*(CS%density_ice * CS%g_Earth) * & + (max(ISS%h_shelf(i,j), CS%min_h_shelf) - Hf)), CS%CF_MinN) + CS%fB_elem(i,j) = CS%alpha_coulomb * & + (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%coulomb_pp_n) + else + CS%fB_elem(i,j) = 0.0 + endif + enddo ; enddo - integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq - integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec, is, js - real :: umid, vmid ! Velocities [L T-1 ~> m s-1] - real :: eps_min ! A minimal strain rate used in the Glens flow law expression [T-1 ~> s-1] - real :: unorm ! The magnitude of the velocity in mks units for use with fractional powers [m s-1] - real :: alpha ! Coulomb coefficient [nondim] - real :: Hf !"floatation thickness" for Coulomb friction [Z ~> m] - real :: fN ! Effective pressure (ice pressure - ocean pressure) for Coulomb friction [R Z L T-2 ~> Pa] - real :: fB !for Coulomb Friction [(T L-1)^CS%CF_PostPeak ~> (s m-1)^CS%CF_PostPeak] - real :: fBuq ! fB * unorm^CF_PostPeak, for Coulomb Newton correction [nondim] - real :: unorm_code2 ! Squared velocity magnitude in code units [L2 T-2 ~> m2 s-2] +end subroutine calc_shelf_basal_prefactors - isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec - iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB - isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed - iegq = G%iegB ; jegq = G%jegB - gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 - giec = G%domain%niglobal+gisc ; gjec = G%domain%njglobal+gjsc - is = iscq - 1 ; js = jscq - 1 +!> Compute area-averaged basal shear stress [R L T-1 ~> Pa s m-1] and return it in basal_tr. +!! Uses CS%u_shelf and CS%v_shelf for velocities and G%US for unit conversions. +subroutine calc_shelf_taub(CS, ISS, G, basal_tr) + type(ice_shelf_dyn_CS), intent(in) :: CS !< Ice shelf dynamics control structure + type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe + !! the ice-shelf state + type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(out) :: basal_tr !< Area-averaged basal traction [R L T-1 ~> Pa s m-1] + + integer :: i, j + real :: umid, vmid ! Cell-center velocity averages [L T-1 ~> m s-1] + real :: eps_min ! Minimal strain rate [T-1 ~> s-1] + real :: unorm ! Velocity magnitude in mks units [m s-1] + real :: alpha ! Coulomb coefficient [nondim] + real :: Hf ! Floatation thickness for Coulomb friction [Z ~> m] + real :: fN ! Effective pressure for Coulomb friction [R Z L T-2 ~> Pa] + real :: fB ! Coulomb friction factor [(T L-1)^CS%CF_PostPeak] + real :: fBuq ! fB * unorm^CF_PostPeak [nondim] + real :: unorm_code2 ! Squared velocity magnitude in code units [L2 T-2 ~> m2 s-2] + real :: basal_trac ! Area-integrated traction coefficient [R Z L2 T-1 ~> kg s-1] eps_min = CS%eps_glen_min if (CS%CoulombFriction) then - if (CS%CF_PostPeak/=1.0) THEN - alpha = (CS%CF_PostPeak-1.0)**(CS%CF_PostPeak-1.0) / CS%CF_PostPeak**CS%CF_PostPeak ![nondim] + if (CS%CF_PostPeak /= 1.0) then + alpha = CS%alpha_coulomb else alpha = 1.0 endif endif - do j=jsd+1,jed - do i=isd+1,ied - CS%newton_drag_coef(i,j) = 0.0 - if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then - umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 - vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 - unorm_code2 = ((umid**2) + (vmid**2)) + (eps_min**2 * ((G%dxT(i,j)**2) + (G%dyT(i,j)**2))) - unorm = US%L_T_to_m_s * sqrt( unorm_code2 ) - - !Coulomb friction (Schoof 2005, Gagliardini et al 2007) - if (CS%CoulombFriction) then - !Effective pressure - Hf = max((CS%density_ocean_avg/CS%density_ice) * CS%bed_elev(i,j), 0.0) - fN = max((US%L_to_Z*(CS%density_ice * CS%g_Earth) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) - Hf)), CS%CF_MinN) - fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric) - fBuq = fB * unorm**CS%CF_PostPeak - - CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * & - (unorm**(CS%n_basal_fric-1.0) / (1.0 + fBuq)**(CS%n_basal_fric))) * & - US%L_T_to_m_s ! Restore the scaling after the fractional power law. - else - !linear (CS%n_basal_fric=1) or "Weertman"/power-law (CS%n_basal_fric /= 1) - fBuq = 0.0 - CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * (unorm**(CS%n_basal_fric-1))) * & - US%L_T_to_m_s ! Rescale after the fractional power law. - endif + basal_tr(:,:) = 0.0 - CS%basal_traction(i,j)=max(CS%basal_traction(i,j), CS%min_basal_traction * G%areaT(i,j)) - - ! Store Newton basal drag data for Newton tangent stiffness correction. - ! newton_drag_coef = 2 * d(basal_trac)/d(|u|^2), - ! where d(tau_b_i)/d(u_j) = basal_trac*delta_ij + newton_drag_coef*u_i*u_j - ! This is the coefficient of the rank-1 correction u_i*(u.delta_u) to the Picard basal stiffness. - ! For Weertman: newton_drag_coef = (m-1) * basal_trac/|u|^2 - ! For Coulomb: newton_drag_coef = basal_trac/|u|^2 * [(m-1) - m*q*fB*|u|^q/(1+fB*|u|^q)] - CS%newton_umid(i,j) = umid - CS%newton_vmid(i,j) = vmid - ! unorm_code2: squared velocity magnitude in code units [L2 T-2], including regularization - ! (same expression as inside the sqrt in unorm, without US%L_T_to_m_s factor) - if (CS%CoulombFriction) then - CS%newton_drag_coef(i,j) = (1.0 / max(unorm_code2, epsilon(unorm_code2))) * & - CS%basal_traction(i,j) * ((CS%n_basal_fric - 1.) - & - CS%n_basal_fric * CS%CF_PostPeak * fBuq / (1. + fBuq)) - else - CS%newton_drag_coef(i,j) = real(CS%n_basal_fric - 1.) * CS%basal_traction(i,j) / & - max(unorm_code2, epsilon(unorm_code2)) - endif + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then + umid = ((CS%u_shelf(I,J) + CS%u_shelf(I-1,J-1)) + (CS%u_shelf(I,J-1) + CS%u_shelf(I-1,J))) * 0.25 + vmid = ((CS%v_shelf(I,J) + CS%v_shelf(I-1,J-1)) + (CS%v_shelf(I,J-1) + CS%v_shelf(I-1,J))) * 0.25 + unorm_code2 = ((umid**2) + (vmid**2)) + (eps_min**2 * ((G%dxT(i,j)**2) + (G%dyT(i,j)**2))) + unorm = G%US%L_T_to_m_s * sqrt(unorm_code2) + + !Coulomb friction (Schoof 2005, Gagliardini et al 2007) + if (CS%CoulombFriction) then + !Effective pressure + Hf = max(CS%rhow_rhoi * CS%bed_elev(i,j), 0.0) + fN = max((G%US%L_to_Z*(CS%density_ice * CS%g_Earth) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) - Hf)), CS%CF_MinN) + fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%coulomb_pp_n) + fBuq = fB * unorm**CS%CF_PostPeak + basal_trac = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * & + (unorm**(CS%n_basal_fric-1.0) / (1.0 + fBuq)**(CS%n_basal_fric))) * & + G%US%L_T_to_m_s ! Restore the scaling after the fractional power law. + else + !linear (CS%n_basal_fric = 1) or "Weertman"/power-law (CS%n_basal_fric /= 1) + basal_trac = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * (unorm**(CS%n_basal_fric-1))) * & + G%US%L_T_to_m_s ! Rescale after the fractional power law. endif - enddo - enddo + + basal_trac = max(basal_trac, CS%min_basal_traction * G%areaT(i,j)) + basal_tr(i,j) = basal_trac * G%IareaT(i,j) * CS%ground_frac(i,j) + endif + enddo ; enddo end subroutine calc_shelf_taub @@ -4345,14 +4710,13 @@ subroutine update_OD_ffrac_uncoupled(CS, G, h_shelf) intent(in) :: h_shelf !< the thickness of the ice shelf [Z ~> m]. integer :: i, j, isd, ied, jsd, jed - real :: rhoi_rhow, OD + real :: OD - rhoi_rhow = CS%density_ice / CS%density_ocean_avg isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed do j=jsd,jed do i=isd,ied - OD = CS%bed_elev(i,j) - rhoi_rhow * max(h_shelf(i,j),CS%min_h_shelf) + OD = CS%bed_elev(i,j) - CS%rhoi_rhow * max(h_shelf(i,j),CS%min_h_shelf) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -4377,9 +4741,8 @@ subroutine change_in_draft(CS, G, h_shelf0, h_shelf1, ddraft) intent(inout) :: ddraft !< the change in shelf draft thickness real :: b0,b1 integer :: i, j, isc, iec, jsc, jec - real :: rhoi_rhow, OD + real :: OD - rhoi_rhow = CS%density_ice / CS%density_ocean_avg isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ddraft = 0.0 @@ -4389,20 +4752,20 @@ subroutine change_in_draft(CS, G, h_shelf0, h_shelf1, ddraft) b0 = 0.0 ; b1 = 0.0 if (h_shelf0(i,j)>0.0) then - OD = CS%bed_elev(i,j) - rhoi_rhow * h_shelf0(i,j) + OD = CS%bed_elev(i,j) - CS%rhoi_rhow * h_shelf0(i,j) if (OD >= 0) then !floating - b0 = rhoi_rhow * h_shelf0(i,j) + b0 = CS%rhoi_rhow * h_shelf0(i,j) else b0 = CS%bed_elev(i,j) endif endif if (h_shelf1(i,j)>0.0) then - OD = CS%bed_elev(i,j) - rhoi_rhow * h_shelf1(i,j) + OD = CS%bed_elev(i,j) - CS%rhoi_rhow * h_shelf1(i,j) if (OD >= 0) then !floating - b1 = rhoi_rhow * h_shelf1(i,j) + b1 = CS%rhoi_rhow * h_shelf1(i,j) else b1 = CS%bed_elev(i,j) endif @@ -4877,8 +5240,8 @@ subroutine ice_shelf_dyn_end(CS) deallocate(CS%ice_visc, CS%AGlen_visc) deallocate(CS%newton_visc_factor, CS%newton_str_ux, CS%newton_str_vy, CS%newton_str_sh) - deallocate(CS%newton_umid, CS%newton_vmid, CS%newton_drag_coef) - deallocate(CS%basal_traction,CS%C_basal_friction) + deallocate(CS%C_basal_friction) + deallocate(CS%coef_prefactor, CS%fB_elem) deallocate(CS%OD_rt, CS%OD_av) deallocate(CS%t_bdry_val, CS%bed_elev) deallocate(CS%ground_frac, CS%ground_frac_rt)