MOM6
Data Types | Functions/Subroutines
mom_ice_shelf_dynamics Module Reference

Detailed Description

Implements a crude placeholder for a later implementation of full ice shelf dynamics.

Data Types

type  ice_shelf_dyn_cs
 The control structure for the ice shelf dynamics. More...
 
type  loop_bounds_type
 A container for loop bounds. More...
 

Functions/Subroutines

real function slope_limiter (num, denom)
 used for flux limiting in advective subroutines Van Leer limiter (source: Wikipedia) The return value is between 0 and 2 [nondim]. More...
 
real function quad_area (X, Y)
 Calculate area of quadrilateral. More...
 
subroutine, public register_ice_shelf_dyn_restarts (G, param_file, CS, restart_CS)
 This subroutine is used to register any fields related to the ice shelf dynamics that should be written to or read from the restart file. More...
 
subroutine, public initialize_ice_shelf_dyn (param_file, Time, ISS, CS, G, US, diag, new_sim, solo_ice_sheet_in)
 Initializes shelf model data, parameters and diagnostics. More...
 
subroutine initialize_diagnostic_fields (CS, ISS, G, US, Time)
 
real function, public ice_time_step_cfl (CS, ISS, G)
 This function returns the global maximum advective timestep that can be taken based on the current ice velocities. Because it involves finding a global minimum, it can be surprisingly expensive. More...
 
subroutine, public update_ice_shelf (CS, ISS, G, US, time_step, Time, ocean_mass, coupled_grounding, must_update_vel)
 This subroutine updates the ice shelf velocities, mass, stresses and properties due to the ice shelf dynamics. More...
 
subroutine ice_shelf_advect (CS, ISS, G, time_step, Time)
 This subroutine takes the velocity (on the Bgrid) and timesteps h_t = - div (uh) once. Additionally, it will update the volume of ice in partially-filled cells, and update hmask accordingly. More...
 
subroutine ice_shelf_solve_outer (CS, ISS, G, US, u_shlf, v_shlf, iters, time)
 
subroutine ice_shelf_solve_inner (CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H_node, float_cond, hmask, conv_flag, iters, time, Phi, Phisub)
 
subroutine ice_shelf_advect_thickness_x (CS, G, LB, time_step, hmask, h0, h_after_uflux, uh_ice)
 
subroutine ice_shelf_advect_thickness_y (CS, G, LB, time_step, hmask, h0, h_after_vflux, vh_ice)
 
subroutine, public shelf_advance_front (CS, ISS, G, hmask, uh_ice, vh_ice)
 
subroutine, public ice_shelf_min_thickness_calve (G, h_shelf, area_shelf_h, hmask, thickness_calve, halo)
 Apply a very simple calving law using a minimum thickness rule. More...
 
subroutine, public calve_to_mask (G, h_shelf, area_shelf_h, hmask, calve_mask)
 
subroutine calc_shelf_driving_stress (CS, ISS, G, US, taudx, taudy, OD)
 
subroutine init_boundary_values (CS, G, time, hmask, input_flux, input_thick, new_sim)
 
subroutine cg_action (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)
 
subroutine cg_action_subgrid_basal (Phisub, H, U, V, bathyT, dens_ratio, Ucontr, Vcontr)
 
subroutine matrix_diagonal (CS, G, US, float_cond, H_node, ice_visc, basal_trac, hmask, dens_ratio, Phisub, u_diagonal, v_diagonal)
 returns the diagonal entries of the matrix for a Jacobi preconditioning More...
 
subroutine cg_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, sub_grnd)
 
subroutine apply_boundary_values (CS, ISS, G, US, time, Phisub, H_node, ice_visc, basal_trac, float_cond, dens_ratio, u_bdry_contr, v_bdry_contr)
 
subroutine calc_shelf_visc (CS, ISS, G, US, u_shlf, v_shlf)
 Update depth integrated viscosity, based on horizontal strain rates, and also update the nonlinear part of the basal traction. More...
 
subroutine update_od_ffrac (CS, G, US, ocean_mass, find_avg)
 
subroutine update_od_ffrac_uncoupled (CS, G, h_shelf)
 
subroutine bilinear_shape_functions (X, Y, Phi, area)
 This subroutine calculates the gradients of bilinear basis elements that that are centered at the vertices of the cell. Values are calculated at points of gaussian quadrature. More...
 
subroutine bilinear_shape_fn_grid (G, i, j, Phi)
 This subroutine calculates the gradients of bilinear basis elements that are centered at the vertices of the cell using a locally orthogoal MOM6 grid. Values are calculated at points of gaussian quadrature. More...
 
subroutine bilinear_shape_functions_subgrid (Phisub, nsub)
 
subroutine update_velocity_masks (CS, G, hmask, umask, vmask, u_face_mask, v_face_mask)
 
subroutine interpolate_h_to_b (G, h_shelf, hmask, H_node)
 Interpolate the ice shelf thickness from tracer point to nodal points, subject to a mask. More...
 
subroutine, public ice_shelf_dyn_end (CS)
 Deallocates all memory associated with the ice shelf dynamics module. More...
 
subroutine ice_shelf_temp (CS, ISS, G, US, time_step, melt_rate, Time)
 This subroutine updates the vertically averaged ice shelf temperature. More...
 
subroutine ice_shelf_advect_temp_x (CS, G, time_step, hmask, h0, h_after_uflux)
 
subroutine ice_shelf_advect_temp_y (CS, G, time_step, hmask, h_after_uflux, h_after_vflux)
 

Function/Subroutine Documentation

◆ apply_boundary_values()

subroutine mom_ice_shelf_dynamics::apply_boundary_values ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(time_type), intent(in)  time,
real, dimension(:,:,:,:,:,:), intent(in)  Phisub,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  H_node,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  ice_visc,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  basal_trac,
real, dimension(szdi_(g),szdj_(g)), intent(in)  float_cond,
real, intent(in)  dens_ratio,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  u_bdry_contr,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  v_bdry_contr 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in]issA structure with elements that describe the ice-shelf state
[in]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]timeThe current model time
[in]phisubQuadrature structure weights at subgridscale
[in]h_nodeThe ice shelf thickness at nodal
[in]ice_viscA field related to the ice viscosity from Glen's
[in]basal_tracA field related to the nonlinear part of the
[in]float_condAn array indicating where the ice
[in]dens_ratioThe density of ice divided by the density of seawater, nondimensional
[in,out]u_bdry_contrZonal force contributions due to the
[in,out]v_bdry_contrMeridional force contributions due to the

Definition at line 2308 of file MOM_ice_shelf_dynamics.F90.

2310 
2311  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
2312  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
2313  !! the ice-shelf state
2314  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2315  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
2316  type(time_type), intent(in) :: Time !< The current model time
2317  real, dimension(:,:,:,:,:,:), &
2318  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2319  !! locations for finite element calculations [nondim]
2320  real, dimension(SZDIB_(G),SZDJB_(G)), &
2321  intent(in) :: H_node !< The ice shelf thickness at nodal
2322  !! (corner) points [Z ~> m].
2323  real, dimension(SZDIB_(G),SZDJB_(G)), &
2324  intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's
2325  !! flow law. The exact form and units depend on the
2326  !! basal law exponent. [R L4 Z T-1 ~> kg m2 s-1].
2327  real, dimension(SZDIB_(G),SZDJB_(G)), &
2328  intent(in) :: basal_trac !< A field related to the nonlinear part of the
2329  !! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1].
2330  real, dimension(SZDI_(G),SZDJ_(G)), &
2331  intent(in) :: float_cond !< An array indicating where the ice
2332  !! shelf is floating: 0 if floating, 1 if not.
2333  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2334  !! of seawater, nondimensional
2335  real, dimension(SZDIB_(G),SZDJB_(G)), &
2336  intent(inout) :: u_bdry_contr !< Zonal force contributions due to the
2337  !! open boundaries [R L3 Z T-2 ~> kg m s-2]
2338  real, dimension(SZDIB_(G),SZDJB_(G)), &
2339  intent(inout) :: v_bdry_contr !< Meridional force contributions due to the
2340  !! open boundaries [R L3 Z T-2 ~> kg m s-2]
2341 
2342 ! this will be a per-setup function. the boundary values of thickness and velocity
2343 ! (and possibly other variables) will be updated in this function
2344 
2345  real, dimension(8,4) :: Phi
2346  real, dimension(2) :: xquad
2347  real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1]
2348  real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1]
2349  real :: area
2350  real, dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr
2351  integer :: i, j, isc, jsc, iec, jec, iq, jq, iphi, jphi, ilq, jlq, Itgt, Jtgt
2352 
2353  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2354 
2355  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2356 
2357  do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (iss%hmask(i,j) == 1) then
2358 
2359  ! process this cell if any corners have umask set to non-dirichlet bdry.
2360  ! NOTE: vmask not considered, probably should be
2361 
2362  if ((cs%umask(i-1,j-1) == 3) .OR. (cs%umask(i,j-1) == 3) .OR. &
2363  (cs%umask(i-1,j) == 3) .OR. (cs%umask(i,j) == 3)) then
2364 
2365  call bilinear_shape_fn_grid(g, i, j, phi)
2366 
2367  ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2368  ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j
2369 
2370  do iq=1,2 ; do jq=1,2
2371 
2372  uq = cs%u_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2373  cs%u_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
2374  cs%u_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
2375  cs%u_bdry_val(i,j) * xquad(iq) * xquad(jq)
2376 
2377  vq = cs%v_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2378  cs%v_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
2379  cs%v_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
2380  cs%v_bdry_val(i,j) * xquad(iq) * xquad(jq)
2381 
2382  ux = cs%u_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
2383  cs%u_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
2384  cs%u_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
2385  cs%u_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
2386 
2387  vx = cs%v_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
2388  cs%v_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
2389  cs%v_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
2390  cs%v_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
2391 
2392  uy = cs%u_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
2393  cs%u_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
2394  cs%u_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
2395  cs%u_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
2396 
2397  vy = cs%v_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
2398  cs%v_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
2399  cs%v_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
2400  cs%v_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
2401 
2402  do iphi=1,2 ; do jphi=1,2 ; itgt = i-2+iphi ; jtgt = j-2-jphi
2403  ilq = 1 ; if (iq == iphi) ilq = 2
2404  jlq = 1 ; if (jq == jphi) jlq = 2
2405 
2406  if (cs%umask(itgt,jtgt) == 1) then
2407  u_bdry_contr(itgt,jtgt) = u_bdry_contr(itgt,jtgt) + &
2408  0.25 * ice_visc(i,j) * ( (4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2409  (uy+vx) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) )
2410 
2411  if (float_cond(i,j) == 0) then
2412  u_bdry_contr(itgt,jtgt) = u_bdry_contr(itgt,jtgt) + &
2413  0.25 * basal_trac(i,j) * uq * xquad(ilq) * xquad(jlq)
2414  endif
2415  endif
2416 
2417  if (cs%vmask(itgt,jtgt) == 1) then
2418  v_bdry_contr(itgt,jtgt) = v_bdry_contr(itgt,jtgt) + &
2419  0.25 * ice_visc(i,j) * ( (uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2420  (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) )
2421 
2422  if (float_cond(i,j) == 0) then
2423  v_bdry_contr(itgt,jtgt) = v_bdry_contr(itgt,jtgt) + &
2424  0.25 * basal_trac(i,j) * vq * xquad(ilq) * xquad(jlq)
2425  endif
2426  endif
2427  enddo ; enddo
2428  enddo ; enddo
2429 
2430  if (float_cond(i,j) == 1) then
2431  ucell(:,:) = cs%u_bdry_val(i-1:i,j-1:j) ; vcell(:,:) = cs%v_bdry_val(i-1:i,j-1:j)
2432  hcell(:,:) = h_node(i-1:i,j-1:j)
2433  call cg_action_subgrid_basal(phisub, hcell, ucell, vcell, g%bathyT(i,j), &
2434  dens_ratio, usubcontr, vsubcontr)
2435 
2436  if (cs%umask(i-1,j-1)==1) u_bdry_contr(i-1,j-1) = u_bdry_contr(i-1,j-1) + usubcontr(1,1) * basal_trac(i,j)
2437  if (cs%umask(i-1,j) == 1) u_bdry_contr(i-1,j) = u_bdry_contr(i-1,j) + usubcontr(1,2) * basal_trac(i,j)
2438  if (cs%umask(i,j-1) == 1) u_bdry_contr(i,j-1) = u_bdry_contr(i,j-1) + usubcontr(2,1) * basal_trac(i,j)
2439  if (cs%umask(i,j) == 1) u_bdry_contr(i,j) = u_bdry_contr(i,j) + usubcontr(2,2) * basal_trac(i,j)
2440 
2441  if (cs%vmask(i-1,j-1)==1) v_bdry_contr(i-1,j-1) = v_bdry_contr(i-1,j-1) + vsubcontr(1,1) * basal_trac(i,j)
2442  if (cs%vmask(i-1,j) == 1) v_bdry_contr(i-1,j) = v_bdry_contr(i-1,j) + vsubcontr(1,2) * basal_trac(i,j)
2443  if (cs%vmask(i,j-1) == 1) v_bdry_contr(i,j-1) = v_bdry_contr(i,j-1) + vsubcontr(2,1) * basal_trac(i,j)
2444  if (cs%vmask(i,j) == 1) v_bdry_contr(i,j) = v_bdry_contr(i,j) + vsubcontr(2,2) * basal_trac(i,j)
2445  endif
2446  endif
2447  endif ; enddo ; enddo
2448 

◆ bilinear_shape_fn_grid()

subroutine mom_ice_shelf_dynamics::bilinear_shape_fn_grid ( type(ocean_grid_type), intent(in)  G,
integer, intent(in)  i,
integer, intent(in)  j,
real, dimension(8,4), intent(inout)  Phi 
)
private

This subroutine calculates the gradients of bilinear basis elements that are centered at the vertices of the cell using a locally orthogoal MOM6 grid. Values are calculated at points of gaussian quadrature.

Parameters
[in]gThe grid structure used by the ice shelf.
[in]iThe i-index in the grid to work on.
[in]jThe j-index in the grid to work on.
[in,out]phiThe gradients of bilinear basis elements at Gaussian quadrature points surrounding the cell vertices [L-1 ~> m-1].

Definition at line 2649 of file MOM_ice_shelf_dynamics.F90.

2650  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2651  integer, intent(in) :: i !< The i-index in the grid to work on.
2652  integer, intent(in) :: j !< The j-index in the grid to work on.
2653  real, dimension(8,4), intent(inout) :: Phi !< The gradients of bilinear basis elements at Gaussian
2654  !! quadrature points surrounding the cell vertices [L-1 ~> m-1].
2655 
2656 ! This subroutine calculates the gradients of bilinear basis elements that
2657 ! that are centered at the vertices of the cell. The values are calculated at
2658 ! points of gaussian quadrature. (in 1D: .5 * (1 +/- sqrt(1/3)) for [0,1])
2659 ! (ordered in same way as vertices)
2660 !
2661 ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2662 ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j
2663 ! Phi_i is equal to 1 at vertex i, and 0 at vertex k /= i, and bilinear
2664 !
2665 ! This should be a one-off; once per nonlinear solve? once per lifetime?
2666 
2667  real, dimension(4) :: xquad, yquad ! [nondim]
2668  real :: a, d ! Interpolated grid spacings [L ~> m]
2669  real :: xexp, yexp ! [nondim]
2670  integer :: node, qpoint, xnode, xq, ynode, yq
2671 
2672  xquad(1:3:2) = .5 * (1-sqrt(1./3)) ; yquad(1:2) = .5 * (1-sqrt(1./3))
2673  xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3))
2674 
2675  do qpoint=1,4
2676  a = g%dxCv(i,j-1) * (1-yquad(qpoint)) + g%dxCv(i,j) * yquad(qpoint) ! d(x)/d(x*)
2677  d = g%dyCu(i-1,j) * (1-xquad(qpoint)) + g%dyCu(i,j) * xquad(qpoint) ! d(y)/d(y*)
2678 
2679  do node=1,4
2680  xnode = 2-mod(node,2) ; ynode = ceiling(real(node)/2)
2681 
2682  if (ynode == 1) then
2683  yexp = 1-yquad(qpoint)
2684  else
2685  yexp = yquad(qpoint)
2686  endif
2687 
2688  if (1 == xnode) then
2689  xexp = 1-xquad(qpoint)
2690  else
2691  xexp = xquad(qpoint)
2692  endif
2693 
2694  phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp ) / (a*d)
2695  phi(2*node,qpoint) = ( a * (2 * ynode - 3) * xexp ) / (a*d)
2696 
2697  enddo
2698  enddo
2699 

◆ bilinear_shape_functions()

subroutine mom_ice_shelf_dynamics::bilinear_shape_functions ( real, dimension(4), intent(in)  X,
real, dimension(4), intent(in)  Y,
real, dimension(8,4), intent(inout)  Phi,
real, intent(out)  area 
)
private

This subroutine calculates the gradients of bilinear basis elements that that are centered at the vertices of the cell. Values are calculated at points of gaussian quadrature.

Parameters
[in]xThe x-positions of the vertices of the quadrilateral [L ~> m].
[in]yThe y-positions of the vertices of the quadrilateral [L ~> m].
[in,out]phiThe gradients of bilinear basis elements at Gaussian quadrature points surrounding the cell vertices [L-1 ~> m-1].
[out]areaThe quadrilateral cell area [L2 ~> m2].

Definition at line 2581 of file MOM_ice_shelf_dynamics.F90.

2582  real, dimension(4), intent(in) :: X !< The x-positions of the vertices of the quadrilateral [L ~> m].
2583  real, dimension(4), intent(in) :: Y !< The y-positions of the vertices of the quadrilateral [L ~> m].
2584  real, dimension(8,4), intent(inout) :: Phi !< The gradients of bilinear basis elements at Gaussian
2585  !! quadrature points surrounding the cell vertices [L-1 ~> m-1].
2586  real, intent(out) :: area !< The quadrilateral cell area [L2 ~> m2].
2587 
2588 ! X and Y must be passed in the form
2589  ! 3 - 4
2590  ! | |
2591  ! 1 - 2
2592 
2593 ! this subroutine calculates the gradients of bilinear basis elements that
2594 ! that are centered at the vertices of the cell. values are calculated at
2595 ! points of gaussian quadrature. (in 1D: .5 * (1 +/- sqrt(1/3)) for [0,1])
2596 ! (ordered in same way as vertices)
2597 !
2598 ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2599 ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j
2600 ! Phi_i is equal to 1 at vertex i, and 0 at vertex k /= i, and bilinear
2601 !
2602 ! This should be a one-off; once per nonlinear solve? once per lifetime?
2603 ! ... will all cells have the same shape and dimension?
2604 
2605  real, dimension(4) :: xquad, yquad ! [nondim]
2606  real :: a,b,c,d ! Various lengths [L ~> m]
2607  real :: xexp, yexp ! [nondim]
2608  integer :: node, qpoint, xnode, xq, ynode, yq
2609 
2610  xquad(1:3:2) = .5 * (1-sqrt(1./3)) ; yquad(1:2) = .5 * (1-sqrt(1./3))
2611  xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3))
2612 
2613  do qpoint=1,4
2614 
2615  a = -x(1)*(1-yquad(qpoint)) + x(2)*(1-yquad(qpoint)) - x(3)*yquad(qpoint) + x(4)*yquad(qpoint) ! d(x)/d(x*)
2616  b = -y(1)*(1-yquad(qpoint)) + y(2)*(1-yquad(qpoint)) - y(3)*yquad(qpoint) + y(4)*yquad(qpoint) ! d(y)/d(x*)
2617  c = -x(1)*(1-xquad(qpoint)) - x(2)*(xquad(qpoint)) + x(3)*(1-xquad(qpoint)) + x(4)*(xquad(qpoint)) ! d(x)/d(y*)
2618  d = -y(1)*(1-xquad(qpoint)) - y(2)*(xquad(qpoint)) + y(3)*(1-xquad(qpoint)) + y(4)*(xquad(qpoint)) ! d(y)/d(y*)
2619 
2620  do node=1,4
2621 
2622  xnode = 2-mod(node,2) ; ynode = ceiling(real(node)/2)
2623 
2624  if (ynode == 1) then
2625  yexp = 1-yquad(qpoint)
2626  else
2627  yexp = yquad(qpoint)
2628  endif
2629 
2630  if (1 == xnode) then
2631  xexp = 1-xquad(qpoint)
2632  else
2633  xexp = xquad(qpoint)
2634  endif
2635 
2636  phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / (a*d-b*c)
2637  phi(2*node,qpoint) = (-c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / (a*d-b*c)
2638 
2639  enddo
2640  enddo
2641 
2642  area = quad_area(x, y)
2643 

◆ bilinear_shape_functions_subgrid()

subroutine mom_ice_shelf_dynamics::bilinear_shape_functions_subgrid ( real, dimension(nsub,nsub,2,2,2,2), intent(inout)  Phisub,
integer, intent(in)  nsub 
)
private
Parameters
[in,out]phisubQuadrature structure weights at subgridscale
[in]nsubThe number of subgridscale quadrature locations in each direction

Definition at line 2703 of file MOM_ice_shelf_dynamics.F90.

2704  real, dimension(nsub,nsub,2,2,2,2), &
2705  intent(inout) :: Phisub !< Quadrature structure weights at subgridscale
2706  !! locations for finite element calculations [nondim]
2707  integer, intent(in) :: nsub !< The number of subgridscale quadrature locations in each direction
2708 
2709  ! this subroutine is a helper for interpolation of floatation condition
2710  ! for the purposes of evaluating the terms \int (u,v) \phi_i dx dy in a cell that is
2711  ! in partial floatation
2712  ! the array Phisub contains the values of \phi_i (where i is a node of the cell)
2713  ! at quad point j
2714  ! i think this general approach may not work for nonrectangular elements...
2715  !
2716 
2717  ! Phisub(i,j,k,l,q1,q2)
2718  ! i: subgrid index in x-direction
2719  ! j: subgrid index in y-direction
2720  ! k: basis function x-index
2721  ! l: basis function y-index
2722  ! q1: quad point x-index
2723  ! q2: quad point y-index
2724 
2725  ! e.g. k=1,l=1 => node 1
2726  ! q1=2,q2=1 => quad point 2
2727 
2728  ! 3 - 4
2729  ! | |
2730  ! 1 - 2
2731 
2732  integer :: i, j, k, l, qx, qy, indx, indy
2733  real,dimension(2) :: xquad
2734  real :: x0, y0, x, y, val, fracx
2735 
2736  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2737  fracx = 1.0/real(nsub)
2738 
2739  do j=1,nsub ; do i=1,nsub
2740  x0 = (i-1) * fracx ; y0 = (j-1) * fracx
2741  do qy=1,2 ; do qx=1,2
2742  x = x0 + fracx*xquad(qx)
2743  y = y0 + fracx*xquad(qy)
2744  phisub(i,j,1,1,qx,qy) = (1.0-x) * (1.0-y)
2745  phisub(i,j,1,2,qx,qy) = (1.0-x) * y
2746  phisub(i,j,2,1,qx,qy) = x * (1.0-y)
2747  phisub(i,j,2,2,qx,qy) = x * y
2748  enddo ; enddo
2749  enddo ; enddo
2750 

◆ calc_shelf_driving_stress()

subroutine mom_ice_shelf_dynamics::calc_shelf_driving_stress ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  taudx,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  taudy,
real, dimension(szdi_(g),szdj_(g)), intent(in)  OD 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]odocean floor depth at tracer points [Z ~> m].
[in,out]taudxX-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
[in,out]taudyY-direction driving stress at q-points [kg L s-2 ~> kg m s-2]

Definition at line 1706 of file MOM_ice_shelf_dynamics.F90.

1707  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1708  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
1709  !! the ice-shelf state
1710  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
1711  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
1712  real, dimension(SZDI_(G),SZDJ_(G)), &
1713  intent(in) :: OD !< ocean floor depth at tracer points [Z ~> m].
1714  real, dimension(SZDIB_(G),SZDJB_(G)), &
1715  intent(inout) :: taudx !< X-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
1716  real, dimension(SZDIB_(G),SZDJB_(G)), &
1717  intent(inout) :: taudy !< Y-direction driving stress at q-points [kg L s-2 ~> kg m s-2]
1718  ! This will become [R L3 Z T-2 ~> kg m s-2]
1719 
1720 ! driving stress!
1721 
1722 ! ! taudx and taudy will hold driving stress in the x- and y- directions when done.
1723 ! they will sit on the BGrid, and so their size depends on whether the grid is symmetric
1724 !
1725 ! Since this is a finite element solve, they will actually have the form \int \Phi_i rho g h \nabla s
1726 !
1727 ! OD -this is important and we do not yet know where (in MOM) it will come from. It represents
1728 ! "average" ocean depth -- and is needed to find surface elevation
1729 ! (it is assumed that base_ice = bed + OD)
1730 
1731  real, dimension(SIZE(OD,1),SIZE(OD,2)) :: S, & ! surface elevation [Z ~> m].
1732  BASE ! basal elevation of shelf/stream [Z ~> m].
1733 
1734 
1735  real :: rho, rhow ! Ice and ocean densities [R ~> kg m-3]
1736  real :: sx, sy ! Ice shelf top slopes [Z L-1 ~> m s-1]
1737  real :: neumann_val ! [R Z L2 T-2 ~> kg s-2]
1738  real :: dxh, dyh ! Local grid spacing [L ~> m]
1739  real :: grav ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
1740 
1741  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
1742  integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
1743  integer :: i_off, j_off
1744 
1745  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
1746  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
1747  isd = g%isd ; jsd = g%jsd
1748  iegq = g%iegB ; jegq = g%jegB
1749  gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
1750  giec = g%domain%niglobal+g%domain%nihalo ; gjec = g%domain%njglobal+g%domain%njhalo
1751  is = iscq - 1; js = jscq - 1
1752  i_off = g%idg_offset ; j_off = g%jdg_offset
1753 
1754  rho = cs%density_ice
1755  rhow = cs%density_ocean_avg
1756  grav = cs%g_Earth
1757 
1758  ! prelim - go through and calculate S
1759 
1760  ! or is this faster?
1761  base(:,:) = -g%bathyT(:,:) + od(:,:)
1762  s(:,:) = base(:,:) + iss%h_shelf(:,:)
1763 
1764  do j=jsc-1,jec+1
1765  do i=isc-1,iec+1
1766  cnt = 0
1767  sx = 0
1768  sy = 0
1769  dxh = g%dxT(i,j)
1770  dyh = g%dyT(i,j)
1771 
1772  if (iss%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell
1773 
1774  ! calculate sx
1775  if ((i+i_off) == gisc) then ! at left computational bdry
1776  if (iss%hmask(i+1,j) == 1) then
1777  sx = (s(i+1,j)-s(i,j))/dxh
1778  else
1779  sx = 0
1780  endif
1781  elseif ((i+i_off) == giec) then ! at east computational bdry
1782  if (iss%hmask(i-1,j) == 1) then
1783  sx = (s(i,j)-s(i-1,j))/dxh
1784  else
1785  sx = 0
1786  endif
1787  else ! interior
1788  if (iss%hmask(i+1,j) == 1) then
1789  cnt = cnt+1
1790  sx = s(i+1,j)
1791  else
1792  sx = s(i,j)
1793  endif
1794  if (iss%hmask(i-1,j) == 1) then
1795  cnt = cnt+1
1796  sx = sx - s(i-1,j)
1797  else
1798  sx = sx - s(i,j)
1799  endif
1800  if (cnt == 0) then
1801  sx = 0
1802  else
1803  sx = sx / (cnt * dxh)
1804  endif
1805  endif
1806 
1807  cnt = 0
1808 
1809  ! calculate sy, similarly
1810  if ((j+j_off) == gjsc) then ! at south computational bdry
1811  if (iss%hmask(i,j+1) == 1) then
1812  sy = (s(i,j+1)-s(i,j))/dyh
1813  else
1814  sy = 0
1815  endif
1816  elseif ((j+j_off) == gjec) then ! at nprth computational bdry
1817  if (iss%hmask(i,j-1) == 1) then
1818  sy = (s(i,j)-s(i,j-1))/dyh
1819  else
1820  sy = 0
1821  endif
1822  else ! interior
1823  if (iss%hmask(i,j+1) == 1) then
1824  cnt = cnt+1
1825  sy = s(i,j+1)
1826  else
1827  sy = s(i,j)
1828  endif
1829  if (iss%hmask(i,j-1) == 1) then
1830  cnt = cnt+1
1831  sy = sy - s(i,j-1)
1832  else
1833  sy = sy - s(i,j)
1834  endif
1835  if (cnt == 0) then
1836  sy = 0
1837  else
1838  sy = sy / (cnt * dyh)
1839  endif
1840  endif
1841 
1842  ! SW vertex
1843  taudx(i-1,j-1) = taudx(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * g%areaT(i,j)
1844  taudy(i-1,j-1) = taudy(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * g%areaT(i,j)
1845 
1846  ! SE vertex
1847  taudx(i,j-1) = taudx(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * g%areaT(i,j)
1848  taudy(i,j-1) = taudy(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * g%areaT(i,j)
1849 
1850  ! NW vertex
1851  taudx(i-1,j) = taudx(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * g%areaT(i,j)
1852  taudy(i-1,j) = taudy(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * g%areaT(i,j)
1853 
1854  ! NE vertex
1855  taudx(i,j) = taudx(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * g%areaT(i,j)
1856  taudy(i,j) = taudy(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * g%areaT(i,j)
1857 
1858  if (cs%ground_frac(i,j) == 1) then
1859  neumann_val = .5 * grav * (rho * iss%h_shelf(i,j)**2 - rhow * g%bathyT(i,j)**2)
1860  else
1861  neumann_val = .5 * grav * (1-rho/rhow) * rho * iss%h_shelf(i,j)**2
1862  endif
1863 
1864  if ((cs%u_face_mask(i-1,j) == 2) .OR. (iss%hmask(i-1,j) == 0) .OR. (iss%hmask(i-1,j) == 2) ) then
1865  ! left face of the cell is at a stress boundary
1866  ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated
1867  ! pressure on either side of the face
1868  ! on the ice side, it is rho g h^2 / 2
1869  ! on the ocean side, it is rhow g (delta OD)^2 / 2
1870  ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation
1871  ! is not above the base of the ice in the current cell
1872 
1873  ! Note the negative sign due to the direction of the normal vector
1874  taudx(i-1,j-1) = taudx(i-1,j-1) - .5 * dyh * neumann_val
1875  taudx(i-1,j) = taudx(i-1,j) - .5 * dyh * neumann_val
1876  endif
1877 
1878  if ((cs%u_face_mask(i,j) == 2) .OR. (iss%hmask(i+1,j) == 0) .OR. (iss%hmask(i+1,j) == 2) ) then
1879  ! east face of the cell is at a stress boundary
1880  taudx(i,j-1) = taudx(i,j-1) + .5 * dyh * neumann_val
1881  taudx(i,j) = taudx(i,j) + .5 * dyh * neumann_val
1882  endif
1883 
1884  if ((cs%v_face_mask(i,j-1) == 2) .OR. (iss%hmask(i,j-1) == 0) .OR. (iss%hmask(i,j-1) == 2) ) then
1885  ! south face of the cell is at a stress boundary
1886  taudy(i-1,j-1) = taudy(i-1,j-1) - .5 * dxh * neumann_val
1887  taudy(i,j-1) = taudy(i,j-1) - .5 * dxh * neumann_val
1888  endif
1889 
1890  if ((cs%v_face_mask(i,j) == 2) .OR. (iss%hmask(i,j+1) == 0) .OR. (iss%hmask(i,j+1) == 2) ) then
1891  ! north face of the cell is at a stress boundary
1892  taudy(i-1,j) = taudy(i-1,j) + .5 * dxh * neumann_val
1893  taudy(i,j) = taudy(i,j) + .5 * dxh * neumann_val
1894  endif
1895 
1896  endif
1897  enddo
1898  enddo
1899 

◆ calc_shelf_visc()

subroutine mom_ice_shelf_dynamics::calc_shelf_visc ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(inout)  u_shlf,
real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(inout)  v_shlf 
)
private

Update depth integrated viscosity, based on horizontal strain rates, and also update the nonlinear part of the basal traction.

Parameters
[in,out]csA pointer to the ice shelf control structure
[in]issA structure with elements that describe the ice-shelf state
[in]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in,out]u_shlfThe zonal ice shelf velocity [L T-1 ~> m s-1].
[in,out]v_shlfThe meridional ice shelf velocity [L T-1 ~> m s-1].

Definition at line 2453 of file MOM_ice_shelf_dynamics.F90.

2454  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
2455  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
2456  !! the ice-shelf state
2457  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2458  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
2459  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2460  intent(inout) :: u_shlf !< The zonal ice shelf velocity [L T-1 ~> m s-1].
2461  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2462  intent(inout) :: v_shlf !< The meridional ice shelf velocity [L T-1 ~> m s-1].
2463 
2464 ! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve
2465 ! so there is an "upper" and "lower" bilinear viscosity
2466 
2467 ! also this subroutine updates the nonlinear part of the basal traction
2468 
2469 ! this may be subject to change later... to make it "hybrid"
2470 
2471  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
2472  integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js
2473  real :: Visc_coef, n_g
2474  real :: ux, uy, vx, vy, eps_min ! Velocity shears [T-1 ~> s-1]
2475  real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1]
2476 
2477  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2478  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
2479  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
2480  iegq = g%iegB ; jegq = g%jegB
2481  gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
2482  giec = g%domain%niglobal+gisc ; gjec = g%domain%njglobal+gjsc
2483  is = iscq - 1; js = jscq - 1
2484 
2485  n_g = cs%n_glen; eps_min = cs%eps_glen_min
2486 
2487  visc_coef = us%kg_m2s_to_RZ_T*us%m_to_L*us%Z_to_L*(cs%A_glen_isothermal)**(1./cs%n_glen)
2488 
2489  do j=jsd+1,jed-1
2490  do i=isd+1,ied-1
2491 
2492  if (iss%hmask(i,j) == 1) then
2493  ux = ((u_shlf(i,j) + u_shlf(i,j-1)) - (u_shlf(i-1,j) + u_shlf(i-1,j-1))) / (2*g%dxT(i,j))
2494  vx = ((v_shlf(i,j) + v_shlf(i,j-1)) - (v_shlf(i-1,j) + v_shlf(i-1,j-1))) / (2*g%dxT(i,j))
2495  uy = ((u_shlf(i,j) + u_shlf(i-1,j)) - (u_shlf(i,j-1) + u_shlf(i-1,j-1))) / (2*g%dyT(i,j))
2496  vy = ((v_shlf(i,j) + v_shlf(i-1,j)) - (v_shlf(i,j-1) + v_shlf(i-1,j-1))) / (2*g%dyT(i,j))
2497  cs%ice_visc(i,j) = 0.5 * visc_coef * (g%areaT(i,j) * iss%h_shelf(i,j)) * &
2498  (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))
2499 
2500  umid = ((u_shlf(i,j) + u_shlf(i-1,j-1)) + (u_shlf(i,j-1) + u_shlf(i-1,j))) * 0.25
2501  vmid = ((v_shlf(i,j) + v_shlf(i-1,j-1)) + (v_shlf(i,j-1) + v_shlf(i-1,j))) * 0.25
2502  unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(g%dxT(i,j)**2 + g%dyT(i,j)**2))
2503  cs%basal_traction(i,j) = g%areaT(i,j) * cs%C_basal_friction * (us%L_T_to_m_s*unorm)**(cs%n_basal_fric-1)
2504  endif
2505  enddo
2506  enddo
2507 

◆ calve_to_mask()

subroutine, public mom_ice_shelf_dynamics::calve_to_mask ( type(ocean_grid_type), intent(in)  G,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  h_shelf,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  area_shelf_h,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  hmask,
real, dimension(szdi_(g),szdj_(g)), intent(in)  calve_mask 
)
Parameters
[in]gThe grid structure used by the ice shelf.
[in,out]h_shelfThe ice shelf thickness [Z ~> m].
[in,out]area_shelf_hThe area per cell covered by the ice shelf [L2 ~> m2].
[in,out]hmaskA mask indicating which tracer points are partly or fully covered by an ice-shelf
[in]calve_maskA mask that indicates where the ice shelf can exist, and where it will calve.

Definition at line 1684 of file MOM_ice_shelf_dynamics.F90.

1685  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1686  real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
1687  real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: area_shelf_h !< The area per cell covered by
1688  !! the ice shelf [L2 ~> m2].
1689  real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: hmask !< A mask indicating which tracer points are
1690  !! partly or fully covered by an ice-shelf
1691  real, dimension(SZDI_(G),SZDJ_(G)), intent(in) :: calve_mask !< A mask that indicates where the ice
1692  !! shelf can exist, and where it will calve.
1693 
1694  integer :: i,j
1695 
1696  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1697  if ((calve_mask(i,j) == 0.0) .and. (hmask(i,j) /= 0.0)) then
1698  h_shelf(i,j) = 0.0
1699  area_shelf_h(i,j) = 0.0
1700  hmask(i,j) = 0.0
1701  endif
1702  enddo ; enddo
1703 

◆ cg_action()

subroutine mom_ice_shelf_dynamics::cg_action ( real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(inout)  uret,
real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(inout)  vret,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  u_shlf,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  v_shlf,
real, dimension(szdi_(g),szdj_(g),8,4), intent(in)  Phi,
real, dimension(:,:,:,:,:,:), intent(in)  Phisub,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  umask,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  vmask,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  H_node,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  ice_visc,
real, dimension(szdi_(g),szdj_(g)), intent(in)  float_cond,
real, dimension(szdi_(g),szdj_(g)), intent(in)  bathyT,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  basal_trac,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  is,
integer, intent(in)  ie,
integer, intent(in)  js,
integer, intent(in)  je,
real, intent(in)  dens_ratio 
)
private
Parameters
[in]gThe grid structure used by the ice shelf.
[in,out]uretThe retarding stresses working at u-points [R L3 Z T-2 ~> kg m s-2].
[in,out]vretThe retarding stresses working at v-points [R L3 Z T-2 ~> kg m s-2].
[in]phiThe gradients of bilinear basis elements at Gaussian
[in]phisubQuadrature structure weights at subgridscale
[in]u_shlfThe zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
[in]v_shlfThe meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
[in]umaskA coded mask indicating the nature of the
[in]vmaskA coded mask indicating the nature of the
[in]h_nodeThe ice shelf thickness at nodal (corner)
[in]hmaskA mask indicating which tracer points are
[in]ice_viscA field related to the ice viscosity from Glen's
[in]float_condAn array indicating where the ice
[in]bathytThe depth of ocean bathymetry at tracer points [Z ~> m].
[in]basal_tracA field related to the nonlinear part of the
[in]dens_ratioThe density of ice divided by the density of seawater, nondimensional
[in]usA structure containing unit conversion factors
[in]isThe starting i-index to work on
[in]ieThe ending i-index to work on
[in]jsThe starting j-index to work on
[in]jeThe ending j-index to work on

Definition at line 1971 of file MOM_ice_shelf_dynamics.F90.

1973 
1974  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1975  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
1976  intent(inout) :: uret !< The retarding stresses working at u-points [R L3 Z T-2 ~> kg m s-2].
1977  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
1978  intent(inout) :: vret !< The retarding stresses working at v-points [R L3 Z T-2 ~> kg m s-2].
1979  real, dimension(SZDI_(G),SZDJ_(G),8,4), &
1980  intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian
1981  !! quadrature points surrounding the cell vertices [L-1 ~> m-1].
1982  real, dimension(:,:,:,:,:,:), &
1983  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
1984  !! locations for finite element calculations [nondim]
1985  real, dimension(SZDIB_(G),SZDJB_(G)), &
1986  intent(in) :: u_shlf !< The zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
1987  real, dimension(SZDIB_(G),SZDJB_(G)), &
1988  intent(in) :: v_shlf !< The meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
1989  real, dimension(SZDIB_(G),SZDJB_(G)), &
1990  intent(in) :: umask !< A coded mask indicating the nature of the
1991  !! zonal flow at the corner point
1992  real, dimension(SZDIB_(G),SZDJB_(G)), &
1993  intent(in) :: vmask !< A coded mask indicating the nature of the
1994  !! meridional flow at the corner point
1995  real, dimension(SZDIB_(G),SZDJB_(G)), &
1996  intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
1997  !! points [Z ~> m].
1998  real, dimension(SZDI_(G),SZDJ_(G)), &
1999  intent(in) :: hmask !< A mask indicating which tracer points are
2000  !! partly or fully covered by an ice-shelf
2001  real, dimension(SZDIB_(G),SZDJB_(G)), &
2002  intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's
2003  !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form
2004  !! and units depend on the basal law exponent.
2005  real, dimension(SZDI_(G),SZDJ_(G)), &
2006  intent(in) :: float_cond !< An array indicating where the ice
2007  !! shelf is floating: 0 if floating, 1 if not.
2008  real, dimension(SZDI_(G),SZDJ_(G)), &
2009  intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
2010  real, dimension(SZDIB_(G),SZDJB_(G)), &
2011  intent(in) :: basal_trac !< A field related to the nonlinear part of the
2012  !! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1].
2013  ! and/or whether flow is "hybridized"
2014  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2015  !! of seawater, nondimensional
2016  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
2017  integer, intent(in) :: is !< The starting i-index to work on
2018  integer, intent(in) :: ie !< The ending i-index to work on
2019  integer, intent(in) :: js !< The starting j-index to work on
2020  integer, intent(in) :: je !< The ending j-index to work on
2021 
2022 ! the linear action of the matrix on (u,v) with bilinear finite elements
2023 ! as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
2024 ! but this may change pursuant to conversations with others
2025 !
2026 ! is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
2027 ! in order to make less frequent halo updates
2028 
2029 ! the linear action of the matrix on (u,v) with bilinear finite elements
2030 ! Phi has the form
2031 ! Phi(k,q,i,j) - applies to cell i,j
2032 
2033  ! 3 - 4
2034  ! | |
2035  ! 1 - 2
2036 
2037 ! Phi(2*k-1,q,i,j) gives d(Phi_k)/dx at quadrature point q
2038 ! Phi(2*k,q,i,j) gives d(Phi_k)/dy at quadrature point q
2039 ! Phi_k is equal to 1 at vertex k, and 0 at vertex l /= k, and bilinear
2040 
2041  real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1]
2042  real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1]
2043  integer :: iq, jq, iphi, jphi, i, j, ilq, jlq, Itgt, Jtgt
2044  real, dimension(2) :: xquad
2045  real, dimension(2,2) :: Ucell, Vcell, Hcell, Usub, Vsub
2046 
2047  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2048 
2049  do j=js,je ; do i=is,ie ; if (hmask(i,j) == 1) then
2050 
2051  do iq=1,2 ; do jq=1,2
2052 
2053  uq = u_shlf(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2054  u_shlf(i,j-1) * xquad(iq) * xquad(3-jq) + &
2055  u_shlf(i-1,j) * xquad(3-iq) * xquad(jq) + &
2056  u_shlf(i,j) * xquad(iq) * xquad(jq)
2057 
2058  vq = v_shlf(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2059  v_shlf(i,j-1) * xquad(iq) * xquad(3-jq) + &
2060  v_shlf(i-1,j) * xquad(3-iq) * xquad(jq) + &
2061  v_shlf(i,j) * xquad(iq) * xquad(jq)
2062 
2063  ux = u_shlf(i-1,j-1) * phi(1,2*(jq-1)+iq,i,j) + &
2064  u_shlf(i,j-1) * phi(3,2*(jq-1)+iq,i,j) + &
2065  u_shlf(i-1,j) * phi(5,2*(jq-1)+iq,i,j) + &
2066  u_shlf(i,j) * phi(7,2*(jq-1)+iq,i,j)
2067 
2068  vx = v_shlf(i-1,j-1) * phi(1,2*(jq-1)+iq,i,j) + &
2069  v_shlf(i,j-1) * phi(3,2*(jq-1)+iq,i,j) + &
2070  v_shlf(i-1,j) * phi(5,2*(jq-1)+iq,i,j) + &
2071  v_shlf(i,j) * phi(7,2*(jq-1)+iq,i,j)
2072 
2073  uy = u_shlf(i-1,j-1) * phi(2,2*(jq-1)+iq,i,j) + &
2074  u_shlf(i,j-1) * phi(4,2*(jq-1)+iq,i,j) + &
2075  u_shlf(i-1,j) * phi(6,2*(jq-1)+iq,i,j) + &
2076  u_shlf(i,j) * phi(8,2*(jq-1)+iq,i,j)
2077 
2078  vy = v_shlf(i-1,j-1) * phi(2,2*(jq-1)+iq,i,j) + &
2079  v_shlf(i,j-1) * phi(4,2*(jq-1)+iq,i,j) + &
2080  v_shlf(i-1,j) * phi(6,2*(jq-1)+iq,i,j) + &
2081  v_shlf(i,j) * phi(8,2*(jq-1)+iq,i,j)
2082 
2083  do iphi=1,2 ; do jphi=1,2 ; itgt = i-2+iphi ; jtgt = j-2-jphi
2084  if (umask(itgt,jtgt) == 1) uret(itgt,jtgt) = uret(itgt,jtgt) + 0.25 * ice_visc(i,j) * &
2085  ((4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + &
2086  (uy+vx) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j))
2087  if (vmask(itgt,jtgt) == 1) vret(itgt,jtgt) = vret(itgt,jtgt) + 0.25 * ice_visc(i,j) * &
2088  ((uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + &
2089  (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j))
2090 
2091  if (float_cond(i,j) == 0) then
2092  ilq = 1 ; if (iq == iphi) ilq = 2
2093  jlq = 1 ; if (jq == jphi) jlq = 2
2094  if (umask(itgt,jtgt) == 1) uret(itgt,jtgt) = uret(itgt,jtgt) + &
2095  0.25 * basal_trac(i,j) * uq * xquad(ilq) * xquad(jlq)
2096  if (vmask(itgt,jtgt) == 1) vret(itgt,jtgt) = vret(itgt,jtgt) + &
2097  0.25 * basal_trac(i,j) * vq * xquad(ilq) * xquad(jlq)
2098  endif
2099  enddo ; enddo
2100  enddo ; enddo
2101 
2102  if (float_cond(i,j) == 1) then
2103  ucell(:,:) = u_shlf(i-1:i,j-1:j) ; vcell(:,:) = v_shlf(i-1:i,j-1:j)
2104  hcell(:,:) = h_node(i-1:i,j-1:j)
2105  call cg_action_subgrid_basal(phisub, hcell, ucell, vcell, bathyt(i,j), dens_ratio, usub, vsub)
2106 
2107  if (umask(i-1,j-1)==1) uret(i-1,j-1) = uret(i-1,j-1) + usub(1,1) * basal_trac(i,j)
2108  if (umask(i-1,j) == 1) uret(i-1,j) = uret(i-1,j) + usub(1,2) * basal_trac(i,j)
2109  if (umask(i,j-1) == 1) uret(i,j-1) = uret(i,j-1) + usub(2,1) * basal_trac(i,j)
2110  if (umask(i,j) == 1) uret(i,j) = uret(i,j) + usub(2,2) * basal_trac(i,j)
2111 
2112  if (vmask(i-1,j-1)==1) vret(i-1,j-1) = vret(i-1,j-1) + vsub(1,1) * basal_trac(i,j)
2113  if (vmask(i-1,j) == 1) vret(i-1,j) = vret(i-1,j) + vsub(1,2) * basal_trac(i,j)
2114  if (vmask(i,j-1) == 1) vret(i,j-1) = vret(i,j-1) + vsub(2,1) * basal_trac(i,j)
2115  if (vmask(i,j) == 1) vret(i,j) = vret(i,j) + vsub(2,2) * basal_trac(i,j)
2116  endif
2117 
2118  endif ; enddo ; enddo
2119 

◆ cg_action_subgrid_basal()

subroutine mom_ice_shelf_dynamics::cg_action_subgrid_basal ( real, dimension(:,:,:,:,:,:), intent(in)  Phisub,
real, dimension(2,2), intent(in)  H,
real, dimension(2,2), intent(in)  U,
real, dimension(2,2), intent(in)  V,
real, intent(in)  bathyT,
real, intent(in)  dens_ratio,
real, dimension(2,2), intent(out)  Ucontr,
real, dimension(2,2), intent(out)  Vcontr 
)
private
Parameters
[in]phisubQuadrature structure weights at subgridscale
[in]hThe ice shelf thickness at nodal (corner) points [Z ~> m].
[in]uThe zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
[in]vThe meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
[in]bathytThe depth of ocean bathymetry at tracer points [Z ~> m].
[in]dens_ratioThe density of ice divided by the density of seawater [nondim]
[out]ucontrThe areal average of u-velocities where the ice shelf is grounded, or 0 where it is floating [L T-1 ~> m s-1].
[out]vcontrThe areal average of v-velocities where the ice shelf is grounded, or 0 where it is floating [L T-1 ~> m s-1].

Definition at line 2122 of file MOM_ice_shelf_dynamics.F90.

2123  real, dimension(:,:,:,:,:,:), &
2124  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2125  !! locations for finite element calculations [nondim]
2126  real, dimension(2,2), intent(in) :: H !< The ice shelf thickness at nodal (corner) points [Z ~> m].
2127  real, dimension(2,2), intent(in) :: U !< The zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
2128  real, dimension(2,2), intent(in) :: V !< The meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
2129  real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
2130  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2131  !! of seawater [nondim]
2132  real, dimension(2,2), intent(out) :: Ucontr !< The areal average of u-velocities where the ice shelf
2133  !! is grounded, or 0 where it is floating [L T-1 ~> m s-1].
2134  real, dimension(2,2), intent(out) :: Vcontr !< The areal average of v-velocities where the ice shelf
2135  !! is grounded, or 0 where it is floating [L T-1 ~> m s-1].
2136 
2137  real :: subarea ! The fractional sub-cell area [nondim]
2138  real :: hloc ! The local sub-cell ice thickness [Z ~> m]
2139  integer :: nsub, i, j, qx, qy, m, n
2140 
2141  nsub = size(phisub,1)
2142  subarea = 1.0 / (nsub**2)
2143 
2144  do n=1,2 ; do m=1,2
2145  ucontr(m,n) = 0.0 ; vcontr(m,n) = 0.0
2146  do qy=1,2 ; do qx=1,2 ; do j=1,nsub ; do i=1,nsub
2147  hloc = (phisub(i,j,1,1,qx,qy)*h(1,1) + phisub(i,j,2,2,qx,qy)*h(2,2)) + &
2148  (phisub(i,j,1,2,qx,qy)*h(1,2) + phisub(i,j,2,1,qx,qy)*h(2,1))
2149  if (dens_ratio * hloc - bathyt > 0) then
2150  ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * &
2151  ((phisub(i,j,1,1,qx,qy) * u(1,1) + phisub(i,j,2,2,qx,qy) * u(2,2)) + &
2152  (phisub(i,j,1,2,qx,qy) * u(1,2) + phisub(i,j,2,1,qx,qy) * u(2,1)))
2153  vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * &
2154  ((phisub(i,j,1,1,qx,qy) * v(1,1) + phisub(i,j,2,2,qx,qy) * v(2,2)) + &
2155  (phisub(i,j,1,2,qx,qy) * v(1,2) + phisub(i,j,2,1,qx,qy) * v(2,1)))
2156  endif
2157  enddo ; enddo ; enddo ; enddo
2158  enddo ; enddo
2159 

◆ cg_diagonal_subgrid_basal()

subroutine mom_ice_shelf_dynamics::cg_diagonal_subgrid_basal ( real, dimension(:,:,:,:,:,:), intent(in)  Phisub,
real, dimension(2,2), intent(in)  H_node,
real, intent(in)  bathyT,
real, intent(in)  dens_ratio,
real, dimension(2,2), intent(out)  sub_grnd 
)
private
Parameters
[in]phisubQuadrature structure weights at subgridscale
[in]h_nodeThe ice shelf thickness at nodal (corner) points [Z ~> m].
[in]bathytThe depth of ocean bathymetry at tracer points [Z ~> m].
[in]dens_ratioThe density of ice divided by the density of seawater [nondim]
[out]sub_grndThe weighted fraction of the sub-cell where the ice shelf is grounded [nondim]

Definition at line 2272 of file MOM_ice_shelf_dynamics.F90.

2273  real, dimension(:,:,:,:,:,:), &
2274  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2275  !! locations for finite element calculations [nondim]
2276  real, dimension(2,2), intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
2277  !! points [Z ~> m].
2278  real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
2279  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2280  !! of seawater [nondim]
2281  real, dimension(2,2), intent(out) :: sub_grnd !< The weighted fraction of the sub-cell where the ice shelf
2282  !! is grounded [nondim]
2283 
2284  ! bathyT = cellwise-constant bed elevation
2285 
2286  real :: subarea ! The fractional sub-cell area [nondim]
2287  real :: hloc ! The local sub-region thickness [Z ~> m]
2288  integer :: nsub, i, j, k, l, qx, qy, m, n
2289 
2290  nsub = size(phisub,1)
2291  subarea = 1.0 / (nsub**2)
2292 
2293  sub_grnd(:,:) = 0.0
2294  do m=1,2 ; do n=1,2 ; do j=1,nsub ; do i=1,nsub ; do qx=1,2 ; do qy = 1,2
2295 
2296  hloc = (phisub(i,j,1,1,qx,qy)*h_node(1,1) + phisub(i,j,2,2,qx,qy)*h_node(2,2)) + &
2297  (phisub(i,j,1,2,qx,qy)*h_node(1,2) + phisub(i,j,2,1,qx,qy)*h_node(2,1))
2298 
2299  if (dens_ratio * hloc - bathyt > 0) then
2300  sub_grnd(m,n) = sub_grnd(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
2301  endif
2302 
2303  enddo ; enddo ; enddo ; enddo ; enddo ; enddo
2304 

◆ ice_shelf_advect()

subroutine mom_ice_shelf_dynamics::ice_shelf_advect ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(inout)  ISS,
type(ocean_grid_type), intent(inout)  G,
real, intent(in)  time_step,
type(time_type), intent(in)  Time 
)
private

This subroutine takes the velocity (on the Bgrid) and timesteps h_t = - div (uh) once. Additionally, it will update the volume of ice in partially-filled cells, and update hmask accordingly.

Parameters
[in,out]csThe ice shelf dynamics control structure
[in,out]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]time_steptime step [T ~> s]
[in]timeThe current model time

Definition at line 694 of file MOM_ice_shelf_dynamics.F90.

695  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
696  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
697  !! the ice-shelf state
698  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
699  real, intent(in) :: time_step !< time step [T ~> s]
700  type(time_type), intent(in) :: Time !< The current model time
701 
702 
703 ! 3/8/11 DNG
704 !
705 ! This subroutine takes the velocity (on the Bgrid) and timesteps h_t = - div (uh) once.
706 ! ADDITIONALLY, it will update the volume of ice in partially-filled cells, and update
707 ! hmask accordingly
708 !
709 ! The flux overflows are included here. That is because they will be used to advect 3D scalars
710 ! into partial cells
711 
712  real, dimension(SZDI_(G),SZDJ_(G)) :: h_after_uflux, h_after_vflux ! Ice thicknesses [Z ~> m].
713  real, dimension(SZDIB_(G),SZDJ_(G)) :: uh_ice ! The accumulated zonal ice volume flux [Z L2 ~> m3]
714  real, dimension(SZDI_(G),SZDJB_(G)) :: vh_ice ! The accumulated meridional ice volume flux [Z L2 ~> m3]
715  type(loop_bounds_type) :: LB
716  integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec, stencil
717 
718  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
719  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
720 
721  uh_ice(:,:) = 0.0
722  vh_ice(:,:) = 0.0
723 
724  h_after_uflux(:,:) = 0.0
725  h_after_vflux(:,:) = 0.0
726  ! call MOM_mesg("MOM_ice_shelf.F90: ice_shelf_advect called")
727 
728  do j=jsd,jed ; do i=isd,ied ; if (cs%thickness_bdry_val(i,j) /= 0.0) then
729  iss%h_shelf(i,j) = cs%thickness_bdry_val(i,j)
730  endif ; enddo ; enddo
731 
732  stencil = 2
733  lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc-stencil ; lb%jeh = g%jec+stencil
734  if (lb%jsh < jsd) call mom_error(fatal, &
735  "ice_shelf_advect: Halo is too small for the ice thickness advection stencil.")
736 
737  call ice_shelf_advect_thickness_x(cs, g, lb, time_step, iss%hmask, iss%h_shelf, h_after_uflux, uh_ice)
738 
739 ! call enable_averages(time_step, Time, CS%diag)
740 ! call pass_var(h_after_uflux, G%domain)
741 ! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag)
742 ! call disable_averaging(CS%diag)
743 
744  lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
745  call ice_shelf_advect_thickness_y(cs, g, lb, time_step, iss%hmask, h_after_uflux, h_after_vflux, vh_ice)
746 
747 ! call enable_averages(time_step, Time, CS%diag)
748 ! call pass_var(h_after_vflux, G%domain)
749 ! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag)
750 ! call disable_averaging(CS%diag)
751 
752  do j=jsd,jed
753  do i=isd,ied
754  if (iss%hmask(i,j) == 1) iss%h_shelf(i,j) = h_after_vflux(i,j)
755  enddo
756  enddo
757 
758  if (cs%moving_shelf_front) then
759  call shelf_advance_front(cs, iss, g, iss%hmask, uh_ice, vh_ice)
760  if (cs%min_thickness_simple_calve > 0.0) then
761  call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
762  cs%min_thickness_simple_calve)
763  endif
764  if (cs%calve_to_mask) then
765  call calve_to_mask(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, cs%calve_mask)
766  endif
767  endif
768 
769  !call enable_averages(time_step, Time, CS%diag)
770  !if (CS%id_h_after_adv > 0) call post_data(CS%id_h_after_adv, ISS%h_shelf, CS%diag)
771  !call disable_averaging(CS%diag)
772 
773  !call change_thickness_using_melt(ISS, G, time_step, fluxes, CS%density_ice)
774 
775  call update_velocity_masks(cs, g, iss%hmask, cs%umask, cs%vmask, cs%u_face_mask, cs%v_face_mask)
776 

◆ ice_shelf_advect_temp_x()

subroutine mom_ice_shelf_dynamics::ice_shelf_advect_temp_x ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(inout)  G,
real, intent(in)  time_step,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, dimension(szdi_(g),szdj_(g)), intent(in)  h0,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  h_after_uflux 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in,out]gThe grid structure used by the ice shelf.
[in]time_stepThe time step for this update [T ~> s].
[in]hmaskA mask indicating which tracer points are
[in]h0The initial ice shelf thicknesses [Z ~> m].
[in,out]h_after_ufluxThe ice shelf thicknesses after

Definition at line 3069 of file MOM_ice_shelf_dynamics.F90.

3070  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
3071  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
3072  real, intent(in) :: time_step !< The time step for this update [T ~> s].
3073  real, dimension(SZDI_(G),SZDJ_(G)), &
3074  intent(in) :: hmask !< A mask indicating which tracer points are
3075  !! partly or fully covered by an ice-shelf
3076  real, dimension(SZDI_(G),SZDJ_(G)), &
3077  intent(in) :: h0 !< The initial ice shelf thicknesses [Z ~> m].
3078  real, dimension(SZDI_(G),SZDJ_(G)), &
3079  intent(inout) :: h_after_uflux !< The ice shelf thicknesses after
3080  !! the zonal mass fluxes [Z ~> m].
3081 
3082  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
3083  ! if there is an input bdry condition, the thickness there will be set in initialization
3084 
3085  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3086  integer :: i_off, j_off
3087  logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
3088  real, dimension(-2:2) :: stencil
3089  real :: u_face ! Zonal velocity at a face, positive if out {L T-1 ~> m s-1]
3090  real :: flux_diff, phi
3091 
3092  character (len=1) :: debug_str
3093 
3094 
3095  is = g%isc-2 ; ie = g%iec+2 ; js = g%jsc ; je = g%jec ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3096  i_off = g%idg_offset ; j_off = g%jdg_offset
3097 
3098  do j=jsd+1,jed-1
3099  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
3100  ((j+j_off) >= g%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries
3101 
3102  stencil(:) = -1
3103 ! if (i+i_off == G%domain%nihalo+G%domain%nihalo)
3104  do i=is,ie
3105 
3106  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
3107  ((i+i_off) >= g%domain%nihalo+1)) then
3108 
3109  if (i+i_off == g%domain%nihalo+1) then
3110  at_west_bdry=.true.
3111  else
3112  at_west_bdry=.false.
3113  endif
3114 
3115  if (i+i_off == g%domain%niglobal+g%domain%nihalo) then
3116  at_east_bdry=.true.
3117  else
3118  at_east_bdry=.false.
3119  endif
3120 
3121  if (hmask(i,j) == 1) then
3122 
3123  h_after_uflux(i,j) = h0(i,j)
3124 
3125  stencil(:) = h0(i-2:i+2,j) ! fine as long has nx_halo >= 2
3126 
3127  flux_diff = 0
3128 
3129  ! 1ST DO LEFT FACE
3130 
3131  if (cs%u_face_mask(i-1,j) == 4.) then
3132 
3133  flux_diff = flux_diff + g%dyCu(i-1,j) * time_step * cs%u_flux_bdry_val(i-1,j) * &
3134  cs%t_bdry_val(i-1,j) / g%areaT(i,j)
3135  else
3136 
3137  ! get u-velocity at center of left face
3138  u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3139 
3140  if (u_face > 0) then !flux is into cell - we need info from h(i-2), h(i-1) if available
3141 
3142  ! i may not cover all the cases.. but i cover the realistic ones
3143 
3144  if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a
3145  ! thickness bdry condition, and the stencil contains it
3146  flux_diff = flux_diff + abs(u_face) * g%dyCu(i-1,j) * time_step * stencil(-1) / g%areaT(i,j)
3147 
3148  elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then ! h(i-2) and h(i-1) are valid
3149  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3150  flux_diff = flux_diff + abs(u_face) * g%dyCu(i-1,j)* time_step / g%areaT(i,j) * &
3151  (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3152 
3153  else ! h(i-1) is valid
3154  ! (o.w. flux would most likely be out of cell)
3155  ! but h(i-2) is not
3156 
3157  flux_diff = flux_diff + abs(u_face) * g%dyCu(i-1,j) * time_step / g%areaT(i,j) * stencil(-1)
3158 
3159  endif
3160 
3161  elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
3162  if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid
3163  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3164  flux_diff = flux_diff - abs(u_face) * g%dyCu(i-1,j) * time_step / g%areaT(i,j) * &
3165  (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3166 
3167  else
3168  flux_diff = flux_diff - abs(u_face) * g%dyCu(i-1,j) * time_step / g%areaT(i,j) * stencil(0)
3169  endif
3170  endif
3171  endif
3172 
3173  ! NEXT DO RIGHT FACE
3174 
3175  ! get u-velocity at center of eastern face
3176 
3177  if (cs%u_face_mask(i,j) == 4.) then
3178 
3179  flux_diff = flux_diff + g%dyCu(i,j) * time_step * cs%u_flux_bdry_val(i,j) *&
3180  cs%t_bdry_val(i+1,j) / g%areaT(i,j)
3181  else
3182 
3183  u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3184 
3185  if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available
3186 
3187  if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a
3188  ! thickness bdry condition, and the stencil contains it
3189 
3190  flux_diff = flux_diff + abs(u_face) * g%dyCu(i,j) * time_step * stencil(1) / g%areaT(i,j)
3191 
3192  elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then ! h(i+2) and h(i+1) are valid
3193 
3194  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3195  flux_diff = flux_diff + abs(u_face) * g%dyCu(i,j) * time_step / g%areaT(i,j) * &
3196  (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3197 
3198  else ! h(i+1) is valid
3199  ! (o.w. flux would most likely be out of cell)
3200  ! but h(i+2) is not
3201 
3202  flux_diff = flux_diff + abs(u_face) * g%dyCu(i,j) * time_step / g%areaT(i,j) * stencil(1)
3203 
3204  endif
3205 
3206  elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
3207 
3208  if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid
3209 
3210  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3211  flux_diff = flux_diff - abs(u_face) * g%dyCu(i,j) * time_step / g%areaT(i,j) * &
3212  (stencil(0) - phi * (stencil(0)-stencil(1))/2)
3213 
3214  else ! h(i+1) is valid (o.w. flux would most likely be out of cell) but h(i+2) is not
3215 
3216  flux_diff = flux_diff - abs(u_face) * g%dyCu(i,j) * time_step / g%areaT(i,j) * stencil(0)
3217 
3218  endif
3219 
3220  endif
3221 
3222  h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff
3223 
3224  endif
3225 
3226  endif
3227 
3228  endif
3229 
3230  enddo ! i loop
3231 
3232  endif
3233 
3234  enddo ! j loop
3235 

◆ ice_shelf_advect_temp_y()

subroutine mom_ice_shelf_dynamics::ice_shelf_advect_temp_y ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
real, intent(in)  time_step,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, dimension(szdi_(g),szdj_(g)), intent(in)  h_after_uflux,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  h_after_vflux 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in]gThe grid structure used by the ice shelf.
[in]time_stepThe time step for this update [T ~> s].
[in]hmaskA mask indicating which tracer points are
[in]h_after_ufluxThe ice shelf thicknesses after
[in,out]h_after_vfluxThe ice shelf thicknesses after

Definition at line 3238 of file MOM_ice_shelf_dynamics.F90.

3239  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
3240  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
3241  real, intent(in) :: time_step !< The time step for this update [T ~> s].
3242  real, dimension(SZDI_(G),SZDJ_(G)), &
3243  intent(in) :: hmask !< A mask indicating which tracer points are
3244  !! partly or fully covered by an ice-shelf
3245  real, dimension(SZDI_(G),SZDJ_(G)), &
3246  intent(in) :: h_after_uflux !< The ice shelf thicknesses after
3247  !! the zonal mass fluxes [Z ~> m].
3248  real, dimension(SZDI_(G),SZDJ_(G)), &
3249  intent(inout) :: h_after_vflux !< The ice shelf thicknesses after
3250  !! the meridional mass fluxes [Z ~> m].
3251 
3252  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
3253  ! if there is an input bdry condition, the thickness there will be set in initialization
3254 
3255  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3256  integer :: i_off, j_off
3257  logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
3258  real, dimension(-2:2) :: stencil
3259  real :: v_face ! Pseudo-meridional velocity at a cell face, positive if out {L T-1 ~> m s-1]
3260  real :: flux_diff, phi
3261  character(len=1) :: debug_str
3262 
3263  is = g%isc ; ie = g%iec ; js = g%jsc-1 ; je = g%jec+1 ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3264  i_off = g%idg_offset ; j_off = g%jdg_offset
3265 
3266  do i=isd+2,ied-2
3267  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
3268  ((i+i_off) >= g%domain%nihalo+1)) then ! based on mehmet's code - only if btw east & west boundaries
3269 
3270  stencil(:) = -1
3271 
3272  do j=js,je
3273 
3274  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
3275  ((j+j_off) >= g%domain%njhalo+1)) then
3276 
3277  if (j+j_off == g%domain%njhalo+1) then
3278  at_south_bdry=.true.
3279  else
3280  at_south_bdry=.false.
3281  endif
3282  if (j+j_off == g%domain%njglobal+g%domain%njhalo) then
3283  at_north_bdry=.true.
3284  else
3285  at_north_bdry=.false.
3286  endif
3287 
3288  if (hmask(i,j) == 1) then
3289  h_after_vflux(i,j) = h_after_uflux(i,j)
3290 
3291  stencil(:) = h_after_uflux(i,j-2:j+2) ! fine as long has ny_halo >= 2
3292  flux_diff = 0
3293 
3294  ! 1ST DO south FACE
3295 
3296  if (cs%v_face_mask(i,j-1) == 4.) then
3297 
3298  flux_diff = flux_diff + g%dxCv(i,j-1) * time_step * cs%v_flux_bdry_val(i,j-1) * &
3299  cs%t_bdry_val(i,j-1)/ g%areaT(i,j)
3300  else
3301 
3302  ! get u-velocity at center of west face
3303  v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
3304 
3305  if (v_face > 0) then !flux is into cell - we need info from h(j-2), h(j-1) if available
3306 
3307  ! i may not cover all the cases.. but i cover the realistic ones
3308 
3309  if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a
3310  ! thickness bdry condition, and the stencil contains it
3311  flux_diff = flux_diff + abs(v_face) * g%dxCv(i,j-1) * time_step * stencil(-1) / g%areaT(i,j)
3312 
3313  elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then ! h(j-2) and h(j-1) are valid
3314 
3315  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3316  flux_diff = flux_diff + abs(v_face) * g%dxCv(i,j-1) * time_step / g%areaT(i,j) * &
3317  (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3318 
3319  else ! h(j-1) is valid
3320  ! (o.w. flux would most likely be out of cell)
3321  ! but h(j-2) is not
3322  flux_diff = flux_diff + abs(v_face) * g%dxCv(i,j-1) * time_step / g%areaT(i,j) * stencil(-1)
3323  endif
3324 
3325  elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
3326 
3327  if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid
3328  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3329  flux_diff = flux_diff - abs(v_face) * g%dxCv(i,j-1) * time_step / g%areaT(i,j) * &
3330  (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3331  else
3332  flux_diff = flux_diff - abs(v_face) * g%dxCv(i,j-1) * time_step / g%areaT(i,j) * stencil(0)
3333  endif
3334 
3335  endif
3336 
3337  endif
3338 
3339  ! NEXT DO north FACE
3340 
3341  if (cs%v_face_mask(i,j) == 4.) then
3342  flux_diff = flux_diff + g%dxCv(i,j) * time_step * cs%v_flux_bdry_val(i,j) *&
3343  cs%t_bdry_val(i,j+1)/ g%areaT(i,j)
3344  else
3345 
3346  ! get u-velocity at center of east face
3347  v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
3348 
3349  if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available
3350 
3351  if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a
3352  ! thickness bdry condition, and the stencil contains it
3353  flux_diff = flux_diff + abs(v_face) * g%dxCv(i,j) * time_step * stencil(1) / g%areaT(i,j)
3354  elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then ! h(j+2) and h(j+1) are valid
3355  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3356  flux_diff = flux_diff + abs(v_face) * g%dxCv(i,j) * time_step / g%areaT(i,j) * &
3357  (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3358  else ! h(j+1) is valid
3359  ! (o.w. flux would most likely be out of cell)
3360  ! but h(j+2) is not
3361  flux_diff = flux_diff + abs(v_face) * g%dxCv(i,j) * time_step / g%areaT(i,j) * stencil(1)
3362  endif
3363 
3364  elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
3365 
3366  if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid
3367  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3368  flux_diff = flux_diff - abs(v_face) * g%dxCv(i,j) * time_step / g%areaT(i,j) * &
3369  (stencil(0) - phi * (stencil(0)-stencil(1))/2)
3370  else ! h(j+1) is valid
3371  ! (o.w. flux would most likely be out of cell)
3372  ! but h(j+2) is not
3373  flux_diff = flux_diff - abs(v_face) * g%dxCv(i,j) * time_step / g%areaT(i,j) * stencil(0)
3374  endif
3375 
3376  endif
3377 
3378  endif
3379 
3380  h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff
3381  endif
3382  endif
3383  enddo ! j loop
3384  endif
3385  enddo ! i loop
3386 

◆ ice_shelf_advect_thickness_x()

subroutine mom_ice_shelf_dynamics::ice_shelf_advect_thickness_x ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB,
real, intent(in)  time_step,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  hmask,
real, dimension(szdi_(g),szdj_(g)), intent(in)  h0,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  h_after_uflux,
real, dimension(szdib_(g),szdj_(g)), intent(inout)  uh_ice 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in]gThe grid structure used by the ice shelf.
[in]lbLoop bounds structure.
[in]time_stepThe time step for this update [T ~> s].
[in,out]hmaskA mask indicating which tracer points are
[in]h0The initial ice shelf thicknesses [Z ~> m].
[in,out]h_after_ufluxThe ice shelf thicknesses after
[in,out]uh_iceThe accumulated zonal ice volume flux [Z L2 ~> m3]

Definition at line 1299 of file MOM_ice_shelf_dynamics.F90.

1300  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1301  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1302  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1303  real, intent(in) :: time_step !< The time step for this update [T ~> s].
1304  real, dimension(SZDI_(G),SZDJ_(G)), &
1305  intent(inout) :: hmask !< A mask indicating which tracer points are
1306  !! partly or fully covered by an ice-shelf
1307  real, dimension(SZDI_(G),SZDJ_(G)), &
1308  intent(in) :: h0 !< The initial ice shelf thicknesses [Z ~> m].
1309  real, dimension(SZDI_(G),SZDJ_(G)), &
1310  intent(inout) :: h_after_uflux !< The ice shelf thicknesses after
1311  !! the zonal mass fluxes [Z ~> m].
1312  real, dimension(SZDIB_(G),SZDJ_(G)), &
1313  intent(inout) :: uh_ice !< The accumulated zonal ice volume flux [Z L2 ~> m3]
1314 
1315  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
1316  ! if there is an input bdry condition, the thickness there will be set in initialization
1317 
1318 
1319  integer :: i, j
1320  integer :: ish, ieh, jsh, jeh
1321  real :: u_face ! Zonal velocity at a face [L Z-1 ~> m s-1]
1322  real :: h_face ! Thickness at a face for transport [Z ~> m]
1323  real :: slope_lim ! The value of the slope limiter, in the range of 0 to 2 [nondim]
1324 
1325 ! is = G%isc-2 ; ie = G%iec+2 ; js = G%jsc ; je = G%jec
1326 ! isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1327 ! i_off = G%idg_offset ; j_off = G%jdg_offset
1328 
1329  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1330 
1331  ! hmask coded values: 1) fully covered; 2) partly covered - no export; 3) Specified boundary condition
1332  ! relevant u_face_mask coded values: 1) Normal interior point; 4) Specified flux BC
1333 
1334  do j=jsh,jeh ; do i=ish-1,ieh
1335  if (cs%u_face_mask(i,j) == 4.) then ! The flux itself is a specified boundary condition.
1336  uh_ice(i,j) = time_step * g%dyCu(i,j) * cs%u_flux_bdry_val(i,j)
1337  elseif ((hmask(i,j)==1) .or. (hmask(i+1,j) == 1)) then
1338  u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
1339  h_face = 0.0 ! This will apply when the source cell is iceless or not fully ice covered.
1340 
1341  if (u_face > 0) then
1342  if (hmask(i,j) == 3) then ! This is a open boundary inflow from the west
1343  h_face = cs%thickness_bdry_val(i,j)
1344  elseif (hmask(i,j) == 1) then ! There can be eastward flow through this face.
1345  if ((hmask(i-1,j) == 1) .and. (hmask(i+1,j) == 1)) then
1346  slope_lim = slope_limiter(h0(i,j)-h0(i-1,j), h0(i+1,j)-h0(i,j))
1347  ! This is a 2nd-order centered scheme with a slope limiter. We could try PPM here.
1348  h_face = h0(i,j) - slope_lim * 0.5 * (h0(i,j)-h0(i+1,j))
1349  else
1350  h_face = h0(i,j)
1351  endif
1352  endif
1353  else
1354  if (hmask(i+1,j) == 3) then ! This is a open boundary inflow from the east
1355  h_face = cs%thickness_bdry_val(i+1,j)
1356  elseif (hmask(i+1,j) == 1) then
1357  if ((hmask(i,j) == 1) .and. (hmask(i+2,j) == 1)) then
1358  slope_lim = slope_limiter(h0(i+1,j)-h0(i,j), h0(i+2,j)-h0(i+1,j))
1359  h_face = h0(i+1,j) - slope_lim * 0.5 * (h0(i+1,j)-h0(i,j))
1360  else
1361  h_face = h0(i+1,j)
1362  endif
1363  endif
1364  endif
1365 
1366  uh_ice(i,j) = time_step * g%dyCu(i,j) * u_face * h_face
1367  else
1368  uh_ice(i,j) = 0.0
1369  endif
1370  enddo ; enddo
1371 
1372  do j=jsh,jeh ; do i=ish,ieh
1373  if (hmask(i,j) /= 3) &
1374  h_after_uflux(i,j) = h0(i,j) + (uh_ice(i-1,j) - uh_ice(i,j)) * g%IareaT(i,j)
1375 
1376  ! Update the masks of cells that have gone from no ice to partial ice.
1377  if ((hmask(i,j) == 0) .and. ((uh_ice(i-1,j) > 0.0) .or. (uh_ice(i,j) < 0.0))) hmask(i,j) = 2
1378  enddo ; enddo
1379 

◆ ice_shelf_advect_thickness_y()

subroutine mom_ice_shelf_dynamics::ice_shelf_advect_thickness_y ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB,
real, intent(in)  time_step,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  hmask,
real, dimension(szdi_(g),szdj_(g)), intent(in)  h0,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  h_after_vflux,
real, dimension(szdi_(g),szdjb_(g)), intent(inout)  vh_ice 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in]gThe grid structure used by the ice shelf.
[in]lbLoop bounds structure.
[in]time_stepThe time step for this update [T ~> s].
[in,out]hmaskA mask indicating which tracer points are
[in]h0The initial ice shelf thicknesses [Z ~> m].
[in,out]h_after_vfluxThe ice shelf thicknesses after
[in,out]vh_iceThe accumulated meridional ice volume flux [Z L2 ~> m3]

Definition at line 1382 of file MOM_ice_shelf_dynamics.F90.

1383  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1384  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1385  type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure.
1386  real, intent(in) :: time_step !< The time step for this update [T ~> s].
1387  real, dimension(SZDI_(G),SZDJ_(G)), &
1388  intent(inout) :: hmask !< A mask indicating which tracer points are
1389  !! partly or fully covered by an ice-shelf
1390  real, dimension(SZDI_(G),SZDJ_(G)), &
1391  intent(in) :: h0 !< The initial ice shelf thicknesses [Z ~> m].
1392  real, dimension(SZDI_(G),SZDJ_(G)), &
1393  intent(inout) :: h_after_vflux !< The ice shelf thicknesses after
1394  !! the meridional mass fluxes [Z ~> m].
1395  real, dimension(SZDI_(G),SZDJB_(G)), &
1396  intent(inout) :: vh_ice !< The accumulated meridional ice volume flux [Z L2 ~> m3]
1397 
1398  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
1399  ! if there is an input bdry condition, the thickness there will be set in initialization
1400 
1401 
1402  integer :: i, j
1403  integer :: ish, ieh, jsh, jeh
1404  real :: v_face ! Pseudo-meridional velocity at a face [L Z-1 ~> m s-1]
1405  real :: h_face ! Thickness at a face for transport [Z ~> m]
1406  real :: slope_lim ! The value of the slope limiter, in the range of 0 to 2 [nondim]
1407 
1408  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1409 
1410  ! hmask coded values: 1) fully covered; 2) partly covered - no export; 3) Specified boundary condition
1411  ! relevant u_face_mask coded values: 1) Normal interior point; 4) Specified flux BC
1412 
1413  do j=jsh-1,jeh ; do i=ish,ieh
1414  if (cs%v_face_mask(i,j) == 4.) then ! The flux itself is a specified boundary condition.
1415  vh_ice(i,j) = time_step * g%dxCv(i,j) * cs%v_flux_bdry_val(i,j)
1416  elseif ((hmask(i,j)==1) .or. (hmask(i,j+1) == 1)) then
1417 
1418  v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
1419  h_face = 0.0 ! This will apply when the source cell is iceless or not fully ice covered.
1420 
1421  if (v_face > 0) then
1422  if (hmask(i,j) == 3) then ! This is a open boundary inflow from the south
1423  h_face = cs%thickness_bdry_val(i,j)
1424  elseif (hmask(i,j) == 1) then ! There can be northtward flow through this face.
1425  if ((hmask(i,j-1) == 1) .and. (hmask(i,j+1) == 1)) then
1426  slope_lim = slope_limiter(h0(i,j)-h0(i,j-1), h0(i,j+1)-h0(i,j))
1427  ! This is a 2nd-order centered scheme with a slope limiter. We could try PPM here.
1428  h_face = h0(i,j) - slope_lim * 0.5 * (h0(i,j)-h0(i,j+1))
1429  else
1430  h_face = h0(i,j)
1431  endif
1432  endif
1433  else
1434  if (hmask(i,j+1) == 3) then ! This is a open boundary inflow from the north
1435  h_face = cs%thickness_bdry_val(i,j+1)
1436  elseif (hmask(i,j+1) == 1) then
1437  if ((hmask(i,j) == 1) .and. (hmask(i,j+2) == 1)) then
1438  slope_lim = slope_limiter(h0(i,j+1)-h0(i,j), h0(i,j+2)-h0(i,j+1))
1439  h_face = h0(i,j+1) - slope_lim * 0.5 * (h0(i,j+1)-h0(i,j))
1440  else
1441  h_face = h0(i,j+1)
1442  endif
1443  endif
1444  endif
1445 
1446  vh_ice(i,j) = time_step * g%dxCv(i,j) * v_face * h_face
1447  else
1448  vh_ice(i,j) = 0.0
1449  endif
1450  enddo ; enddo
1451 
1452  do j=jsh,jeh ; do i=ish,ieh
1453  if (hmask(i,j) /= 3) &
1454  h_after_vflux(i,j) = h0(i,j) + (vh_ice(i,j-1) - vh_ice(i,j)) * g%IareaT(i,j)
1455 
1456  ! Update the masks of cells that have gone from no ice to partial ice.
1457  if ((hmask(i,j) == 0) .and. ((vh_ice(i,j-1) > 0.0) .or. (vh_ice(i,j) < 0.0))) hmask(i,j) = 2
1458  enddo ; enddo
1459 

◆ ice_shelf_dyn_end()

subroutine, public mom_ice_shelf_dynamics::ice_shelf_dyn_end ( type(ice_shelf_dyn_cs), pointer  CS)

Deallocates all memory associated with the ice shelf dynamics module.

Parameters
csA pointer to the ice shelf dynamics control structure

Definition at line 2954 of file MOM_ice_shelf_dynamics.F90.

2955  type(ice_shelf_dyn_CS), pointer :: CS !< A pointer to the ice shelf dynamics control structure
2956 
2957  if (.not.associated(cs)) return
2958 
2959  deallocate(cs%u_shelf, cs%v_shelf)
2960  deallocate(cs%t_shelf, cs%tmask)
2961  deallocate(cs%u_bdry_val, cs%v_bdry_val, cs%t_bdry_val)
2962  deallocate(cs%u_face_mask, cs%v_face_mask)
2963  deallocate(cs%umask, cs%vmask)
2964 
2965  deallocate(cs%ice_visc, cs%basal_traction)
2966  deallocate(cs%OD_rt, cs%OD_av)
2967  deallocate(cs%ground_frac, cs%ground_frac_rt)
2968 
2969  deallocate(cs)
2970 

◆ ice_shelf_min_thickness_calve()

subroutine, public mom_ice_shelf_dynamics::ice_shelf_min_thickness_calve ( type(ocean_grid_type), intent(in)  G,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  h_shelf,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  area_shelf_h,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  hmask,
real, intent(in)  thickness_calve,
integer, intent(in), optional  halo 
)

Apply a very simple calving law using a minimum thickness rule.

Parameters
[in]gThe grid structure used by the ice shelf.
[in,out]h_shelfThe ice shelf thickness [Z ~> m].
[in,out]area_shelf_hThe area per cell covered by the ice shelf [L2 ~> m2].
[in,out]hmaskA mask indicating which tracer points are partly or fully covered by an ice-shelf
[in]thickness_calveThe thickness at which to trigger calving [Z ~> m].
[in]haloThe number of halo points to use. If not present, work on the entire data domain.

Definition at line 1654 of file MOM_ice_shelf_dynamics.F90.

1655  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1656  real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
1657  real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: area_shelf_h !< The area per cell covered by
1658  !! the ice shelf [L2 ~> m2].
1659  real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: hmask !< A mask indicating which tracer points are
1660  !! partly or fully covered by an ice-shelf
1661  real, intent(in) :: thickness_calve !< The thickness at which to trigger calving [Z ~> m].
1662  integer, optional, intent(in) :: halo !< The number of halo points to use. If not present,
1663  !! work on the entire data domain.
1664  integer :: i, j, is, ie, js, je
1665 
1666  if (present(halo)) then
1667  is = g%isc - halo ; ie = g%iec + halo ; js = g%jsc - halo ; je = g%jec + halo
1668  else
1669  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
1670  endif
1671 
1672  do j=js,je ; do i=is,ie
1673 ! if ((h_shelf(i,j) < CS%thickness_calve) .and. (hmask(i,j) == 1) .and. &
1674 ! (CS%ground_frac(i,j) == 0.0)) then
1675  if ((h_shelf(i,j) < thickness_calve) .and. (area_shelf_h(i,j) > 0.)) then
1676  h_shelf(i,j) = 0.0
1677  area_shelf_h(i,j) = 0.0
1678  hmask(i,j) = 0.0
1679  endif
1680  enddo ; enddo
1681 

◆ ice_shelf_solve_inner()

subroutine mom_ice_shelf_dynamics::ice_shelf_solve_inner ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  u_shlf,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  v_shlf,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  taudx,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  taudy,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  H_node,
real, dimension(szdi_(g),szdj_(g)), intent(in)  float_cond,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
integer, intent(out)  conv_flag,
integer, intent(out)  iters,
type(time_type), intent(in)  time,
real, dimension(8,4,szdi_(g),szdj_(g)), intent(in)  Phi,
real, dimension(:,:,:,:,:,:), intent(in)  Phisub 
)
private
Parameters
[in]csA pointer to the ice shelf control structure
[in]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in,out]u_shlfThe zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
[in,out]v_shlfThe meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
[in]taudxThe x-direction driving stress [R L3 Z T-2 ~> kg m s-2]
[in]taudyThe y-direction driving stress [R L3 Z T-2 ~> kg m s-2]
[in]h_nodeThe ice shelf thickness at nodal (corner)
[in]float_condAn array indicating where the ice
[in]hmaskA mask indicating which tracer points are
[out]conv_flagA flag indicating whether (1) or not (0) the iterations have converged to the specified tolerance
[out]itersThe number of iterations used in the solver.
[in]timeThe current model time
[in]phiThe gradients of bilinear basis elements at Gaussian
[in]phisubQuadrature structure weights at subgridscale

Definition at line 1001 of file MOM_ice_shelf_dynamics.F90.

1003  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1004  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
1005  !! the ice-shelf state
1006  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
1007  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
1008  real, dimension(SZDIB_(G),SZDJB_(G)), &
1009  intent(inout) :: u_shlf !< The zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
1010  real, dimension(SZDIB_(G),SZDJB_(G)), &
1011  intent(inout) :: v_shlf !< The meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
1012  real, dimension(SZDIB_(G),SZDJB_(G)), &
1013  intent(in) :: taudx !< The x-direction driving stress [R L3 Z T-2 ~> kg m s-2]
1014  real, dimension(SZDIB_(G),SZDJB_(G)), &
1015  intent(in) :: taudy !< The y-direction driving stress [R L3 Z T-2 ~> kg m s-2]
1016  real, dimension(SZDIB_(G),SZDJB_(G)), &
1017  intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
1018  !! points [Z ~> m].
1019  real, dimension(SZDI_(G),SZDJ_(G)), &
1020  intent(in) :: float_cond !< An array indicating where the ice
1021  !! shelf is floating: 0 if floating, 1 if not.
1022  real, dimension(SZDI_(G),SZDJ_(G)), &
1023  intent(in) :: hmask !< A mask indicating which tracer points are
1024  !! partly or fully covered by an ice-shelf
1025  integer, intent(out) :: conv_flag !< A flag indicating whether (1) or not (0) the
1026  !! iterations have converged to the specified tolerance
1027  integer, intent(out) :: iters !< The number of iterations used in the solver.
1028  type(time_type), intent(in) :: Time !< The current model time
1029  real, dimension(8,4,SZDI_(G),SZDJ_(G)), &
1030  intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian
1031  !! quadrature points surrounding the cell vertices [L-1 ~> m-1].
1032  real, dimension(:,:,:,:,:,:), &
1033  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
1034  !! locations for finite element calculations [nondim]
1035 ! one linear solve (nonlinear iteration) of the solution for velocity
1036 
1037 ! in this subroutine:
1038 ! boundary contributions are added to taud to get the RHS
1039 ! diagonal of matrix is found (for Jacobi precondition)
1040 ! CG iteration is carried out for max. iterations or until convergence
1041 
1042 ! assumed - u, v, taud, visc, basal_traction are valid on the halo
1043 
1044  real, dimension(SZDIB_(G),SZDJB_(G)) :: &
1045  Ru, Rv, & ! Residuals in the stress calculations [R L3 Z T-2 ~> m kg s-2]
1046  Ru_old, Rv_old, & ! Previous values of Ru and Rv [R L3 Z T-2 ~> m kg s-2]
1047  Zu, Zv, & ! Contributions to velocity changes [L T-1 ~> m s-1]
1048  Zu_old, Zv_old, & ! Previous values of Zu and Zv [L T-1 ~> m s-1]
1049  DIAGu, DIAGv, & ! Diagonals with units like Ru/Zu [R L2 Z T-1 ~> kg s-1]
1050  RHSu, RHSv, & ! Right hand side of the stress balance [R L3 Z T-2 ~> m kg s-2]
1051  ubd, vbd, & ! Boundary stress contributions [R L3 Z T-2 ~> kg m s-2]
1052  Au, Av, & ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2]
1053  Du, Dv, & ! Velocity changes [L T-1 ~> m s-1]
1054  sum_vec, sum_vec_2
1055  real :: tol, beta_k, area, dot_p1, resid0, cg_halo
1056  real :: num, denom
1057  real :: alpha_k ! A scaling factor for iterative corrections [nondim]
1058  real :: resid_scale ! A scaling factor for redimensionalizing the global residuals [m2 L-2 ~> 1]
1059  ! [m2 L-2 ~> 1] [R L3 Z T-2 ~> m kg s-2]
1060  real :: resid2_scale ! A scaling factor for redimensionalizing the global squared residuals
1061  ! [m2 L-2 ~> 1] [R L3 Z T-2 ~> m kg s-2]
1062  real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim]
1063  integer :: iter, i, j, isd, ied, jsd, jed, isc, iec, jsc, jec, is, js, ie, je
1064  integer :: Is_sum, Js_sum, Ie_sum, Je_sum ! Loop bounds for global sums or arrays starting at 1.
1065  integer :: isdq, iedq, jsdq, jedq, iscq, iecq, jscq, jecq, nx_halo, ny_halo
1066 
1067  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
1068  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
1069  ny_halo = g%domain%njhalo ; nx_halo = g%domain%nihalo
1070  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1071  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1072 
1073  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
1074 
1075  zu(:,:) = 0 ; zv(:,:) = 0 ; diagu(:,:) = 0 ; diagv(:,:) = 0
1076  ru(:,:) = 0 ; rv(:,:) = 0 ; au(:,:) = 0 ; av(:,:) = 0
1077  du(:,:) = 0 ; dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0
1078  dot_p1 = 0
1079 
1080  ! Determine the loop limits for sums, bearing in mind that the arrays will be starting at 1.
1081  is_sum = g%isc + (1-g%IsdB)
1082  ie_sum = g%iecB + (1-g%IsdB)
1083  ! Include the edge if tile is at the western bdry; Should add a test to avoid this if reentrant.
1084  if (g%isc+g%idg_offset==g%isg) is_sum = g%IscB + (1-g%IsdB)
1085 
1086  js_sum = g%jsc + (1-g%JsdB)
1087  je_sum = g%jecB + (1-g%JsdB)
1088  ! Include the edge if tile is at the southern bdry; Should add a test to avoid this if reentrant.
1089  if (g%jsc+g%jdg_offset==g%jsg) js_sum = g%JscB + (1-g%JsdB)
1090 
1091  call apply_boundary_values(cs, iss, g, us, time, phisub, h_node, cs%ice_visc, &
1092  cs%basal_traction, float_cond, rhoi_rhow, ubd, vbd)
1093 
1094  rhsu(:,:) = taudx(:,:) - ubd(:,:)
1095  rhsv(:,:) = taudy(:,:) - vbd(:,:)
1096 
1097  call pass_vector(rhsu, rhsv, g%domain, to_all, bgrid_ne)
1098 
1099  call matrix_diagonal(cs, g, us, float_cond, h_node, cs%ice_visc, cs%basal_traction, &
1100  hmask, rhoi_rhow, phisub, diagu, diagv)
1101 
1102  call pass_vector(diagu, diagv, g%domain, to_all, bgrid_ne)
1103 
1104  call cg_action(au, av, u_shlf, v_shlf, phi, phisub, cs%umask, cs%vmask, hmask, &
1105  h_node, cs%ice_visc, float_cond, g%bathyT, cs%basal_traction, &
1106  g, us, isc-1, iec+1, jsc-1, jec+1, rhoi_rhow)
1107 
1108  call pass_vector(au, av, g%domain, to_all, bgrid_ne)
1109 
1110  ru(:,:) = (rhsu(:,:) - au(:,:))
1111  rv(:,:) = (rhsv(:,:) - av(:,:))
1112 
1113  resid_scale = us%L_to_m**2*us%s_to_T*us%RZ_to_kg_m2*us%L_T_to_m_s**2
1114  resid2_scale = (us%RZ_to_kg_m2*us%L_to_m*us%L_T_to_m_s**2)**2
1115 
1116  sum_vec(:,:) = 0.0
1117  do j=jscq,jecq ; do i=iscq,iecq
1118  if (cs%umask(i,j) == 1) sum_vec(i,j) = resid2_scale*ru(i,j)**2
1119  if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + resid2_scale*rv(i,j)**2
1120  enddo ; enddo
1121 
1122  dot_p1 = reproducing_sum( sum_vec, js_sum, ie_sum, js_sum, je_sum )
1123 
1124  resid0 = sqrt(dot_p1)
1125 
1126  do j=jsdq,jedq
1127  do i=isdq,iedq
1128  if (cs%umask(i,j) == 1) zu(i,j) = ru(i,j) / diagu(i,j)
1129  if (cs%vmask(i,j) == 1) zv(i,j) = rv(i,j) / diagv(i,j)
1130  enddo
1131  enddo
1132 
1133  du(:,:) = zu(:,:) ; dv(:,:) = zv(:,:)
1134 
1135  cg_halo = 3
1136  conv_flag = 0
1137 
1138  !!!!!!!!!!!!!!!!!!
1139  !! !!
1140  !! MAIN CG LOOP !!
1141  !! !!
1142  !!!!!!!!!!!!!!!!!!
1143 
1144  ! initially, c-grid data is valid up to 3 halo nodes out
1145 
1146  do iter = 1,cs%cg_max_iterations
1147 
1148  ! assume asymmetry
1149  ! thus we can never assume that any arrays are legit more than 3 vertices past
1150  ! the computational domain - this is their state in the initial iteration
1151 
1152 
1153  is = isc - cg_halo ; ie = iecq + cg_halo
1154  js = jscq - cg_halo ; je = jecq + cg_halo
1155 
1156  au(:,:) = 0 ; av(:,:) = 0
1157 
1158  call cg_action(au, av, du, dv, phi, phisub, cs%umask, cs%vmask, hmask, &
1159  h_node, cs%ice_visc, float_cond, g%bathyT, cs%basal_traction, &
1160  g, us, is, ie, js, je, rhoi_rhow)
1161 
1162  ! Au, Av valid region moves in by 1
1163 
1164 
1165  sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
1166 
1167  do j=jscq,jecq ; do i=iscq,iecq
1168  if (cs%umask(i,j) == 1) then
1169  sum_vec(i,j) = resid_scale * zu(i,j) * ru(i,j)
1170  sum_vec_2(i,j) = resid_scale * du(i,j) * au(i,j)
1171  endif
1172  if (cs%vmask(i,j) == 1) then
1173  sum_vec(i,j) = sum_vec(i,j) + resid_scale * zv(i,j) * rv(i,j)
1174  sum_vec_2(i,j) = sum_vec_2(i,j) + resid_scale * dv(i,j) * av(i,j)
1175  endif
1176  enddo ; enddo
1177 
1178  alpha_k = reproducing_sum( sum_vec, is_sum, ie_sum, js_sum, je_sum ) / &
1179  reproducing_sum( sum_vec_2, is_sum, ie_sum, js_sum, je_sum )
1180 
1181 
1182  do j=jsd,jed ; do i=isd,ied
1183  if (cs%umask(i,j) == 1) u_shlf(i,j) = u_shlf(i,j) + alpha_k * du(i,j)
1184  if (cs%vmask(i,j) == 1) v_shlf(i,j) = v_shlf(i,j) + alpha_k * dv(i,j)
1185  enddo ; enddo
1186 
1187  do j=jsd,jed ; do i=isd,ied
1188  if (cs%umask(i,j) == 1) then
1189  ru_old(i,j) = ru(i,j) ; zu_old(i,j) = zu(i,j)
1190  endif
1191  if (cs%vmask(i,j) == 1) then
1192  rv_old(i,j) = rv(i,j) ; zv_old(i,j) = zv(i,j)
1193  endif
1194  enddo ; enddo
1195 
1196 ! Ru(:,:) = Ru(:,:) - alpha_k * Au(:,:)
1197 ! Rv(:,:) = Rv(:,:) - alpha_k * Av(:,:)
1198 
1199  do j=jsd,jed
1200  do i=isd,ied
1201  if (cs%umask(i,j) == 1) ru(i,j) = ru(i,j) - alpha_k * au(i,j)
1202  if (cs%vmask(i,j) == 1) rv(i,j) = rv(i,j) - alpha_k * av(i,j)
1203  enddo
1204  enddo
1205 
1206  do j=jsdq,jedq
1207  do i=isdq,iedq
1208  if (cs%umask(i,j) == 1) then
1209  zu(i,j) = ru(i,j) / diagu(i,j)
1210  endif
1211  if (cs%vmask(i,j) == 1) then
1212  zv(i,j) = rv(i,j) / diagv(i,j)
1213  endif
1214  enddo
1215  enddo
1216 
1217  ! R,u,v,Z valid region moves in by 1
1218 
1219  ! beta_k = (Z \dot R) / (Zold \dot Rold}
1220  sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
1221 
1222  do j=jscq,jecq ; do i=iscq,iecq
1223  if (cs%umask(i,j) == 1) then
1224  sum_vec(i,j) = resid_scale * zu(i,j) * ru(i,j)
1225  sum_vec_2(i,j) = resid_scale * zu_old(i,j) * ru_old(i,j)
1226  endif
1227  if (cs%vmask(i,j) == 1) then
1228  sum_vec(i,j) = sum_vec(i,j) + resid_scale * zv(i,j) * rv(i,j)
1229  sum_vec_2(i,j) = sum_vec_2(i,j) + resid_scale * zv_old(i,j) * rv_old(i,j)
1230  endif
1231  enddo ; enddo
1232 
1233  beta_k = reproducing_sum(sum_vec, is_sum, ie_sum, js_sum, je_sum ) / &
1234  reproducing_sum(sum_vec_2, is_sum, ie_sum, js_sum, je_sum )
1235 
1236 ! Du(:,:) = Zu(:,:) + beta_k * Du(:,:)
1237 ! Dv(:,:) = Zv(:,:) + beta_k * Dv(:,:)
1238 
1239  do j=jsd,jed
1240  do i=isd,ied
1241  if (cs%umask(i,j) == 1) du(i,j) = zu(i,j) + beta_k * du(i,j)
1242  if (cs%vmask(i,j) == 1) dv(i,j) = zv(i,j) + beta_k * dv(i,j)
1243  enddo
1244  enddo
1245 
1246  ! D valid region moves in by 1
1247 
1248  sum_vec(:,:) = 0.0
1249  do j=jscq,jecq ; do i=iscq,iecq
1250  if (cs%umask(i,j) == 1) sum_vec(i,j) = resid2_scale*ru(i,j)**2
1251  if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + resid2_scale*rv(i,j)**2
1252  enddo ; enddo
1253 
1254  dot_p1 = reproducing_sum( sum_vec, is_sum, ie_sum, js_sum, je_sum )
1255  dot_p1 = sqrt(dot_p1)
1256 
1257  if (dot_p1 <= cs%cg_tolerance * resid0) then
1258  iters = iter
1259  conv_flag = 1
1260  exit
1261  endif
1262 
1263  cg_halo = cg_halo - 1
1264 
1265  if (cg_halo == 0) then
1266  ! pass vectors
1267  call pass_vector(du, dv, g%domain, to_all, bgrid_ne)
1268  call pass_vector(u_shlf, v_shlf, g%domain, to_all, bgrid_ne)
1269  call pass_vector(ru, rv, g%domain, to_all, bgrid_ne)
1270  cg_halo = 3
1271  endif
1272 
1273  enddo ! end of CG loop
1274 
1275  do j=jsdq,jedq
1276  do i=isdq,iedq
1277  if (cs%umask(i,j) == 3) then
1278  u_shlf(i,j) = cs%u_bdry_val(i,j)
1279  elseif (cs%umask(i,j) == 0) then
1280  u_shlf(i,j) = 0
1281  endif
1282 
1283  if (cs%vmask(i,j) == 3) then
1284  v_shlf(i,j) = cs%v_bdry_val(i,j)
1285  elseif (cs%vmask(i,j) == 0) then
1286  v_shlf(i,j) = 0
1287  endif
1288  enddo
1289  enddo
1290 
1291  call pass_vector(u_shlf, v_shlf, g%domain, to_all, bgrid_ne)
1292 
1293  if (conv_flag == 0) then
1294  iters = cs%cg_max_iterations
1295  endif
1296 

◆ ice_shelf_solve_outer()

subroutine mom_ice_shelf_dynamics::ice_shelf_solve_outer ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  u_shlf,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  v_shlf,
integer, intent(out)  iters,
type(time_type), intent(in)  time 
)
private
Parameters
[in,out]csThe ice shelf dynamics control structure
[in]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in,out]u_shlfThe zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
[in,out]v_shlfThe meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
[out]itersThe number of iterations used in the solver.
[in]timeThe current model time

Definition at line 779 of file MOM_ice_shelf_dynamics.F90.

780  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
781  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
782  !! the ice-shelf state
783  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
784  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
785  real, dimension(SZDIB_(G),SZDJB_(G)), &
786  intent(inout) :: u_shlf !< The zonal ice shelf velocity at vertices [L T-1 ~> m s-1]
787  real, dimension(SZDIB_(G),SZDJB_(G)), &
788  intent(inout) :: v_shlf !< The meridional ice shelf velocity at vertices [L T-1 ~> m s-1]
789  integer, intent(out) :: iters !< The number of iterations used in the solver.
790  type(time_type), intent(in) :: Time !< The current model time
791 
792  real, dimension(SZDIB_(G),SZDJB_(G)) :: taudx, taudy ! Driving stresses at q-points [R L3 Z T-2 ~> kg m s-2]
793  real, dimension(SZDIB_(G),SZDJB_(G)) :: u_bdry_cont ! Boundary u-stress contribution [R L3 Z T-2 ~> kg m s-2]
794  real, dimension(SZDIB_(G),SZDJB_(G)) :: v_bdry_cont ! Boundary v-stress contribution [R L3 Z T-2 ~> kg m s-2]
795  real, dimension(SZDIB_(G),SZDJB_(G)) :: Au, Av ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2]
796  real, dimension(SZDIB_(G),SZDJB_(G)) :: err_u, err_v
797  real, dimension(SZDIB_(G),SZDJB_(G)) :: u_last, v_last ! Previous velocities [L T-1 ~> m s-1]
798  real, dimension(SZDIB_(G),SZDJB_(G)) :: H_node ! Ice shelf thickness at corners [Z ~> m].
799  real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! An array indicating where the ice
800  ! shelf is floating: 0 if floating, 1 if not.
801  character(len=160) :: mesg ! The text of an error message
802  integer :: conv_flag, i, j, k,l, iter
803  integer :: isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, nodefloat, nsub
804  real :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv
805  real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim]
806  real, pointer, dimension(:,:,:,:) :: Phi => null() ! The gradients of bilinear basis elements at Gaussian
807  ! quadrature points surrounding the cell vertices [m-1].
808  real, pointer, dimension(:,:,:,:,:,:) :: Phisub => null() ! Quadrature structure weights at subgridscale
809  ! locations for finite element calculations [nondim]
810  character(2) :: iternum
811  character(2) :: numproc
812 
813  ! for GL interpolation
814  nsub = cs%n_sub_regularize
815 
816  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
817  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
818  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
819 
820  taudx(:,:) = 0.0 ; taudy(:,:) = 0.0
821  u_bdry_cont(:,:) = 0.0 ; v_bdry_cont(:,:) = 0.0
822  au(:,:) = 0.0 ; av(:,:) = 0.0
823 
824  ! need to make these conditional on GL interpolation
825  float_cond(:,:) = 0.0 ; h_node(:,:) = 0.0
826  allocate(phisub(nsub,nsub,2,2,2,2)) ; phisub(:,:,:,:,:,:) = 0.0
827 
828  call calc_shelf_driving_stress(cs, iss, g, us, taudx, taudy, cs%OD_av)
829 
830  ! This is to determine which cells contain the grounding line, the criterion being that the cell
831  ! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by
832  ! assuming topography is cellwise constant and H is bilinear in a cell; floating where
833  ! rho_i/rho_w * H_node - D is negative
834 
835  ! need to make this conditional on GL interp
836 
837  if (cs%GL_regularize) then
838 
839  call interpolate_h_to_b(g, iss%h_shelf, iss%hmask, h_node)
840 
841  do j=g%jsc,g%jec ; do i=g%isc,g%iec
842  nodefloat = 0
843 
844  do l=0,1 ; do k=0,1
845  if ((iss%hmask(i,j) == 1) .and. &
846  (rhoi_rhow * h_node(i-1+k,j-1+l) - g%bathyT(i,j) <= 0)) then
847  nodefloat = nodefloat + 1
848  endif
849  enddo ; enddo
850  if ((nodefloat > 0) .and. (nodefloat < 4)) then
851  float_cond(i,j) = 1.0
852  cs%ground_frac(i,j) = 1.0
853  endif
854  enddo ; enddo
855 
856  call pass_var(float_cond, g%Domain)
857 
858  call bilinear_shape_functions_subgrid(phisub, nsub)
859 
860  endif
861 
862  ! must prepare Phi
863  allocate(phi(1:8,1:4,isd:ied,jsd:jed)) ; phi(:,:,:,:) = 0.0
864 
865  do j=jsd,jed ; do i=isd,ied
866  call bilinear_shape_fn_grid(g, i, j, phi(:,:,i,j))
867  enddo ; enddo
868 
869  call calc_shelf_visc(cs, iss, g, us, u_shlf, v_shlf)
870 
871  call pass_var(cs%ice_visc, g%domain)
872  call pass_var(cs%basal_traction, g%domain)
873 
874  ! This makes sure basal stress is only applied when it is supposed to be
875  do j=g%jsd,g%jed ; do i=g%isd,g%ied
876  cs%basal_traction(i,j) = cs%basal_traction(i,j) * cs%ground_frac(i,j)
877  enddo ; enddo
878 
879  call apply_boundary_values(cs, iss, g, us, time, phisub, h_node, cs%ice_visc, &
880  cs%basal_traction, float_cond, rhoi_rhow, u_bdry_cont, v_bdry_cont)
881 
882  au(:,:) = 0.0 ; av(:,:) = 0.0
883 
884  call cg_action(au, av, u_shlf, v_shlf, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
885  cs%ice_visc, float_cond, g%bathyT, cs%basal_traction, &
886  g, us, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
887 
888  if (cs%nonlin_solve_err_mode == 1) then
889  err_init = 0 ; err_tempu = 0 ; err_tempv = 0
890  do j=g%IscB,g%JecB ; do i=g%IscB,g%IecB
891  if (cs%umask(i,j) == 1) then
892  err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
893  if (err_tempu >= err_init) err_init = err_tempu
894  endif
895  if (cs%vmask(i,j) == 1) then
896  err_tempv = abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j))
897  if (err_tempv >= err_init) err_init = err_tempv
898  endif
899  enddo ; enddo
900 
901  call max_across_pes(err_init)
902  endif
903 
904  u_last(:,:) = u_shlf(:,:) ; v_last(:,:) = v_shlf(:,:)
905 
906  !! begin loop
907 
908  do iter=1,100
909 
910  call ice_shelf_solve_inner(cs, iss, g, us, u_shlf, v_shlf, taudx, taudy, h_node, float_cond, &
911  iss%hmask, conv_flag, iters, time, phi, phisub)
912 
913  if (cs%debug) then
914  call qchksum(u_shlf, "u shelf", g%HI, haloshift=2, scale=us%L_T_to_m_s)
915  call qchksum(v_shlf, "v shelf", g%HI, haloshift=2, scale=us%L_T_to_m_s)
916  endif
917 
918  write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations"
919  call mom_mesg(mesg, 5)
920 
921  call calc_shelf_visc(cs, iss, g, us, u_shlf, v_shlf)
922  call pass_var(cs%ice_visc, g%domain)
923  call pass_var(cs%basal_traction, g%domain)
924 
925  ! makes sure basal stress is only applied when it is supposed to be
926  do j=g%jsd,g%jed ; do i=g%isd,g%ied
927  cs%basal_traction(i,j) = cs%basal_traction(i,j) * cs%ground_frac(i,j)
928  enddo ; enddo
929 
930  u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0
931 
932  call apply_boundary_values(cs, iss, g, us, time, phisub, h_node, cs%ice_visc, &
933  cs%basal_traction, float_cond, rhoi_rhow, u_bdry_cont, v_bdry_cont)
934 
935  au(:,:) = 0 ; av(:,:) = 0
936 
937  call cg_action(au, av, u_shlf, v_shlf, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
938  cs%ice_visc, float_cond, g%bathyT, cs%basal_traction, &
939  g, us, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
940 
941  err_max = 0
942 
943  if (cs%nonlin_solve_err_mode == 1) then
944 
945  do j=g%jscB,g%jecB ; do i=g%jscB,g%iecB
946  if (cs%umask(i,j) == 1) then
947  err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
948  if (err_tempu >= err_max) err_max = err_tempu
949  endif
950  if (cs%vmask(i,j) == 1) then
951  err_tempv = abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j))
952  if (err_tempv >= err_max) err_max = err_tempv
953  endif
954  enddo ; enddo
955 
956  call max_across_pes(err_max)
957 
958  elseif (cs%nonlin_solve_err_mode == 2) then
959 
960  max_vel = 0 ; tempu = 0 ; tempv = 0
961  do j=g%jscB,g%jecB ; do i=g%iscB,g%iecB
962  if (cs%umask(i,j) == 1) then
963  err_tempu = abs(u_last(i,j)-u_shlf(i,j))
964  if (err_tempu >= err_max) err_max = err_tempu
965  tempu = u_shlf(i,j)
966  else
967  tempu = 0.0
968  endif
969  if (cs%vmask(i,j) == 1) then
970  err_tempv = max(abs(v_last(i,j)-v_shlf(i,j)), err_tempu)
971  if (err_tempv >= err_max) err_max = err_tempv
972  tempv = sqrt(v_shlf(i,j)**2 + tempu**2)
973  endif
974  if (tempv >= max_vel) max_vel = tempv
975  enddo ; enddo
976 
977  u_last(:,:) = u_shlf(:,:)
978  v_last(:,:) = v_shlf(:,:)
979 
980  call max_across_pes(max_vel)
981  call max_across_pes(err_max)
982  err_init = max_vel
983  endif
984 
985  write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init
986  call mom_mesg(mesg, 5)
987 
988  if (err_max <= cs%nonlinear_tolerance * err_init) then
989  write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations"
990  call mom_mesg(mesg, 5)
991  exit
992  endif
993 
994  enddo
995 
996  deallocate(phi)
997  deallocate(phisub)
998 

◆ ice_shelf_temp()

subroutine mom_ice_shelf_dynamics::ice_shelf_temp ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, intent(in)  time_step,
real, dimension(szdi_(g),szdj_(g)), intent(in)  melt_rate,
type(time_type), intent(in)  Time 
)
private

This subroutine updates the vertically averaged ice shelf temperature.

Parameters
[in,out]csA pointer to the ice shelf control structure
[in]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]time_stepThe time step for this update [T ~> s].
[in]melt_ratebasal melt rate [R Z T-1 ~> kg m-2 s-1]
[in]timeThe current model time

Definition at line 2975 of file MOM_ice_shelf_dynamics.F90.

2976  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
2977  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
2978  !! the ice-shelf state
2979  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
2980  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
2981  real, intent(in) :: time_step !< The time step for this update [T ~> s].
2982  real, dimension(SZDI_(G),SZDJ_(G)), &
2983  intent(in) :: melt_rate !< basal melt rate [R Z T-1 ~> kg m-2 s-1]
2984  type(time_type), intent(in) :: Time !< The current model time
2985 
2986 ! 5/23/12 OVS
2987 ! This subroutine takes the velocity (on the Bgrid) and timesteps
2988 ! (HT)_t = - div (uHT) + (adot Tsurf -bdot Tbot) once and then calculates T=HT/H
2989 !
2990 ! The flux overflows are included here. That is because they will be used to advect 3D scalars
2991 ! into partial cells
2992 
2993  real, dimension(SZDI_(G),SZDJ_(G)) :: th_after_uflux, th_after_vflux, TH
2994  integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
2995  real :: Tsurf ! Surface air temperature. This is hard coded but should be an input argument.
2996  real :: adot ! A surface heat exchange coefficient divided by the heat capacity of
2997  ! ice [R Z T-1 degC-1 ~> kg m-2 s-1 degC-1].
2998 
2999 
3000  ! For now adot and Tsurf are defined here adot=surf acc 0.1m/yr, Tsurf=-20oC, vary them later
3001  adot = (0.1/(365.0*86400.0))*us%m_to_Z*us%T_to_s * cs%density_ice
3002  tsurf = -20.0
3003 
3004  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3005  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
3006 
3007  th_after_uflux(:,:) = 0.0
3008  th_after_vflux(:,:) = 0.0
3009 
3010  do j=jsd,jed ; do i=isd,ied
3011 ! if (ISS%hmask(i,j) > 1) then
3012  if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2)) then
3013  cs%t_shelf(i,j) = cs%t_bdry_val(i,j)
3014  endif
3015  enddo ; enddo
3016 
3017  do j=jsd,jed ; do i=isd,ied
3018  ! Convert the averge temperature to a depth integrated temperature.
3019  th(i,j) = cs%t_shelf(i,j)*iss%h_shelf(i,j)
3020  enddo ; enddo
3021 
3022 ! call enable_averages(time_step, Time, CS%diag)
3023 ! call pass_var(h_after_uflux, G%domain)
3024 ! call pass_var(h_after_vflux, G%domain)
3025 ! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag)
3026 ! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag)
3027 ! call disable_averaging(CS%diag)
3028 
3029  call ice_shelf_advect_temp_x(cs, g, time_step, iss%hmask, th, th_after_uflux)
3030  call ice_shelf_advect_temp_y(cs, g, time_step, iss%hmask, th_after_uflux, th_after_vflux)
3031 
3032  do j=jsc,jec ; do i=isc,iec
3033  ! Convert the integrated temperature back to the average temperature.
3034 ! if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 2)) then
3035  if (iss%h_shelf(i,j) > 0.0) then
3036  cs%t_shelf(i,j) = th_after_vflux(i,j) / iss%h_shelf(i,j)
3037  else
3038  cs%t_shelf(i,j) = -10.0
3039  endif
3040 ! endif
3041 
3042  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
3043  if (iss%h_shelf(i,j) > 0.0) then
3044  cs%t_shelf(i,j) = cs%t_shelf(i,j) + &
3045  time_step*(adot*tsurf - melt_rate(i,j)*iss%tfreeze(i,j))/(cs%density_ice*iss%h_shelf(i,j))
3046  else
3047  ! the ice is about to melt away in this case set thickness, area, and mask to zero
3048  ! NOTE: not mass conservative, should maybe scale salt & heat flux for this cell
3049  cs%t_shelf(i,j) = -10.0
3050  cs%tmask(i,j) = 0.0
3051  endif
3052  elseif (iss%hmask(i,j) == 0) then
3053  cs%t_shelf(i,j) = -10.0
3054  elseif ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2)) then
3055  cs%t_shelf(i,j) = cs%t_bdry_val(i,j)
3056  endif
3057  enddo ; enddo
3058 
3059  call pass_var(cs%t_shelf, g%domain)
3060  call pass_var(cs%tmask, g%domain)
3061 
3062  if (cs%debug) then
3063  call hchksum(cs%t_shelf, "temp after front", g%HI, haloshift=3)
3064  endif
3065 

◆ ice_time_step_cfl()

real function, public mom_ice_shelf_dynamics::ice_time_step_cfl ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(inout)  ISS,
type(ocean_grid_type), intent(inout)  G 
)

This function returns the global maximum advective timestep that can be taken based on the current ice velocities. Because it involves finding a global minimum, it can be surprisingly expensive.

Parameters
[in,out]csThe ice shelf dynamics control structure
[in,out]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
Returns
The maximum permitted timestep based on the ice velocities [T ~> s].

Definition at line 600 of file MOM_ice_shelf_dynamics.F90.

601  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
602  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
603  !! the ice-shelf state
604  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
605  real :: ice_time_step_CFL !< The maximum permitted timestep based on the ice velocities [T ~> s].
606 
607  real :: dt_local, min_dt ! These should be the minimum stable timesteps at a CFL of 1 [T ~> s]
608  real :: min_vel ! A minimal velocity for estimating a timestep [L T-1 ~> m s-1]
609  integer :: i, j
610 
611  min_dt = 5.0e17*g%US%s_to_T ! The starting maximum is roughly the lifetime of the universe.
612  min_vel = (1.0e-12/(365.0*86400.0)) * g%US%m_s_to_L_T
613  do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (iss%hmask(i,j) == 1.0) then
614  dt_local = 2.0*g%areaT(i,j) / &
615  ((g%dyCu(i,j) * max(abs(cs%u_shelf(i,j) + cs%u_shelf(i,j-1)), min_vel) + &
616  g%dyCu(i-1,j)* max(abs(cs%u_shelf(i-1,j)+ cs%u_shelf(i-1,j-1)), min_vel)) + &
617  (g%dxCv(i,j) * max(abs(cs%v_shelf(i,j) + cs%v_shelf(i-1,j)), min_vel) + &
618  g%dxCv(i,j-1)* max(abs(cs%v_shelf(i,j-1)+ cs%v_shelf(i-1,j-1)), min_vel)))
619 
620  min_dt = min(min_dt, dt_local)
621  endif ; enddo ; enddo ! i- and j- loops
622 
623  call min_across_pes(min_dt)
624 
625  ice_time_step_cfl = cs%CFL_factor * min_dt
626 

◆ init_boundary_values()

subroutine mom_ice_shelf_dynamics::init_boundary_values ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ocean_grid_type), intent(inout)  G,
type(time_type), intent(in)  time,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, intent(in)  input_flux,
real, intent(in)  input_thick,
logical, intent(in), optional  new_sim 
)
private
Parameters
[in,out]csA pointer to the ice shelf control structure
[in,out]gThe grid structure used by the ice shelf.
[in]timeThe current model time
[in]hmaskA mask indicating which tracer points are
[in]input_fluxThe integrated inward ice thickness flux per unit face length [Z L T-1 ~> m2 s-1]
[in]input_thickThe ice thickness at boundaries [Z ~> m].
[in]new_simIf present and false, this run is being restarted

Definition at line 1902 of file MOM_ice_shelf_dynamics.F90.

1903  type(ice_shelf_dyn_CS),intent(inout) :: CS !< A pointer to the ice shelf control structure
1904  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
1905  type(time_type), intent(in) :: Time !< The current model time
1906  real, dimension(SZDI_(G),SZDJ_(G)), &
1907  intent(in) :: hmask !< A mask indicating which tracer points are
1908  !! partly or fully covered by an ice-shelf
1909  real, intent(in) :: input_flux !< The integrated inward ice thickness flux per
1910  !! unit face length [Z L T-1 ~> m2 s-1]
1911  real, intent(in) :: input_thick !< The ice thickness at boundaries [Z ~> m].
1912  logical, optional, intent(in) :: new_sim !< If present and false, this run is being restarted
1913 
1914 ! this will be a per-setup function. the boundary values of thickness and velocity
1915 ! (and possibly other variables) will be updated in this function
1916 
1917 ! FOR RESTARTING PURPOSES: if grid is not symmetric and the model is restarted, we will
1918 ! need to update those velocity points not *technically* in any
1919 ! computational domain -- if this function gets moves to another module,
1920 ! DO NOT TAKE THE RESTARTING BIT WITH IT
1921  integer :: i, j , isd, jsd, ied, jed
1922  integer :: gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
1923  integer :: i_off, j_off
1924 
1925  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
1926  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1927  i_off = g%idg_offset ; j_off = g%jdg_offset
1928 
1929  ! this loop results in some values being set twice but... eh.
1930 
1931  do j=jsd,jed
1932  do i=isd,ied
1933 
1934  if (hmask(i,j) == 3) then
1935  cs%thickness_bdry_val(i,j) = input_thick
1936  endif
1937 
1938  if ((hmask(i,j) == 0) .or. (hmask(i,j) == 1) .or. (hmask(i,j) == 2)) then
1939  if ((i <= iec).and.(i >= isc)) then
1940  if (cs%u_face_mask(i-1,j) == 3) then
1941  cs%u_bdry_val(i-1,j-1) = (1 - ((g%geoLatBu(i-1,j-1) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
1942  1.5 * input_flux / input_thick
1943  cs%u_bdry_val(i-1,j) = (1 - ((g%geoLatBu(i-1,j) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
1944  1.5 * input_flux / input_thick
1945  endif
1946  endif
1947  endif
1948 
1949  if (.not.(new_sim)) then
1950  if (.not. g%symmetric) then
1951  if (((i+i_off) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3)) then
1952  cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
1953  cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
1954  cs%v_shelf(i-1,j-1) = cs%v_bdry_val(i-1,j-1)
1955  cs%v_shelf(i-1,j) = cs%v_bdry_val(i-1,j)
1956  endif
1957  if (((j+j_off) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3)) then
1958  cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
1959  cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
1960  cs%v_shelf(i-1,j-1) = cs%v_bdry_val(i-1,j-1)
1961  cs%v_shelf(i,j-1) = cs%v_bdry_val(i,j-1)
1962  endif
1963  endif
1964  endif
1965  enddo
1966  enddo
1967 

◆ initialize_diagnostic_fields()

subroutine mom_ice_shelf_dynamics::initialize_diagnostic_fields ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(in)  ISS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(time_type), intent(in)  Time 
)
private
Parameters
[in,out]csA pointer to the ice shelf control structure
[in]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]timeThe current model time

Definition at line 563 of file MOM_ice_shelf_dynamics.F90.

564  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
565  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
566  !! the ice-shelf state
567  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
568  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
569  type(time_type), intent(in) :: Time !< The current model time
570 
571  integer :: i, j, iters, isd, ied, jsd, jed
572  real :: rhoi_rhow
573  real :: OD ! Depth of open water below the ice shelf [Z ~> m]
574  type(time_type) :: dummy_time
575 
576  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
577  dummy_time = set_time(0,0)
578  isd=g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
579 
580  do j=jsd,jed
581  do i=isd,ied
582  od = g%bathyT(i,j) - rhoi_rhow * iss%h_shelf(i,j)
583  if (od >= 0) then
584  ! ice thickness does not take up whole ocean column -> floating
585  cs%OD_av(i,j) = od
586  cs%ground_frac(i,j) = 0.
587  else
588  cs%OD_av(i,j) = 0.
589  cs%ground_frac(i,j) = 1.
590  endif
591  enddo
592  enddo
593 
594  call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, dummy_time)
595 

◆ initialize_ice_shelf_dyn()

subroutine, public mom_ice_shelf_dynamics::initialize_ice_shelf_dyn ( type(param_file_type), intent(in)  param_file,
type(time_type), intent(inout)  Time,
type(ice_shelf_state), intent(in)  ISS,
type(ice_shelf_dyn_cs), pointer  CS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(diag_ctrl), intent(in), target  diag,
logical, intent(in)  new_sim,
logical, intent(in), optional  solo_ice_sheet_in 
)

Initializes shelf model data, parameters and diagnostics.

Parameters
[in]param_fileA structure to parse for run-time parameters
[in,out]timeThe clock that that will indicate the model time
[in]issA structure with elements that describe the ice-shelf state
csA pointer to the ice shelf dynamics control structure
[in,out]gThe grid type describing the ice shelf grid.
[in]usA structure containing unit conversion factors
[in]diagA structure that is used to regulate the diagnostic output.
[in]new_simIf true this is a new simulation, otherwise has been started from a restart file.
[in]solo_ice_sheet_inIf present, this indicates whether a solo ice-sheet driver.

Definition at line 273 of file MOM_ice_shelf_dynamics.F90.

274  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
275  type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure
276  type(time_type), intent(inout) :: Time !< The clock that that will indicate the model time
277  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
278  !! the ice-shelf state
279  type(ice_shelf_dyn_CS), pointer :: CS !< A pointer to the ice shelf dynamics control structure
280  type(ocean_grid_type), intent(inout) :: G !< The grid type describing the ice shelf grid.
281  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
282  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate the diagnostic output.
283  logical, intent(in) :: new_sim !< If true this is a new simulation, otherwise
284  !! has been started from a restart file.
285  logical, optional, intent(in) :: solo_ice_sheet_in !< If present, this indicates whether
286  !! a solo ice-sheet driver.
287 
288  ! Local variables
289  real :: Z_rescale ! A rescaling factor for heights from the representation in
290  ! a restart file to the internal representation in this run.
291  real :: vel_rescale ! A rescaling factor for horizontal velocities from the representation
292  ! in a restart file to the internal representation in this run.
293  !This include declares and sets the variable "version".
294 # include "version_variable.h"
295  character(len=200) :: config
296  character(len=200) :: IC_file,filename,inputdir
297  character(len=40) :: var_name
298  character(len=40) :: mdl = "MOM_ice_shelf_dyn" ! This module's name.
299  logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
300  logical :: debug
301  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq, iters
302 
303  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
304  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
305 
306  if (.not.associated(cs)) then
307  call mom_error(fatal, "MOM_ice_shelf_dyn.F90, initialize_ice_shelf_dyn: "// &
308  "called with an associated control structure.")
309  return
310  endif
311  if (cs%module_is_initialized) then
312  call mom_error(warning, "MOM_ice_shelf_dyn.F90, initialize_ice_shelf_dyn was "//&
313  "called with a control structure that has already been initialized.")
314  endif
315  cs%module_is_initialized = .true.
316 
317  cs%diag => diag ! ; CS%Time => Time
318 
319  ! Read all relevant parameters and write them to the model log.
320  call log_version(param_file, mdl, version, "")
321  call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
322  call get_param(param_file, mdl, "DEBUG_IS", cs%debug, &
323  "If true, write verbose debugging messages for the ice shelf.", &
324  default=debug)
325  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
326  "If true, the ice sheet mass can evolve with time.", &
327  default=.false.)
328  override_shelf_movement = .false. ; active_shelf_dynamics = .false.
329  if (shelf_mass_is_dynamic) then
330  call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", override_shelf_movement, &
331  "If true, user provided code specifies the ice-shelf "//&
332  "movement instead of the dynamic ice model.", default=.false., do_not_log=.true.)
333  active_shelf_dynamics = .not.override_shelf_movement
334 
335  call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
336  "If true, regularize the floatation condition at the "//&
337  "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
338  call get_param(param_file, mdl, "GROUNDING_LINE_INTERP_SUBGRID_N", cs%n_sub_regularize, &
339  "The number of sub-partitions of each cell over which to "//&
340  "integrate for the interpolated grounding line. Each cell "//&
341  "is divided into NxN equally-sized rectangles, over which the "//&
342  "basal contribution is integrated by iterative quadrature.", &
343  default=0)
344  call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", cs%GL_couple, &
345  "If true, let the floatation condition be determined by "//&
346  "ocean column thickness. This means that update_OD_ffrac "//&
347  "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", &
348  default=.false., do_not_log=cs%GL_regularize)
349  if (cs%GL_regularize) cs%GL_couple = .false.
350  if (cs%GL_regularize .and. (cs%n_sub_regularize == 0)) call mom_error (fatal, &
351  "GROUNDING_LINE_INTERP_SUBGRID_N must be a positive integer if GL regularization is used")
352  call get_param(param_file, mdl, "ICE_SHELF_CFL_FACTOR", cs%CFL_factor, &
353  "A factor used to limit timestep as CFL_FACTOR * min (\Delta x / u). "//&
354  "This is only used with an ice-only model.", default=0.25)
355  endif
356  call get_param(param_file, mdl, "RHO_0", cs%density_ocean_avg, &
357  "avg ocean density used in floatation cond", &
358  units="kg m-3", default=1035., scale=us%kg_m3_to_R)
359  if (active_shelf_dynamics) then
360  call get_param(param_file, mdl, "ICE_VELOCITY_TIMESTEP", cs%velocity_update_time_step, &
361  "seconds between ice velocity calcs", units="s", scale=us%s_to_T, &
362  fail_if_missing=.true.)
363  call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
364  "The gravitational acceleration of the Earth.", &
365  units="m s-2", default = 9.80, scale=us%m_s_to_L_T**2*us%Z_to_m)
366 
367  call get_param(param_file, mdl, "A_GLEN_ISOTHERM", cs%A_glen_isothermal, &
368  "Ice viscosity parameter in Glen's Law", &
369  units="Pa-3 yr-1", default=9.461e-18, scale=1.0/(365.0*86400.0))
370  ! This default is equivalent to 3.0001e-25 Pa-3 s-1, appropriate at about -10 C.
371  call get_param(param_file, mdl, "GLEN_EXPONENT", cs%n_glen, &
372  "nonlinearity exponent in Glen's Law", &
373  units="none", default=3.)
374  call get_param(param_file, mdl, "MIN_STRAIN_RATE_GLEN", cs%eps_glen_min, &
375  "min. strain rate to avoid infinite Glen's law viscosity", &
376  units="a-1", default=1.e-12, scale=us%T_to_s/(365.0*86400.0))
377  call get_param(param_file, mdl, "BASAL_FRICTION_EXP", cs%n_basal_fric, &
378  "Exponent in sliding law \tau_b = C u^(n_basal_fric)", &
379  units="none", fail_if_missing=.true.)
380  call get_param(param_file, mdl, "BASAL_FRICTION_COEFF", cs%C_basal_friction, &
381  "Coefficient in sliding law \tau_b = C u^(n_basal_fric)", &
382  units="Pa (m yr-1)-(n_basal_fric)", scale=us%kg_m2s_to_RZ_T*((365.0*86400.0)**cs%n_basal_fric), &
383  fail_if_missing=.true.)
384  call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
385  "A typical density of ice.", units="kg m-3", default=917.0, scale=us%kg_m3_to_R)
386  call get_param(param_file, mdl, "CONJUGATE_GRADIENT_TOLERANCE", cs%cg_tolerance, &
387  "tolerance in CG solver, relative to initial residual", default=1.e-6)
388  call get_param(param_file, mdl, "ICE_NONLINEAR_TOLERANCE", cs%nonlinear_tolerance, &
389  "nonlin tolerance in iterative velocity solve",default=1.e-6)
390  call get_param(param_file, mdl, "CONJUGATE_GRADIENT_MAXIT", cs%cg_max_iterations, &
391  "max iteratiions in CG solver", default=2000)
392  call get_param(param_file, mdl, "THRESH_FLOAT_COL_DEPTH", cs%thresh_float_col_depth, &
393  "min ocean thickness to consider ice *floating*; "//&
394  "will only be important with use of tides", &
395  units="m", default=1.e-3, scale=us%m_to_Z)
396  call get_param(param_file, mdl, "NONLIN_SOLVE_ERR_MODE", cs%nonlin_solve_err_mode, &
397  "Choose whether nonlin error in vel solve is based on nonlinear "//&
398  "residual (1) or relative change since last iteration (2)", default=1)
399 
400  call get_param(param_file, mdl, "SHELF_MOVING_FRONT", cs%moving_shelf_front, &
401  "Specify whether to advance shelf front (and calve).", &
402  default=.true.)
403  call get_param(param_file, mdl, "CALVE_TO_MASK", cs%calve_to_mask, &
404  "If true, do not allow an ice shelf where prohibited by a mask.", &
405  default=.false.)
406  endif
407  call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", cs%min_thickness_simple_calve, &
408  "Min thickness rule for the VERY simple calving law",&
409  units="m", default=0.0, scale=us%m_to_Z)
410 
411  ! Allocate memory in the ice shelf dynamics control structure that was not
412  ! previously allocated for registration for restarts.
413  ! OVS vertically integrated Temperature
414 
415  if (active_shelf_dynamics) then
416  ! DNG
417  allocate( cs%u_bdry_val(isdq:iedq,jsdq:jedq) ) ; cs%u_bdry_val(:,:) = 0.0
418  allocate( cs%v_bdry_val(isdq:iedq,jsdq:jedq) ) ; cs%v_bdry_val(:,:) = 0.0
419  allocate( cs%t_bdry_val(isd:ied,jsd:jed) ) ; cs%t_bdry_val(:,:) = -15.0
420  allocate( cs%h_bdry_val(isd:ied,jsd:jed) ) ; cs%h_bdry_val(:,:) = 0.0
421  allocate( cs%thickness_bdry_val(isd:ied,jsd:jed) ) ; cs%thickness_bdry_val(:,:) = 0.0
422  allocate( cs%u_face_mask(isdq:iedq,jsd:jed) ) ; cs%u_face_mask(:,:) = 0.0
423  allocate( cs%v_face_mask(isd:ied,jsdq:jedq) ) ; cs%v_face_mask(:,:) = 0.0
424  allocate( cs%u_face_mask_bdry(isdq:iedq,jsd:jed) ) ; cs%u_face_mask_bdry(:,:) = -2.0
425  allocate( cs%v_face_mask_bdry(isd:ied,jsdq:jedq) ) ; cs%v_face_mask_bdry(:,:) = -2.0
426  allocate( cs%u_flux_bdry_val(isdq:iedq,jsd:jed) ) ; cs%u_flux_bdry_val(:,:) = 0.0
427  allocate( cs%v_flux_bdry_val(isd:ied,jsdq:jedq) ) ; cs%v_flux_bdry_val(:,:) = 0.0
428  allocate( cs%umask(isdq:iedq,jsdq:jedq) ) ; cs%umask(:,:) = -1.0
429  allocate( cs%vmask(isdq:iedq,jsdq:jedq) ) ; cs%vmask(:,:) = -1.0
430  allocate( cs%tmask(isdq:iedq,jsdq:jedq) ) ; cs%tmask(:,:) = -1.0
431 
432  cs%OD_rt_counter = 0
433  allocate( cs%OD_rt(isd:ied,jsd:jed) ) ; cs%OD_rt(:,:) = 0.0
434  allocate( cs%ground_frac_rt(isd:ied,jsd:jed) ) ; cs%ground_frac_rt(:,:) = 0.0
435 
436  if (cs%calve_to_mask) then
437  allocate( cs%calve_mask(isd:ied,jsd:jed) ) ; cs%calve_mask(:,:) = 0.0
438  endif
439 
440  cs%elapsed_velocity_time = 0.0
441 
442  call update_velocity_masks(cs, g, iss%hmask, cs%umask, cs%vmask, cs%u_face_mask, cs%v_face_mask)
443  endif
444 
445  ! Take additional initialization steps, for example of dependent variables.
446  if (active_shelf_dynamics .and. .not.new_sim) then
447  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) then
448  z_rescale = us%m_to_Z / us%m_to_Z_restart
449  do j=g%jsc,g%jec ; do i=g%isc,g%iec
450  cs%OD_av(i,j) = z_rescale * cs%OD_av(i,j)
451  enddo ; enddo
452  endif
453 
454  if ((us%m_to_L_restart*us%s_to_T_restart /= 0.0) .and. &
455  (us%m_to_L_restart /= us%m_s_to_L_T*us%s_to_T_restart)) then
456  vel_rescale = us%m_s_to_L_T*us%s_to_T_restart / us%m_to_L_restart
457  do j=g%jsc-1,g%jec ; do i=g%isc-1,g%iec
458  cs%u_shelf(i,j) = vel_rescale * cs%u_shelf(i,j)
459  cs%v_shelf(i,j) = vel_rescale * cs%v_shelf(i,j)
460  enddo ; enddo
461  endif
462 
463  ! this is unfortunately necessary; if grid is not symmetric the boundary values
464  ! of u and v are otherwise not set till the end of the first linear solve, and so
465  ! viscosity is not calculated correctly.
466  ! This has to occur after init_boundary_values or some of the arrays on the
467  ! right hand side have not been set up yet.
468  if (.not. g%symmetric) then
469  do j=g%jsd,g%jed ; do i=g%isd,g%ied
470  if (((i+g%idg_offset) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3)) then
471  cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
472  cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
473  cs%v_shelf(i-1,j-1) = cs%v_bdry_val(i-1,j-1)
474  cs%v_shelf(i-1,j) = cs%v_bdry_val(i-1,j)
475  endif
476  if (((j+g%jdg_offset) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3)) then
477  cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
478  cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
479  cs%v_shelf(i-1,j-1) = cs%v_bdry_val(i-1,j-1)
480  cs%v_shelf(i,j-1) = cs%v_bdry_val(i,j-1)
481  endif
482  enddo ; enddo
483  endif
484 
485  call pass_var(cs%OD_av,g%domain)
486  call pass_var(cs%ground_frac,g%domain)
487  call pass_var(cs%ice_visc,g%domain)
488  call pass_var(cs%basal_traction, g%domain)
489  call pass_vector(cs%u_shelf, cs%v_shelf, g%domain, to_all, bgrid_ne)
490  endif
491 
492  if (active_shelf_dynamics) then
493  ! If we are calving to a mask, i.e. if a mask exists where a shelf cannot, read the mask from a file.
494  if (cs%calve_to_mask) then
495  call mom_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading calving_mask")
496 
497  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
498  inputdir = slasher(inputdir)
499  call get_param(param_file, mdl, "CALVING_MASK_FILE", ic_file, &
500  "The file with a mask for where calving might occur.", &
501  default="ice_shelf_h.nc")
502  call get_param(param_file, mdl, "CALVING_MASK_VARNAME", var_name, &
503  "The variable to use in masking calving.", &
504  default="area_shelf_h")
505 
506  filename = trim(inputdir)//trim(ic_file)
507  call log_param(param_file, mdl, "INPUTDIR/CALVING_MASK_FILE", filename)
508  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
509  " calving mask file: Unable to open "//trim(filename))
510 
511  call mom_read_data(filename,trim(var_name),cs%calve_mask,g%Domain)
512  do j=g%jsc,g%jec ; do i=g%isc,g%iec
513  if (cs%calve_mask(i,j) > 0.0) cs%calve_mask(i,j) = 1.0
514  enddo ; enddo
515  call pass_var(cs%calve_mask,g%domain)
516  endif
517 
518 ! call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim)
519 
520  if (new_sim) then
521  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.")
522  call update_od_ffrac_uncoupled(cs, g, iss%h_shelf(:,:))
523  call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, time)
524 
525  if (cs%id_u_shelf > 0) call post_data(cs%id_u_shelf, cs%u_shelf, cs%diag)
526  if (cs%id_v_shelf > 0) call post_data(cs%id_v_shelf, cs%v_shelf,cs%diag)
527  endif
528 
529  ! Register diagnostics.
530  cs%id_u_shelf = register_diag_field('ocean_model','u_shelf',cs%diag%axesCu1, time, &
531  'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*us%L_T_to_m_s)
532  cs%id_v_shelf = register_diag_field('ocean_model','v_shelf',cs%diag%axesCv1, time, &
533  'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*us%L_T_to_m_s)
534  cs%id_u_mask = register_diag_field('ocean_model','u_mask',cs%diag%axesCu1, time, &
535  'mask for u-nodes', 'none')
536  cs%id_v_mask = register_diag_field('ocean_model','v_mask',cs%diag%axesCv1, time, &
537  'mask for v-nodes', 'none')
538 ! CS%id_surf_elev = register_diag_field('ocean_model','ice_surf',CS%diag%axesT1, Time, &
539 ! 'ice surf elev', 'm')
540  cs%id_ground_frac = register_diag_field('ocean_model','ice_ground_frac',cs%diag%axesT1, time, &
541  'fraction of cell that is grounded', 'none')
542  cs%id_col_thick = register_diag_field('ocean_model','col_thick',cs%diag%axesT1, time, &
543  'ocean column thickness passed to ice model', 'm', conversion=us%Z_to_m)
544  cs%id_OD_av = register_diag_field('ocean_model','OD_av',cs%diag%axesT1, time, &
545  'intermediate ocean column thickness passed to ice model', 'm', conversion=us%Z_to_m)
546  !CS%id_h_after_uflux = register_diag_field('ocean_model','h_after_uflux',CS%diag%axesh1, Time, &
547  ! 'thickness after u flux ', 'none')
548  !CS%id_h_after_vflux = register_diag_field('ocean_model','h_after_vflux',CS%diag%axesh1, Time, &
549  ! 'thickness after v flux ', 'none')
550  !CS%id_h_after_adv = register_diag_field('ocean_model','h_after_adv',CS%diag%axesh1, Time, &
551  ! 'thickness after front adv ', 'none')
552 
553 !!! OVS vertically integrated temperature
554  cs%id_t_shelf = register_diag_field('ocean_model','t_shelf',cs%diag%axesT1, time, &
555  'T of ice', 'oC')
556  cs%id_t_mask = register_diag_field('ocean_model','tmask',cs%diag%axesT1, time, &
557  'mask for T-nodes', 'none')
558  endif
559 

◆ interpolate_h_to_b()

subroutine mom_ice_shelf_dynamics::interpolate_h_to_b ( type(ocean_grid_type), intent(inout)  G,
real, dimension(szdi_(g),szdj_(g)), intent(in)  h_shelf,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  H_node 
)
private

Interpolate the ice shelf thickness from tracer point to nodal points, subject to a mask.

Parameters
[in,out]gThe grid structure used by the ice shelf.
[in]h_shelfThe ice shelf thickness at tracer points [Z ~> m].
[in]hmaskA mask indicating which tracer points are
[in,out]h_nodeThe ice shelf thickness at nodal (corner)

Definition at line 2910 of file MOM_ice_shelf_dynamics.F90.

2911  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
2912  real, dimension(SZDI_(G),SZDJ_(G)), &
2913  intent(in) :: h_shelf !< The ice shelf thickness at tracer points [Z ~> m].
2914  real, dimension(SZDI_(G),SZDJ_(G)), &
2915  intent(in) :: hmask !< A mask indicating which tracer points are
2916  !! partly or fully covered by an ice-shelf
2917  real, dimension(SZDIB_(G),SZDJB_(G)), &
2918  intent(inout) :: H_node !< The ice shelf thickness at nodal (corner)
2919  !! points [Z ~> m].
2920 
2921  integer :: i, j, isc, iec, jsc, jec, num_h, k, l
2922  real :: summ
2923 
2924  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2925 
2926  h_node(:,:) = 0.0
2927 
2928  ! H_node is node-centered; average over all cells that share that node
2929  ! if no (active) cells share the node then its value there is irrelevant
2930 
2931  do j=jsc-1,jec
2932  do i=isc-1,iec
2933  summ = 0.0
2934  num_h = 0
2935  do k=0,1
2936  do l=0,1
2937  if (hmask(i+k,j+l) == 1.0) then
2938  summ = summ + h_shelf(i+k,j+l)
2939  num_h = num_h + 1
2940  endif
2941  enddo
2942  enddo
2943  if (num_h > 0) then
2944  h_node(i,j) = summ / num_h
2945  endif
2946  enddo
2947  enddo
2948 
2949  call pass_var(h_node, g%domain, position=corner)
2950 

◆ matrix_diagonal()

subroutine mom_ice_shelf_dynamics::matrix_diagonal ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
real, dimension(szdi_(g),szdj_(g)), intent(in)  float_cond,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  H_node,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  ice_visc,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  basal_trac,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, intent(in)  dens_ratio,
real, dimension(:,:,:,:,:,:), intent(in)  Phisub,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  u_diagonal,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  v_diagonal 
)
private

returns the diagonal entries of the matrix for a Jacobi preconditioning

Parameters
[in]csA pointer to the ice shelf control structure
[in]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]float_condAn array indicating where the ice
[in]h_nodeThe ice shelf thickness at nodal
[in]ice_viscA field related to the ice viscosity from Glen's
[in]basal_tracA field related to the nonlinear part of the
[in]hmaskA mask indicating which tracer points are
[in]dens_ratioThe density of ice divided by the density of seawater [nondim]
[in]phisubQuadrature structure weights at subgridscale locations for finite element calculations [nondim]
[in,out]u_diagonalThe diagonal elements of the u-velocity
[in,out]v_diagonalThe diagonal elements of the v-velocity

Definition at line 2163 of file MOM_ice_shelf_dynamics.F90.

2165 
2166  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
2167  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2168  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
2169  real, dimension(SZDI_(G),SZDJ_(G)), &
2170  intent(in) :: float_cond !< An array indicating where the ice
2171  !! shelf is floating: 0 if floating, 1 if not.
2172  real, dimension(SZDIB_(G),SZDJB_(G)), &
2173  intent(in) :: H_node !< The ice shelf thickness at nodal
2174  !! (corner) points [Z ~> m].
2175  real, dimension(SZDIB_(G),SZDJB_(G)), &
2176  intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's
2177  !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form
2178  !! and units depend on the basal law exponent.
2179  real, dimension(SZDIB_(G),SZDJB_(G)), &
2180  intent(in) :: basal_trac !< A field related to the nonlinear part of the
2181  !! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1].
2182  real, dimension(SZDI_(G),SZDJ_(G)), &
2183  intent(in) :: hmask !< A mask indicating which tracer points are
2184  !! partly or fully covered by an ice-shelf
2185  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2186  !! of seawater [nondim]
2187  real, dimension(:,:,:,:,:,:), intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2188  !! locations for finite element calculations [nondim]
2189  real, dimension(SZDIB_(G),SZDJB_(G)), &
2190  intent(inout) :: u_diagonal !< The diagonal elements of the u-velocity
2191  !! matrix from the left-hand side of the solver [R L2 Z T-1 ~> kg s-1]
2192  real, dimension(SZDIB_(G),SZDJB_(G)), &
2193  intent(inout) :: v_diagonal !< The diagonal elements of the v-velocity
2194  !! matrix from the left-hand side of the solver [R L2 Z T-1 ~> kg s-1]
2195 
2196 
2197 ! returns the diagonal entries of the matrix for a Jacobi preconditioning
2198 
2199  real :: ux, uy, vx, vy ! Interpolated weight gradients [L-1 ~> m-1]
2200  real :: uq, vq
2201  real, dimension(8,4) :: Phi ! Weight gradients [L-1 ~> m-1]
2202  real, dimension(2) :: xquad
2203  real, dimension(2,2) :: Hcell, sub_ground
2204  integer :: i, j, is, js, cnt, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq, Itgt, Jtgt
2205 
2206  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2207 
2208  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2209 
2210  do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1) then
2211 
2212  call bilinear_shape_fn_grid(g, i, j, phi)
2213 
2214  ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2215  ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j
2216 
2217  do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; itgt = i-2+iphi ; jtgt = j-2-jphi
2218  ilq = 1 ; if (iq == iphi) ilq = 2
2219  jlq = 1 ; if (jq == jphi) jlq = 2
2220 
2221  if (cs%umask(itgt,jtgt) == 1) then
2222 
2223  ux = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
2224  uy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
2225  vx = 0.
2226  vy = 0.
2227 
2228  u_diagonal(itgt,jtgt) = u_diagonal(itgt,jtgt) + &
2229  0.25 * ice_visc(i,j) * ((4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2230  (uy+vy) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2231 
2232  if (float_cond(i,j) == 0) then
2233  uq = xquad(ilq) * xquad(jlq)
2234  u_diagonal(itgt,jtgt) = u_diagonal(itgt,jtgt) + &
2235  0.25 * basal_trac(i,j) * uq * xquad(ilq) * xquad(jlq)
2236  endif
2237  endif
2238 
2239  if (cs%vmask(itgt,jtgt) == 1) then
2240 
2241  vx = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
2242  vy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
2243  ux = 0.
2244  uy = 0.
2245 
2246  v_diagonal(itgt,jtgt) = v_diagonal(itgt,jtgt) + &
2247  0.25 * ice_visc(i,j) * ((uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2248  (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2249 
2250  if (float_cond(i,j) == 0) then
2251  vq = xquad(ilq) * xquad(jlq)
2252  v_diagonal(itgt,jtgt) = v_diagonal(itgt,jtgt) + &
2253  0.25 * basal_trac(i,j) * vq * xquad(ilq) * xquad(jlq)
2254  endif
2255  endif
2256  enddo ; enddo ; enddo ; enddo
2257 
2258  if (float_cond(i,j) == 1) then
2259  hcell(:,:) = h_node(i-1:i,j-1:j)
2260  call cg_diagonal_subgrid_basal(phisub, hcell, g%bathyT(i,j), dens_ratio, sub_ground)
2261  do iphi=1,2 ; do jphi=1,2 ; itgt = i-2+iphi ; jtgt = j-2-jphi
2262  if (cs%umask(itgt,jtgt) == 1) then
2263  u_diagonal(itgt,jtgt) = u_diagonal(itgt,jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j)
2264  v_diagonal(itgt,jtgt) = v_diagonal(itgt,jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j)
2265  endif
2266  enddo ; enddo
2267  endif
2268  endif ; enddo ; enddo
2269 

◆ quad_area()

real function mom_ice_shelf_dynamics::quad_area ( real, dimension(4), intent(in)  X,
real, dimension(4), intent(in)  Y 
)
private

Calculate area of quadrilateral.

Parameters
[in]xThe x-positions of the vertices of the quadrilateral [L ~> m].
[in]yThe y-positions of the vertices of the quadrilateral [L ~> m].

Definition at line 193 of file MOM_ice_shelf_dynamics.F90.

194  real, dimension(4), intent(in) :: X !< The x-positions of the vertices of the quadrilateral [L ~> m].
195  real, dimension(4), intent(in) :: Y !< The y-positions of the vertices of the quadrilateral [L ~> m].
196  real :: quad_area ! Computed area [L2 ~> m2]
197  real :: p2, q2, a2, c2, b2, d2
198 
199 ! X and Y must be passed in the form
200  ! 3 - 4
201  ! | |
202  ! 1 - 2
203 
204  p2 = (x(4)-x(1))**2 + (y(4)-y(1))**2 ; q2 = (x(3)-x(2))**2 + (y(3)-y(2))**2
205  a2 = (x(3)-x(4))**2 + (y(3)-y(4))**2 ; c2 = (x(1)-x(2))**2 + (y(1)-y(2))**2
206  b2 = (x(2)-x(4))**2 + (y(2)-y(4))**2 ; d2 = (x(3)-x(1))**2 + (y(3)-y(1))**2
207  quad_area = .25 * sqrt(4*p2*q2-(b2+d2-a2-c2)**2)
208 

◆ register_ice_shelf_dyn_restarts()

subroutine, public mom_ice_shelf_dynamics::register_ice_shelf_dyn_restarts ( type(ocean_grid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(ice_shelf_dyn_cs), pointer  CS,
type(mom_restart_cs), pointer  restart_CS 
)

This subroutine is used to register any fields related to the ice shelf dynamics that should be written to or read from the restart file.

Parameters
[in,out]gThe grid type describing the ice shelf grid.
[in]param_fileA structure to parse for run-time parameters
csA pointer to the ice shelf dynamics control structure
restart_csA pointer to the restart control structure.

Definition at line 213 of file MOM_ice_shelf_dynamics.F90.

214  type(ocean_grid_type), intent(inout) :: G !< The grid type describing the ice shelf grid.
215  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
216  type(ice_shelf_dyn_CS), pointer :: CS !< A pointer to the ice shelf dynamics control structure
217  type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure.
218 
219  logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
220  character(len=40) :: mdl = "MOM_ice_shelf_dyn" ! This module's name.
221  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
222 
223  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
224  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
225 
226  if (associated(cs)) then
227  call mom_error(fatal, "MOM_ice_shelf_dyn.F90, register_ice_shelf_dyn_restarts: "// &
228  "called with an associated control structure.")
229  return
230  endif
231  allocate(cs)
232 
233  override_shelf_movement = .false. ; active_shelf_dynamics = .false.
234  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
235  "If true, the ice sheet mass can evolve with time.", &
236  default=.false., do_not_log=.true.)
237  if (shelf_mass_is_dynamic) then
238  call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", override_shelf_movement, &
239  "If true, user provided code specifies the ice-shelf "//&
240  "movement instead of the dynamic ice model.", default=.false., do_not_log=.true.)
241  active_shelf_dynamics = .not.override_shelf_movement
242  endif
243 
244  if (active_shelf_dynamics) then
245  allocate( cs%u_shelf(isdb:iedb,jsdb:jedb) ) ; cs%u_shelf(:,:) = 0.0
246  allocate( cs%v_shelf(isdb:iedb,jsdb:jedb) ) ; cs%v_shelf(:,:) = 0.0
247  allocate( cs%t_shelf(isd:ied,jsd:jed) ) ; cs%t_shelf(:,:) = -10.0
248  allocate( cs%ice_visc(isd:ied,jsd:jed) ) ; cs%ice_visc(:,:) = 0.0
249  allocate( cs%basal_traction(isd:ied,jsd:jed) ) ; cs%basal_traction(:,:) = 0.0
250  allocate( cs%OD_av(isd:ied,jsd:jed) ) ; cs%OD_av(:,:) = 0.0
251  allocate( cs%ground_frac(isd:ied,jsd:jed) ) ; cs%ground_frac(:,:) = 0.0
252 
253  ! additional restarts for ice shelf state
254  call register_restart_field(cs%u_shelf, "u_shelf", .false., restart_cs, &
255  "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu')
256  call register_restart_field(cs%v_shelf, "v_shelf", .false., restart_cs, &
257  "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu')
258  call register_restart_field(cs%t_shelf, "t_shelf", .true., restart_cs, &
259  "ice sheet/shelf vertically averaged temperature", "deg C")
260  call register_restart_field(cs%OD_av, "OD_av", .true., restart_cs, &
261  "Average open ocean depth in a cell","m")
262  call register_restart_field(cs%ground_frac, "ground_frac", .true., restart_cs, &
263  "fractional degree of grounding", "nondim")
264  call register_restart_field(cs%ice_visc, "viscosity", .true., restart_cs, &
265  "Volume integrated Glens law ice viscosity", "kg m2 s-1")
266  call register_restart_field(cs%basal_traction, "tau_b_beta", .true., restart_cs, &
267  "The area integrated basal traction coefficient", "kg s-1")
268  endif
269 

◆ shelf_advance_front()

subroutine, public mom_ice_shelf_dynamics::shelf_advance_front ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ice_shelf_state), intent(inout)  ISS,
type(ocean_grid_type), intent(in)  G,
real, dimension(szdi_(g),szdj_(g)), intent(inout)  hmask,
real, dimension(szdib_(g),szdj_(g)), intent(inout)  uh_ice,
real, dimension(szdi_(g),szdjb_(g)), intent(inout)  vh_ice 
)
Parameters
[in]csA pointer to the ice shelf control structure
[in,out]issA structure with elements that describe the ice-shelf state
[in]gThe grid structure used by the ice shelf.
[in,out]hmaskA mask indicating which tracer points are
[in,out]uh_iceThe accumulated zonal ice volume flux [Z L2 ~> m3]
[in,out]vh_iceThe accumulated meridional ice volume flux [Z L2 ~> m3]

Definition at line 1462 of file MOM_ice_shelf_dynamics.F90.

1463  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1464  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
1465  !! the ice-shelf state
1466  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1467  real, dimension(SZDI_(G),SZDJ_(G)), &
1468  intent(inout) :: hmask !< A mask indicating which tracer points are
1469  !! partly or fully covered by an ice-shelf
1470  real, dimension(SZDIB_(G),SZDJ_(G)), &
1471  intent(inout) :: uh_ice !< The accumulated zonal ice volume flux [Z L2 ~> m3]
1472  real, dimension(SZDI_(G),SZDJB_(G)), &
1473  intent(inout) :: vh_ice !< The accumulated meridional ice volume flux [Z L2 ~> m3]
1474 
1475  ! in this subroutine we go through the computational cells only and, if they are empty or partial cells,
1476  ! we find the reference thickness and update the shelf mass and partial area fraction and the hmask if necessary
1477 
1478  ! if any cells go from partial to complete, we then must set the thickness, update hmask accordingly,
1479  ! and divide the overflow across the adjacent EMPTY (not partly-covered) cells.
1480  ! (it is highly unlikely there will not be any; in which case this will need to be rethought.)
1481 
1482  ! most likely there will only be one "overflow". If not, though, a pass_var of all relevant variables
1483  ! is done; there will therefore be a loop which, in practice, will hopefully not have to go through
1484  ! many iterations
1485 
1486  ! when 3d advected scalars are introduced, they will be impacted by what is done here
1487 
1488  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
1489  !
1490  ! from eastern neighbor: flux_enter(:,:,1)
1491  ! from western neighbor: flux_enter(:,:,2)
1492  ! from southern neighbor: flux_enter(:,:,3)
1493  ! from northern neighbor: flux_enter(:,:,4)
1494  !
1495  ! o--- (4) ---o
1496  ! | |
1497  ! (1) (2)
1498  ! | |
1499  ! o--- (3) ---o
1500  !
1501 
1502  integer :: i, j, isc, iec, jsc, jec, n_flux, k, l, iter_count
1503  integer :: i_off, j_off
1504  integer :: iter_flag
1505 
1506  real :: h_reference ! A reference thicknesss based on neighboring cells [Z ~> m]
1507  real :: tot_flux ! The total ice mass flux [Z L2 ~> m3]
1508  real :: partial_vol ! The volume covered by ice shelf [Z L2 ~> m3]
1509  real :: dxdyh ! Cell area [L2 ~> m2]
1510  character(len=160) :: mesg ! The text of an error message
1511  integer, dimension(4) :: mapi, mapj, new_partial
1512  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter ! The ice volume flux into the
1513  ! cell through the 4 cell boundaries [Z L2 ~> m3].
1514  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter_replace ! An updated ice volume flux into the
1515  ! cell through the 4 cell boundaries [Z L2 ~> m3].
1516 
1517  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1518  i_off = g%idg_offset ; j_off = g%jdg_offset
1519  iter_count = 0 ; iter_flag = 1
1520 
1521  flux_enter(:,:,:) = 0.0
1522  do j=jsc-1,jec+1 ; do i=isc-1,iec+1
1523  if ((hmask(i,j) == 0) .or. (hmask(i,j) == 2)) then
1524  flux_enter(i,j,1) = max(uh_ice(i-1,j), 0.0)
1525  flux_enter(i,j,2) = max(-uh_ice(i,j), 0.0)
1526  flux_enter(i,j,3) = max(vh_ice(i,j-1), 0.0)
1527  flux_enter(i,j,4) = max(-vh_ice(i,j), 0.0)
1528  endif
1529  enddo ; enddo
1530 
1531  mapi(1) = -1 ; mapi(2) = 1 ; mapi(3:4) = 0
1532  mapj(3) = -1 ; mapj(4) = 1 ; mapj(1:2) = 0
1533 
1534  do while (iter_flag == 1)
1535 
1536  iter_flag = 0
1537 
1538  if (iter_count > 0) then
1539  flux_enter(:,:,:) = flux_enter_replace(:,:,:)
1540  endif
1541  flux_enter_replace(:,:,:) = 0.0
1542 
1543  iter_count = iter_count + 1
1544 
1545  ! if iter_count >= 3 then some halo updates need to be done...
1546 
1547  do j=jsc-1,jec+1
1548 
1549  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1550  ((j+j_off) >= g%domain%njhalo+1)) then
1551 
1552  do i=isc-1,iec+1
1553 
1554  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1555  ((i+i_off) >= g%domain%nihalo+1)) then
1556  ! first get reference thickness by averaging over cells that are fluxing into this cell
1557  n_flux = 0
1558  h_reference = 0.0
1559  tot_flux = 0.0
1560 
1561  do k=1,2
1562  if (flux_enter(i,j,k) > 0) then
1563  n_flux = n_flux + 1
1564  h_reference = h_reference + iss%h_shelf(i+2*k-3,j)
1565  tot_flux = tot_flux + flux_enter(i,j,k)
1566  flux_enter(i,j,k) = 0.0
1567  endif
1568  enddo
1569 
1570  do k=1,2
1571  if (flux_enter(i,j,k+2) > 0) then
1572  n_flux = n_flux + 1
1573  h_reference = h_reference + iss%h_shelf(i,j+2*k-3)
1574  tot_flux = tot_flux + flux_enter(i,j,k+2)
1575  flux_enter(i,j,k+2) = 0.0
1576  endif
1577  enddo
1578 
1579  if (n_flux > 0) then
1580  dxdyh = g%areaT(i,j)
1581  h_reference = h_reference / real(n_flux)
1582  partial_vol = iss%h_shelf(i,j) * iss%area_shelf_h(i,j) + tot_flux
1583 
1584  if ((partial_vol / g%areaT(i,j)) == h_reference) then ! cell is exactly covered, no overflow
1585  iss%hmask(i,j) = 1
1586  iss%h_shelf(i,j) = h_reference
1587  iss%area_shelf_h(i,j) = g%areaT(i,j)
1588  elseif ((partial_vol / g%areaT(i,j)) < h_reference) then
1589  iss%hmask(i,j) = 2
1590  ! ISS%mass_shelf(i,j) = partial_vol * CS%density_ice
1591  iss%area_shelf_h(i,j) = partial_vol / h_reference
1592  iss%h_shelf(i,j) = h_reference
1593  else
1594 
1595  iss%hmask(i,j) = 1
1596  iss%area_shelf_h(i,j) = g%areaT(i,j)
1597  !h_temp(i,j) = h_reference
1598  partial_vol = partial_vol - h_reference * g%areaT(i,j)
1599 
1600  iter_flag = 1
1601 
1602  n_flux = 0 ; new_partial(:) = 0
1603 
1604  do k=1,2
1605  if (cs%u_face_mask(i-2+k,j) == 2) then
1606  n_flux = n_flux + 1
1607  elseif (iss%hmask(i+2*k-3,j) == 0) then
1608  n_flux = n_flux + 1
1609  new_partial(k) = 1
1610  endif
1611  if (cs%v_face_mask(i,j-2+k) == 2) then
1612  n_flux = n_flux + 1
1613  elseif (iss%hmask(i,j+2*k-3) == 0) then
1614  n_flux = n_flux + 1
1615  new_partial(k+2) = 1
1616  endif
1617  enddo
1618 
1619  if (n_flux == 0) then ! there is nowhere to put the extra ice!
1620  iss%h_shelf(i,j) = h_reference + partial_vol / g%areaT(i,j)
1621  else
1622  iss%h_shelf(i,j) = h_reference
1623 
1624  do k=1,2
1625  if (new_partial(k) == 1) &
1626  flux_enter_replace(i+2*k-3,j,3-k) = partial_vol / real(n_flux)
1627  if (new_partial(k+2) == 1) &
1628  flux_enter_replace(i,j+2*k-3,5-k) = partial_vol / real(n_flux)
1629  enddo
1630  endif
1631 
1632  endif ! Parital_vol test.
1633  endif ! n_flux gt 0 test.
1634 
1635  endif
1636  enddo ! j-loop
1637  endif
1638  enddo
1639 
1640  ! call max_across_PEs(iter_flag)
1641 
1642  enddo ! End of do while(iter_flag) loop
1643 
1644  call max_across_pes(iter_count)
1645 
1646  if (is_root_pe() .and. (iter_count > 1)) then
1647  write(mesg,*) "shelf_advance_front: ", iter_count, " max iterations"
1648  call mom_mesg(mesg, 5)
1649  endif
1650 

◆ slope_limiter()

real function mom_ice_shelf_dynamics::slope_limiter ( real, intent(in)  num,
real, intent(in)  denom 
)
private

used for flux limiting in advective subroutines Van Leer limiter (source: Wikipedia) The return value is between 0 and 2 [nondim].

Parameters
[in]numThe numerator of the ratio used in the Van Leer slope limiter
[in]denomThe denominator of the ratio used in the Van Leer slope limiter

Definition at line 175 of file MOM_ice_shelf_dynamics.F90.

176  real, intent(in) :: num !< The numerator of the ratio used in the Van Leer slope limiter
177  real, intent(in) :: denom !< The denominator of the ratio used in the Van Leer slope limiter
178  real :: slope_limiter ! The slope limiter value, between 0 and 2 [nondim].
179  real :: r ! The ratio of num/denom [nondim]
180 
181  if (denom == 0) then
182  slope_limiter = 0
183  elseif (num*denom <= 0) then
184  slope_limiter = 0
185  else
186  r = num/denom
187  slope_limiter = (r+abs(r))/(1+abs(r))
188  endif
189 

◆ update_ice_shelf()

subroutine, public mom_ice_shelf_dynamics::update_ice_shelf ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ice_shelf_state), intent(inout)  ISS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, intent(in)  time_step,
type(time_type), intent(in)  Time,
real, dimension(szdi_(g),szdj_(g)), intent(in), optional  ocean_mass,
logical, intent(in), optional  coupled_grounding,
logical, intent(in), optional  must_update_vel 
)

This subroutine updates the ice shelf velocities, mass, stresses and properties due to the ice shelf dynamics.

Parameters
[in,out]csThe ice shelf dynamics control structure
[in,out]issA structure with elements that describe the ice-shelf state
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]time_steptime step [T ~> s]
[in]timeThe current model time
[in]ocean_massIf present this is the mass per unit area
[in]coupled_groundingIf true, the grounding line is determined by coupled ice-ocean dynamics
[in]must_update_velAlways update the ice velocities if true.

Definition at line 631 of file MOM_ice_shelf_dynamics.F90.

632  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
633  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
634  !! the ice-shelf state
635  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
636  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
637  real, intent(in) :: time_step !< time step [T ~> s]
638  type(time_type), intent(in) :: Time !< The current model time
639  real, dimension(SZDI_(G),SZDJ_(G)), &
640  optional, intent(in) :: ocean_mass !< If present this is the mass per unit area
641  !! of the ocean [R Z ~> kg m-2].
642  logical, optional, intent(in) :: coupled_grounding !< If true, the grounding line is
643  !! determined by coupled ice-ocean dynamics
644  logical, optional, intent(in) :: must_update_vel !< Always update the ice velocities if true.
645 
646  integer :: iters
647  logical :: update_ice_vel, coupled_GL
648 
649  update_ice_vel = .false.
650  if (present(must_update_vel)) update_ice_vel = must_update_vel
651 
652  coupled_gl = .false.
653  if (present(ocean_mass) .and. present(coupled_grounding)) coupled_gl = coupled_grounding
654 
655  call ice_shelf_advect(cs, iss, g, time_step, time)
656  cs%elapsed_velocity_time = cs%elapsed_velocity_time + time_step
657  if (cs%elapsed_velocity_time >= cs%velocity_update_time_step) update_ice_vel = .true.
658 
659  if (coupled_gl) then
660  call update_od_ffrac(cs, g, us, ocean_mass, update_ice_vel)
661  elseif (update_ice_vel) then
662  call update_od_ffrac_uncoupled(cs, g, iss%h_shelf(:,:))
663  endif
664 
665  if (update_ice_vel) then
666  call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, time)
667  endif
668 
669  call ice_shelf_temp(cs, iss, g, us, time_step, iss%water_flux, time)
670 
671  if (update_ice_vel) then
672  call enable_averages(cs%elapsed_velocity_time, time, cs%diag)
673  if (cs%id_col_thick > 0) call post_data(cs%id_col_thick, cs%OD_av, cs%diag)
674  if (cs%id_u_shelf > 0) call post_data(cs%id_u_shelf, cs%u_shelf, cs%diag)
675  if (cs%id_v_shelf > 0) call post_data(cs%id_v_shelf, cs%v_shelf, cs%diag)
676  if (cs%id_t_shelf > 0) call post_data(cs%id_t_shelf,cs%t_shelf,cs%diag)
677  if (cs%id_ground_frac > 0) call post_data(cs%id_ground_frac, cs%ground_frac,cs%diag)
678  if (cs%id_OD_av >0) call post_data(cs%id_OD_av, cs%OD_av,cs%diag)
679 
680  if (cs%id_u_mask > 0) call post_data(cs%id_u_mask,cs%umask,cs%diag)
681  if (cs%id_v_mask > 0) call post_data(cs%id_v_mask,cs%vmask,cs%diag)
682  if (cs%id_t_mask > 0) call post_data(cs%id_t_mask,cs%tmask,cs%diag)
683 
684  call disable_averaging(cs%diag)
685 
686  cs%elapsed_velocity_time = 0.0
687  endif
688 

◆ update_od_ffrac()

subroutine mom_ice_shelf_dynamics::update_od_ffrac ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
real, dimension(szdi_(g),szdj_(g)), intent(in)  ocean_mass,
logical, intent(in)  find_avg 
)
private
Parameters
[in,out]csA pointer to the ice shelf control structure
[in,out]gThe grid structure used by the ice shelf.
[in]usA structure containing unit conversion factors
[in]ocean_massThe mass per unit area of the ocean [kg m-2].
[in]find_avgIf true, find the average of OD and ffrac, and reset the underlying running sums to 0.

Definition at line 2510 of file MOM_ice_shelf_dynamics.F90.

2511  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
2512  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
2513  type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
2514  real, dimension(SZDI_(G),SZDJ_(G)), &
2515  intent(in) :: ocean_mass !< The mass per unit area of the ocean [kg m-2].
2516  logical, intent(in) :: find_avg !< If true, find the average of OD and ffrac, and
2517  !! reset the underlying running sums to 0.
2518 
2519  integer :: isc, iec, jsc, jec, i, j
2520  real :: I_rho_ocean ! A typical specific volume of the ocean [R-1 ~> m3 kg-1]
2521  real :: I_counter
2522 
2523  i_rho_ocean = 1.0 / cs%density_ocean_avg
2524 
2525  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2526 
2527  do j=jsc,jec ; do i=isc,iec
2528  cs%OD_rt(i,j) = cs%OD_rt(i,j) + ocean_mass(i,j)*i_rho_ocean
2529  if (ocean_mass(i,j)*i_rho_ocean > cs%thresh_float_col_depth) then
2530  cs%ground_frac_rt(i,j) = cs%ground_frac_rt(i,j) + 1.0
2531  endif
2532  enddo ; enddo
2533  cs%OD_rt_counter = cs%OD_rt_counter + 1
2534 
2535  if (find_avg) then
2536  i_counter = 1.0 / real(cs%OD_rt_counter)
2537  do j=jsc,jec ; do i=isc,iec
2538  cs%ground_frac(i,j) = 1.0 - (cs%ground_frac_rt(i,j) * i_counter)
2539  cs%OD_av(i,j) = cs%OD_rt(i,j) * i_counter
2540 
2541  cs%OD_rt(i,j) = 0.0 ; cs%ground_frac_rt(i,j) = 0.0
2542  enddo ; enddo
2543 
2544  call pass_var(cs%ground_frac, g%domain)
2545  call pass_var(cs%OD_av, g%domain)
2546  endif
2547 

◆ update_od_ffrac_uncoupled()

subroutine mom_ice_shelf_dynamics::update_od_ffrac_uncoupled ( type(ice_shelf_dyn_cs), intent(inout)  CS,
type(ocean_grid_type), intent(in)  G,
real, dimension(szdi_(g),szdj_(g)), intent(in)  h_shelf 
)
private
Parameters
[in,out]csA pointer to the ice shelf control structure
[in]gThe grid structure used by the ice shelf.
[in]h_shelfthe thickness of the ice shelf [Z ~> m].

Definition at line 2550 of file MOM_ice_shelf_dynamics.F90.

2551  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
2552  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2553  real, dimension(SZDI_(G),SZDJ_(G)), &
2554  intent(in) :: h_shelf !< the thickness of the ice shelf [Z ~> m].
2555 
2556  integer :: i, j, iters, isd, ied, jsd, jed
2557  real :: rhoi_rhow, OD
2558 
2559  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
2560  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2561 
2562  do j=jsd,jed
2563  do i=isd,ied
2564  od = g%bathyT(i,j) - rhoi_rhow * h_shelf(i,j)
2565  if (od >= 0) then
2566  ! ice thickness does not take up whole ocean column -> floating
2567  cs%OD_av(i,j) = od
2568  cs%ground_frac(i,j) = 0.
2569  else
2570  cs%OD_av(i,j) = 0.
2571  cs%ground_frac(i,j) = 1.
2572  endif
2573  enddo
2574  enddo
2575 

◆ update_velocity_masks()

subroutine mom_ice_shelf_dynamics::update_velocity_masks ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(inout)  G,
real, dimension(szdi_(g),szdj_(g)), intent(in)  hmask,
real, dimension(szdib_(g),szdjb_(g)), intent(out)  umask,
real, dimension(szdib_(g),szdjb_(g)), intent(out)  vmask,
real, dimension(szdib_(g),szdj_(g)), intent(out)  u_face_mask,
real, dimension(szdi_(g),szdjb_(g)), intent(out)  v_face_mask 
)
private
Parameters
[in]csA pointer to the ice shelf dynamics control structure
[in,out]gThe grid structure used by the ice shelf.
[in]hmaskA mask indicating which tracer points are
[out]umaskA coded mask indicating the nature of the
[out]vmaskA coded mask indicating the nature of the
[out]u_face_maskA coded mask for velocities at the C-grid u-face
[out]v_face_maskA coded mask for velocities at the C-grid v-face

Definition at line 2754 of file MOM_ice_shelf_dynamics.F90.

2755  type(ice_shelf_dyn_CS),intent(in) :: CS !< A pointer to the ice shelf dynamics control structure
2756  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
2757  real, dimension(SZDI_(G),SZDJ_(G)), &
2758  intent(in) :: hmask !< A mask indicating which tracer points are
2759  !! partly or fully covered by an ice-shelf
2760  real, dimension(SZDIB_(G),SZDJB_(G)), &
2761  intent(out) :: umask !< A coded mask indicating the nature of the
2762  !! zonal flow at the corner point
2763  real, dimension(SZDIB_(G),SZDJB_(G)), &
2764  intent(out) :: vmask !< A coded mask indicating the nature of the
2765  !! meridional flow at the corner point
2766  real, dimension(SZDIB_(G),SZDJ_(G)), &
2767  intent(out) :: u_face_mask !< A coded mask for velocities at the C-grid u-face
2768  real, dimension(SZDI_(G),SZDJB_(G)), &
2769  intent(out) :: v_face_mask !< A coded mask for velocities at the C-grid v-face
2770  ! sets masks for velocity solve
2771  ! ignores the fact that their might be ice-free cells - this only considers the computational boundary
2772 
2773  ! !!!IMPORTANT!!! relies on thickness mask - assumed that this is called after hmask has been updated & halo-updated
2774 
2775  integer :: i, j, k, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
2776  integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec
2777  integer :: i_off, j_off
2778 
2779  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2780  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
2781  i_off = g%idg_offset ; j_off = g%jdg_offset
2782  isd = g%isd ; jsd = g%jsd
2783  iegq = g%iegB ; jegq = g%jegB
2784  gisc = g%Domain%nihalo ; gjsc = g%Domain%njhalo
2785  giec = g%Domain%niglobal+gisc ; gjec = g%Domain%njglobal+gjsc
2786 
2787  umask(:,:) = 0 ; vmask(:,:) = 0
2788  u_face_mask(:,:) = 0 ; v_face_mask(:,:) = 0
2789 
2790  if (g%symmetric) then
2791  is = isd ; js = jsd
2792  else
2793  is = isd+1 ; js = jsd+1
2794  endif
2795 
2796  do j=js,g%jed
2797  do i=is,g%ied
2798 
2799  if (hmask(i,j) == 1) then
2800 
2801  umask(i-1:i,j-1:j) = 1.
2802  vmask(i-1:i,j-1:j) = 1.
2803 
2804  do k=0,1
2805 
2806  select case (int(cs%u_face_mask_bdry(i-1+k,j)))
2807  case (3)
2808  umask(i-1+k,j-1:j)=3.
2809  vmask(i-1+k,j-1:j)=0.
2810  u_face_mask(i-1+k,j)=3.
2811  case (2)
2812  u_face_mask(i-1+k,j)=2.
2813  case (4)
2814  umask(i-1+k,j-1:j)=0.
2815  vmask(i-1+k,j-1:j)=0.
2816  u_face_mask(i-1+k,j)=4.
2817  case (0)
2818  umask(i-1+k,j-1:j)=0.
2819  vmask(i-1+k,j-1:j)=0.
2820  u_face_mask(i-1+k,j)=0.
2821  case (1) ! stress free x-boundary
2822  umask(i-1+k,j-1:j)=0.
2823  case default
2824  end select
2825  enddo
2826 
2827  do k=0,1
2828 
2829  select case (int(cs%v_face_mask_bdry(i,j-1+k)))
2830  case (3)
2831  vmask(i-1:i,j-1+k)=3.
2832  umask(i-1:i,j-1+k)=0.
2833  v_face_mask(i,j-1+k)=3.
2834  case (2)
2835  v_face_mask(i,j-1+k)=2.
2836  case (4)
2837  umask(i-1:i,j-1+k)=0.
2838  vmask(i-1:i,j-1+k)=0.
2839  v_face_mask(i,j-1+k)=4.
2840  case (0)
2841  umask(i-1:i,j-1+k)=0.
2842  vmask(i-1:i,j-1+k)=0.
2843  v_face_mask(i,j-1+k)=0.
2844  case (1) ! stress free y-boundary
2845  vmask(i-1:i,j-1+k)=0.
2846  case default
2847  end select
2848  enddo
2849 
2850  !if (CS%u_face_mask_bdry(I-1,j) >= 0) then ! Western boundary
2851  ! u_face_mask(I-1,j) = CS%u_face_mask_bdry(I-1,j)
2852  ! umask(I-1,J-1:J) = 3.
2853  ! vmask(I-1,J-1:J) = 0.
2854  !endif
2855 
2856  !if (j_off+j == gjsc+1) then ! SoutherN boundary
2857  ! v_face_mask(i,J-1) = 0.
2858  ! umask(I-1:I,J-1) = 0.
2859  ! vmask(I-1:I,J-1) = 0.
2860  !elseif (j_off+j == gjec) then ! Northern boundary
2861  ! v_face_mask(i,J) = 0.
2862  ! umask(I-1:I,J) = 0.
2863  ! vmask(I-1:I,J) = 0.
2864  !endif
2865 
2866  if (i < g%ied) then
2867  if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then
2868  ! east boundary or adjacent to unfilled cell
2869  u_face_mask(i,j) = 2.
2870  endif
2871  endif
2872 
2873  if (i > g%isd) then
2874  if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
2875  !adjacent to unfilled cell
2876  u_face_mask(i-1,j) = 2.
2877  endif
2878  endif
2879 
2880  if (j > g%jsd) then
2881  if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
2882  !adjacent to unfilled cell
2883  v_face_mask(i,j-1) = 2.
2884  endif
2885  endif
2886 
2887  if (j < g%jed) then
2888  if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
2889  !adjacent to unfilled cell
2890  v_face_mask(i,j) = 2.
2891  endif
2892  endif
2893 
2894 
2895  endif
2896 
2897  enddo
2898  enddo
2899 
2900  ! note: if the grid is nonsymmetric, there is a part that will not be transferred with a halo update
2901  ! so this subroutine must update its own symmetric part of the halo
2902 
2903  call pass_vector(u_face_mask, v_face_mask, g%domain, to_all, cgrid_ne)
2904  call pass_vector(umask, vmask, g%domain, to_all, bgrid_ne)
2905 
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of a distributed 2d or 3d field.
Definition: MOM_coms.F90:53