MOM6
Data Types
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...
 

Functions/Subroutines

real function slope_limiter (num, denom)
 used for flux limiting in advective subroutines Van Leer limiter (source: Wikipedia) 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 timestep that can be taken based on the current ice velocities. Because it involves finding a global minimum, it can be suprisingly 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, v, iters, time)
 
subroutine ice_shelf_solve_inner (CS, ISS, G, u, v, taudx, taudy, H_node, float_cond, hmask, conv_flag, iters, time, Phi, Phisub)
 
subroutine ice_shelf_advect_thickness_x (CS, G, time_step, hmask, h0, h_after_uflux, flux_enter)
 
subroutine ice_shelf_advect_thickness_y (CS, G, time_step, hmask, h_after_uflux, h_after_vflux, flux_enter)
 
subroutine, public shelf_advance_front (CS, ISS, G, flux_enter)
 
subroutine, public ice_shelf_min_thickness_calve (G, h_shelf, area_shelf_h, hmask, thickness_calve)
 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, TAUD_X, TAUD_Y, OD)
 
subroutine init_boundary_values (CS, G, time, hmask, input_flux, input_thick, new_sim)
 
subroutine cg_action (uret, vret, u, v, Phi, Phisub, umask, vmask, hmask, H_node, nu, float_cond, bathyT, beta, dxdyh, G, is, ie, js, je, dens_ratio)
 
subroutine cg_action_subgrid_basal (Phisub, H, U, V, DXDYH, bathyT, dens_ratio, Ucontr, Vcontr)
 
subroutine matrix_diagonal (CS, G, float_cond, H_node, nu, beta, 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, DXDYH, bathyT, dens_ratio, Ucontr, Vcontr)
 
subroutine apply_boundary_values (CS, ISS, G, time, Phisub, H_node, nu, beta, float_cond, dens_ratio, u_bdry_contr, v_bdry_contr)
 
subroutine calc_shelf_visc (CS, ISS, G, US, u, v)
 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_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, flux_enter)
 
subroutine ice_shelf_advect_temp_y (CS, G, time_step, hmask, h_after_uflux, h_after_vflux, flux_enter)
 

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(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)  nu,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(in)  beta,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  float_cond,
real, intent(in)  dens_ratio,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(inout)  u_bdry_contr,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), 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]timeThe current model time
[in]phisubQuadrature structure weights at subgridscale
[in]h_nodeThe ice shelf thickness at nodal
[in]nuA field related to the ice viscosity from Glen's
[in]betaA 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_contrContributions to the zonal ice
[in,out]v_bdry_contrContributions to the zonal ice

Definition at line 2809 of file MOM_ice_shelf_dynamics.F90.

2809 
2810  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
2811  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
2812  !! the ice-shelf state
2813  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2814  type(time_type), intent(in) :: Time !< The current model time
2815  real, dimension(:,:,:,:,:,:), &
2816  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2817  !! locations for finite element calculations
2818  real, dimension(SZDIB_(G),SZDJB_(G)), &
2819  intent(in) :: H_node !< The ice shelf thickness at nodal
2820  !! (corner) points [Z ~> m].
2821  real, dimension(SZDIB_(G),SZDJB_(G)), &
2822  intent(in) :: nu !< A field related to the ice viscosity from Glen's
2823  !! flow law. The exact form and units depend on the
2824  !! basal law exponent.
2825  real, dimension(SZDIB_(G),SZDJB_(G)), &
2826  intent(in) :: beta !< A field related to the nonlinear part of the
2827  !! "linearized" basal stress. The exact form and
2828  !! units depend on the basal law exponent
2829  real, dimension(SZDI_(G),SZDJ_(G)), &
2830  intent(in) :: float_cond !< An array indicating where the ice
2831  !! shelf is floating: 0 if floating, 1 if not.
2832  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2833  !! of seawater, nondimensional
2834  real, dimension(SZDIB_(G),SZDJB_(G)), &
2835  intent(inout) :: u_bdry_contr !< Contributions to the zonal ice
2836  !! velocities due to the open boundaries
2837  real, dimension(SZDIB_(G),SZDJB_(G)), &
2838  intent(inout) :: v_bdry_contr !< Contributions to the zonal ice
2839  !! velocities due to the open boundaries
2840 
2841 ! this will be a per-setup function. the boundary values of thickness and velocity
2842 ! (and possibly other variables) will be updated in this function
2843 
2844  real, dimension(8,4) :: Phi
2845  real, dimension(4) :: X, Y
2846  real, dimension(2) :: xquad
2847  integer :: i, j, isc, jsc, iec, jec, iq, jq, iphi, jphi, ilq, jlq
2848  real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, uq, vq, area, basel
2849  real, dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr
2850 
2851 
2852  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2853 
2854  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2855 
2856 ! X and Y must be passed in the form
2857  ! 3 - 4
2858  ! | |
2859  ! 1 - 2
2860 ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2861 ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
2862 
2863  do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (iss%hmask(i,j) == 1) then
2864 
2865  ! process this cell if any corners have umask set to non-dirichlet bdry.
2866  ! NOTE: vmask not considered, probably should be
2867 
2868  if ((cs%umask(i-1,j-1) == 3) .OR. (cs%umask(i,j-1) == 3) .OR. &
2869  (cs%umask(i-1,j) == 3) .OR. (cs%umask(i,j) == 3)) then
2870 
2871  dxh = g%US%L_to_m*g%dxT(i,j)
2872  dyh = g%US%L_to_m*g%dyT(i,j)
2873  dxdyh = g%US%L_to_m**2*g%areaT(i,j)
2874 
2875  x(1:2) = g%geoLonBu(i-1:i,j-1)*1000
2876  x(3:4) = g%geoLonBu(i-1:i,j)*1000
2877  y(1:2) = g%geoLatBu(i-1:i,j-1)*1000
2878  y(3:4) = g%geoLatBu(i-1:i,j)*1000
2879 
2880  call bilinear_shape_functions(x, y, phi, area)
2881 
2882  ! X and Y must be passed in the form
2883  ! 3 - 4
2884  ! | |
2885  ! 1 - 2
2886  ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2887  ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
2888 
2889  do iq=1,2 ; do jq=1,2
2890 
2891  uq = cs%u_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2892  cs%u_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
2893  cs%u_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
2894  cs%u_bdry_val(i,j) * xquad(iq) * xquad(jq)
2895 
2896  vq = cs%v_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2897  cs%v_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
2898  cs%v_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
2899  cs%v_bdry_val(i,j) * xquad(iq) * xquad(jq)
2900 
2901  ux = cs%u_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
2902  cs%u_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
2903  cs%u_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
2904  cs%u_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
2905 
2906  vx = cs%v_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
2907  cs%v_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
2908  cs%v_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
2909  cs%v_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
2910 
2911  uy = cs%u_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
2912  cs%u_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
2913  cs%u_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
2914  cs%u_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
2915 
2916  vy = cs%v_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
2917  cs%v_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
2918  cs%v_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
2919  cs%v_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
2920 
2921  do iphi=1,2 ; do jphi=1,2
2922 
2923  if (iq == iphi) then
2924  ilq = 2
2925  else
2926  ilq = 1
2927  endif
2928 
2929  if (jq == jphi) then
2930  jlq = 2
2931  else
2932  jlq = 1
2933  endif
2934 
2935  if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
2936 
2937 
2938  u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
2939  .25 * dxdyh * nu(i,j) * ( (4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2940  (uy+vx) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) )
2941 
2942  if (float_cond(i,j) == 0) then
2943  u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
2944  .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
2945  endif
2946 
2947  endif
2948 
2949  if (cs%vmask(i-2+iphi,j-2+jphi) == 1) then
2950 
2951 
2952  v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
2953  .25 * dxdyh * nu(i,j) * ( (uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2954  (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2955 
2956  if (float_cond(i,j) == 0) then
2957  v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
2958  .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
2959  endif
2960 
2961  endif
2962  enddo ; enddo
2963  enddo ; enddo
2964 
2965  if (float_cond(i,j) == 1) then
2966  usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
2967  ucell(:,:) = cs%u_bdry_val(i-1:i,j-1:j) ; vcell(:,:) = cs%v_bdry_val(i-1:i,j-1:j)
2968  hcell(:,:) = h_node(i-1:i,j-1:j)
2969  call cg_action_subgrid_basal(phisub, hcell, ucell, vcell, dxdyh, basel, &
2970  dens_ratio, usubcontr, vsubcontr)
2971  do iphi=1,2 ; do jphi = 1,2
2972  if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
2973  u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
2974  usubcontr(iphi,jphi) * beta(i,j)
2975  endif
2976  if (cs%vmask(i-2+iphi,j-2+jphi) == 1) then
2977  v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
2978  vsubcontr(iphi,jphi) * beta(i,j)
2979  endif
2980  enddo ; enddo
2981  endif
2982  endif
2983  endif ; enddo ; enddo
2984 

References bilinear_shape_functions(), and cg_action_subgrid_basal().

Referenced by ice_shelf_solve_inner(), and ice_shelf_solve_outer().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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.
[in]yThe y-positions of the vertices of the quadrilateral.
[in,out]phiThe gradients of bilinear basis elements at Gaussian quadrature points surrounding the cell verticies.
[out]areaThe quadrilateral cell area [m2].

Definition at line 3121 of file MOM_ice_shelf_dynamics.F90.

3121  real, dimension(4), intent(in) :: X !< The x-positions of the vertices of the quadrilateral.
3122  real, dimension(4), intent(in) :: Y !< The y-positions of the vertices of the quadrilateral.
3123  real, dimension(8,4), intent(inout) :: Phi !< The gradients of bilinear basis elements at Gaussian
3124  !! quadrature points surrounding the cell verticies.
3125  real, intent(out) :: area !< The quadrilateral cell area [m2].
3126 
3127 ! X and Y must be passed in the form
3128  ! 3 - 4
3129  ! | |
3130  ! 1 - 2
3131 
3132 ! this subroutine calculates the gradients of bilinear basis elements that
3133 ! that are centered at the vertices of the cell. values are calculated at
3134 ! points of gaussian quadrature. (in 1D: .5 * (1 +/- sqrt(1/3)) for [0,1])
3135 ! (ordered in same way as vertices)
3136 !
3137 ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
3138 ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
3139 ! Phi_i is equal to 1 at vertex i, and 0 at vertex k /= i, and bilinear
3140 !
3141 ! This should be a one-off; once per nonlinear solve? once per lifetime?
3142 ! ... will all cells have the same shape and dimension?
3143 
3144  real, dimension(4) :: xquad, yquad
3145  integer :: node, qpoint, xnode, xq, ynode, yq
3146  real :: a,b,c,d,e,f,xexp,yexp
3147 
3148  xquad(1:3:2) = .5 * (1-sqrt(1./3)) ; yquad(1:2) = .5 * (1-sqrt(1./3))
3149  xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3))
3150 
3151  do qpoint=1,4
3152 
3153  a = -x(1)*(1-yquad(qpoint)) + x(2)*(1-yquad(qpoint)) - x(3)*yquad(qpoint) + x(4)*yquad(qpoint) ! d(x)/d(x*)
3154  b = -y(1)*(1-yquad(qpoint)) + y(2)*(1-yquad(qpoint)) - y(3)*yquad(qpoint) + y(4)*yquad(qpoint) ! d(y)/d(x*)
3155  c = -x(1)*(1-xquad(qpoint)) - x(2)*(xquad(qpoint)) + x(3)*(1-xquad(qpoint)) + x(4)*(xquad(qpoint)) ! d(x)/d(y*)
3156  d = -y(1)*(1-xquad(qpoint)) - y(2)*(xquad(qpoint)) + y(3)*(1-xquad(qpoint)) + y(4)*(xquad(qpoint)) ! d(y)/d(y*)
3157 
3158  do node=1,4
3159 
3160  xnode = 2-mod(node,2) ; ynode = ceiling(real(node)/2)
3161 
3162  if (ynode == 1) then
3163  yexp = 1-yquad(qpoint)
3164  else
3165  yexp = yquad(qpoint)
3166  endif
3167 
3168  if (1 == xnode) then
3169  xexp = 1-xquad(qpoint)
3170  else
3171  xexp = xquad(qpoint)
3172  endif
3173 
3174  phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / (a*d-b*c)
3175  phi(2*node,qpoint) = ( -c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / (a*d-b*c)
3176 
3177  enddo
3178  enddo
3179 
3180  area = quad_area(x, y)
3181 

References quad_area().

Referenced by apply_boundary_values(), ice_shelf_solve_outer(), and matrix_diagonal().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 nubmer of subgridscale quadrature locations in each direction

Definition at line 3186 of file MOM_ice_shelf_dynamics.F90.

3186  real, dimension(nsub,nsub,2,2,2,2), &
3187  intent(inout) :: Phisub !< Quadrature structure weights at subgridscale
3188  !! locations for finite element calculations
3189  integer, intent(in) :: nsub !< The nubmer of subgridscale quadrature locations in each direction
3190 
3191  ! this subroutine is a helper for interpolation of floatation condition
3192  ! for the purposes of evaluating the terms \int (u,v) \phi_i dx dy in a cell that is
3193  ! in partial floatation
3194  ! the array Phisub contains the values of \phi_i (where i is a node of the cell)
3195  ! at quad point j
3196  ! i think this general approach may not work for nonrectangular elements...
3197  !
3198 
3199  ! Phisub(i,j,k,l,q1,q2)
3200  ! i: subgrid index in x-direction
3201  ! j: subgrid index in y-direction
3202  ! k: basis function x-index
3203  ! l: basis function y-index
3204  ! q1: quad point x-index
3205  ! q2: quad point y-index
3206 
3207  ! e.g. k=1,l=1 => node 1
3208  ! q1=2,q2=1 => quad point 2
3209 
3210  ! 3 - 4
3211  ! | |
3212  ! 1 - 2
3213 
3214  integer :: i, j, k, l, qx, qy, indx, indy
3215  real,dimension(2) :: xquad
3216  real :: x0, y0, x, y, val, fracx
3217 
3218  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
3219  fracx = 1.0/real(nsub)
3220 
3221  do j=1,nsub
3222  do i=1,nsub
3223  x0 = (i-1) * fracx ; y0 = (j-1) * fracx
3224  do qx=1,2
3225  do qy=1,2
3226  x = x0 + fracx*xquad(qx)
3227  y = y0 + fracx*xquad(qy)
3228  do k=1,2
3229  do l=1,2
3230  val = 1.0
3231  if (k == 1) then
3232  val = val * (1.0-x)
3233  else
3234  val = val * x
3235  endif
3236  if (l == 1) then
3237  val = val * (1.0-y)
3238  else
3239  val = val * y
3240  endif
3241  phisub(i,j,k,l,qx,qy) = val
3242  enddo
3243  enddo
3244  enddo
3245  enddo
3246  enddo
3247  enddo
3248 

Referenced by ice_shelf_solve_outer().

Here is the caller graph for this function:

◆ 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)  TAUD_X,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  TAUD_Y,
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]usPointer to a structure containing unit conversion factors
[in]odocean floor depth at tracer points [Z ~> m].
[in,out]taud_xX-direction driving stress at q-points
[in,out]taud_yY-direction driving stress at q-points

Definition at line 2091 of file MOM_ice_shelf_dynamics.F90.

2091  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
2092  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
2093  !! the ice-shelf state
2094  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
2095  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
2096  real, dimension(SZDI_(G),SZDJ_(G)), &
2097  intent(in) :: OD !< ocean floor depth at tracer points [Z ~> m].
2098  real, dimension(SZDIB_(G),SZDJB_(G)), &
2099  intent(inout) :: TAUD_X !< X-direction driving stress at q-points
2100  real, dimension(SZDIB_(G),SZDJB_(G)), &
2101  intent(inout) :: TAUD_Y !< Y-direction driving stress at q-points
2102 
2103 ! driving stress!
2104 
2105 ! ! TAUD_X and TAUD_Y will hold driving stress in the x- and y- directions when done.
2106 ! they will sit on the BGrid, and so their size depends on whether the grid is symmetric
2107 !
2108 ! Since this is a finite element solve, they will actually have the form \int \phi_i rho g h \nabla s
2109 !
2110 ! OD -this is important and we do not yet know where (in MOM) it will come from. It represents
2111 ! "average" ocean depth -- and is needed to find surface elevation
2112 ! (it is assumed that base_ice = bed + OD)
2113 
2114  real, dimension(SIZE(OD,1),SIZE(OD,2)) :: S, & ! surface elevation [Z ~> m].
2115  BASE ! basal elevation of shelf/stream [Z ~> m].
2116 
2117 
2118  real :: rho, rhow, sx, sy, neumann_val, dxh, dyh, dxdyh, grav
2119 
2120  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
2121  integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
2122  integer :: i_off, j_off
2123 
2124  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2125  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
2126  isd = g%isd ; jsd = g%jsd
2127  iegq = g%iegB ; jegq = g%jegB
2128  gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
2129  giec = g%domain%niglobal+g%domain%nihalo ; gjec = g%domain%njglobal+g%domain%njhalo
2130  is = iscq - 1; js = jscq - 1
2131  i_off = g%idg_offset ; j_off = g%jdg_offset
2132 
2133  rho = cs%density_ice
2134  rhow = cs%density_ocean_avg
2135  grav = us%Z_to_m**2 * cs%g_Earth
2136 
2137  ! prelim - go through and calculate S
2138 
2139  ! or is this faster?
2140  base(:,:) = -g%bathyT(:,:) + od(:,:)
2141  s(:,:) = base(:,:) + iss%h_shelf(:,:)
2142 
2143  do j=jsc-1,jec+1
2144  do i=isc-1,iec+1
2145  cnt = 0
2146  sx = 0
2147  sy = 0
2148  dxh = us%L_to_m*g%dxT(i,j)
2149  dyh = us%L_to_m*g%dyT(i,j)
2150  dxdyh = us%L_to_m**2*g%areaT(i,j)
2151 
2152  if (iss%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell
2153 
2154  ! calculate sx
2155  if ((i+i_off) == gisc) then ! at left computational bdry
2156  if (iss%hmask(i+1,j) == 1) then
2157  sx = (s(i+1,j)-s(i,j))/dxh
2158  else
2159  sx = 0
2160  endif
2161  elseif ((i+i_off) == giec) then ! at right computational bdry
2162  if (iss%hmask(i-1,j) == 1) then
2163  sx = (s(i,j)-s(i-1,j))/dxh
2164  else
2165  sx = 0
2166  endif
2167  else ! interior
2168  if (iss%hmask(i+1,j) == 1) then
2169  cnt = cnt+1
2170  sx = s(i+1,j)
2171  else
2172  sx = s(i,j)
2173  endif
2174  if (iss%hmask(i-1,j) == 1) then
2175  cnt = cnt+1
2176  sx = sx - s(i-1,j)
2177  else
2178  sx = sx - s(i,j)
2179  endif
2180  if (cnt == 0) then
2181  sx = 0
2182  else
2183  sx = sx / (cnt * dxh)
2184  endif
2185  endif
2186 
2187  cnt = 0
2188 
2189  ! calculate sy, similarly
2190  if ((j+j_off) == gjsc) then ! at south computational bdry
2191  if (iss%hmask(i,j+1) == 1) then
2192  sy = (s(i,j+1)-s(i,j))/dyh
2193  else
2194  sy = 0
2195  endif
2196  elseif ((j+j_off) == gjec) then ! at nprth computational bdry
2197  if (iss%hmask(i,j-1) == 1) then
2198  sy = (s(i,j)-s(i,j-1))/dyh
2199  else
2200  sy = 0
2201  endif
2202  else ! interior
2203  if (iss%hmask(i,j+1) == 1) then
2204  cnt = cnt+1
2205  sy = s(i,j+1)
2206  else
2207  sy = s(i,j)
2208  endif
2209  if (iss%hmask(i,j-1) == 1) then
2210  cnt = cnt+1
2211  sy = sy - s(i,j-1)
2212  else
2213  sy = sy - s(i,j)
2214  endif
2215  if (cnt == 0) then
2216  sy = 0
2217  else
2218  sy = sy / (cnt * dyh)
2219  endif
2220  endif
2221 
2222  ! SW vertex
2223  taud_x(i-1,j-1) = taud_x(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2224  taud_y(i-1,j-1) = taud_y(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2225 
2226  ! SE vertex
2227  taud_x(i,j-1) = taud_x(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2228  taud_y(i,j-1) = taud_y(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2229 
2230  ! NW vertex
2231  taud_x(i-1,j) = taud_x(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2232  taud_y(i-1,j) = taud_y(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2233 
2234  ! NE vertex
2235  taud_x(i,j) = taud_x(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2236  taud_y(i,j) = taud_y(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2237 
2238  if (cs%float_frac(i,j) == 1) then
2239  neumann_val = .5 * grav * (rho * iss%h_shelf(i,j)**2 - rhow * g%bathyT(i,j)**2)
2240  else
2241  neumann_val = .5 * grav * (1-rho/rhow) * rho * iss%h_shelf(i,j)**2
2242  endif
2243 
2244 
2245  if ((cs%u_face_mask(i-1,j) == 2) .OR. (iss%hmask(i-1,j) == 0) .OR. (iss%hmask(i-1,j) == 2) ) then
2246  ! left face of the cell is at a stress boundary
2247  ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated
2248  ! pressure on either side of the face
2249  ! on the ice side, it is rho g h^2 / 2
2250  ! on the ocean side, it is rhow g (delta OD)^2 / 2
2251  ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation
2252  ! is not above the base of the ice in the current cell
2253 
2254  ! note negative sign due to direction of normal vector
2255  taud_x(i-1,j-1) = taud_x(i-1,j-1) - .5 * dyh * neumann_val
2256  taud_x(i-1,j) = taud_x(i-1,j) - .5 * dyh * neumann_val
2257  endif
2258 
2259  if ((cs%u_face_mask(i,j) == 2) .OR. (iss%hmask(i+1,j) == 0) .OR. (iss%hmask(i+1,j) == 2) ) then
2260  ! right face of the cell is at a stress boundary
2261  taud_x(i,j-1) = taud_x(i,j-1) + .5 * dyh * neumann_val
2262  taud_x(i,j) = taud_x(i,j) + .5 * dyh * neumann_val
2263  endif
2264 
2265  if ((cs%v_face_mask(i,j-1) == 2) .OR. (iss%hmask(i,j-1) == 0) .OR. (iss%hmask(i,j-1) == 2) ) then
2266  ! south face of the cell is at a stress boundary
2267  taud_y(i-1,j-1) = taud_y(i-1,j-1) - .5 * dxh * neumann_val
2268  taud_y(i,j-1) = taud_y(i,j-1) - .5 * dxh * neumann_val
2269  endif
2270 
2271  if ((cs%v_face_mask(i,j) == 2) .OR. (iss%hmask(i,j+1) == 0) .OR. (iss%hmask(i,j+1) == 2) ) then
2272  ! north face of the cell is at a stress boundary
2273  taud_y(i-1,j) = taud_y(i-1,j) + .5 * dxh * neumann_val ! note negative sign due to direction of normal vector
2274  taud_y(i,j) = taud_y(i,j) + .5 * dxh * neumann_val
2275  endif
2276 
2277  endif
2278  enddo
2279  enddo
2280 

Referenced by ice_shelf_solve_outer().

Here is the caller graph for this function:

◆ 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,
real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(inout)  v 
)
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]usPointer to a structure containing unit conversion factors
[in,out]uThe zonal ice shelf velocity [m year-1].
[in,out]vThe meridional ice shelf velocity [m year-1].

Definition at line 2990 of file MOM_ice_shelf_dynamics.F90.

2990  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
2991  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
2992  !! the ice-shelf state
2993  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2994  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
2995  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2996  intent(inout) :: u !< The zonal ice shelf velocity [m year-1].
2997  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2998  intent(inout) :: v !< The meridional ice shelf velocity [m year-1].
2999 
3000 ! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve
3001 ! so there is an "upper" and "lower" bilinear viscosity
3002 
3003 ! also this subroutine updates the nonlinear part of the basal traction
3004 
3005 ! this may be subject to change later... to make it "hybrid"
3006 
3007  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
3008  integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js
3009  real :: A, n, ux, uy, vx, vy, eps_min, umid, vmid, unorm, C_basal_friction, n_basal_friction, dxh, dyh, dxdyh
3010 
3011  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3012  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
3013  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
3014  iegq = g%iegB ; jegq = g%jegB
3015  gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
3016  giec = g%domain%niglobal+gisc ; gjec = g%domain%njglobal+gjsc
3017  is = iscq - 1; js = jscq - 1
3018 
3019  a = cs%A_glen_isothermal ; n = cs%n_glen; eps_min = cs%eps_glen_min
3020  c_basal_friction = cs%C_basal_friction ; n_basal_friction = cs%n_basal_friction
3021 
3022  do j=jsd+1,jed-1
3023  do i=isd+1,ied-1
3024 
3025  dxh = us%L_to_m*g%dxT(i,j)
3026  dyh = us%L_to_m*g%dyT(i,j)
3027  dxdyh = us%L_to_m**2*g%areaT(i,j)
3028 
3029  if (iss%hmask(i,j) == 1) then
3030  ux = (u(i,j) + u(i,j-1) - u(i-1,j) - u(i-1,j-1)) / (2*dxh)
3031  vx = (v(i,j) + v(i,j-1) - v(i-1,j) - v(i-1,j-1)) / (2*dxh)
3032  uy = (u(i,j) - u(i,j-1) + u(i-1,j) - u(i-1,j-1)) / (2*dyh)
3033  vy = (v(i,j) - v(i,j-1) + v(i-1,j) - v(i-1,j-1)) / (2*dyh)
3034 
3035  cs%ice_visc(i,j) = .5 * a**(-1/n) * &
3036  (ux**2+vy**2+ux*vy+0.25*(uy+vx)**2+eps_min**2) ** ((1-n)/(2*n)) * &
3037  us%Z_to_m*iss%h_shelf(i,j)
3038 
3039  umid = (u(i,j) + u(i,j-1) + u(i-1,j) + u(i-1,j-1))/4
3040  vmid = (v(i,j) + v(i,j-1) + v(i-1,j) + v(i-1,j-1))/4
3041  unorm = sqrt(umid**2+vmid**2+(eps_min*dxh)**2)
3042  cs%taub_beta_eff(i,j) = c_basal_friction * unorm ** (n_basal_friction-1)
3043  endif
3044  enddo
3045  enddo
3046 

Referenced by ice_shelf_solve_outer().

Here is the caller graph for this function:

◆ 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 [m2].
[in,out]hmaskA mask indicating which tracer points are
[in]calve_maskA mask that indicates where the ice shelf

Definition at line 2066 of file MOM_ice_shelf_dynamics.F90.

2066  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2067  real, dimension(SZDI_(G),SZDJ_(G)), &
2068  intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
2069  real, dimension(SZDI_(G),SZDJ_(G)), &
2070  intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [m2].
2071  real, dimension(SZDI_(G),SZDJ_(G)), &
2072  intent(inout) :: hmask !< A mask indicating which tracer points are
2073  !! partly or fully covered by an ice-shelf
2074  real, dimension(SZDI_(G),SZDJ_(G)), &
2075  intent(in) :: calve_mask !< A mask that indicates where the ice shelf
2076  !! can exist, and where it will calve.
2077 
2078  integer :: i,j
2079 
2080  do j=g%jsc,g%jec ; do i=g%isc,g%iec
2081  if ((calve_mask(i,j) == 0.0) .and. (hmask(i,j) /= 0.0)) then
2082  h_shelf(i,j) = 0.0
2083  area_shelf_h(i,j) = 0.0
2084  hmask(i,j) = 0.0
2085  endif
2086  enddo ; enddo
2087 

Referenced by ice_shelf_advect().

Here is the caller graph for this function:

◆ 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,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  v,
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)  nu,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  float_cond,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  bathyT,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(in)  beta,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  dxdyh,
type(ocean_grid_type), intent(in)  G,
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.
[in,out]vretThe retarding stresses working at v-points.
[in]phiThe gradients of bilinear basis elements at Gaussian
[in]phisubQuadrature structure weights at subgridscale
[in]uThe zonal ice shelf velocity at vertices [m year-1]
[in]vThe meridional ice shelf velocity at vertices [m year-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]nuA 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]betaA field related to the nonlinear part of the
[in]dxdyhThe tracer cell area [m2]
[in]dens_ratioThe density of ice divided by the density of seawater, nondimensional
[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 2355 of file MOM_ice_shelf_dynamics.F90.

2355 
2356  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2357  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2358  intent(inout) :: uret !< The retarding stresses working at u-points.
2359  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2360  intent(inout) :: vret !< The retarding stresses working at v-points.
2361  real, dimension(SZDI_(G),SZDJ_(G),8,4), &
2362  intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian
2363  !! quadrature points surrounding the cell verticies.
2364  real, dimension(:,:,:,:,:,:), &
2365  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2366  !! locations for finite element calculations
2367  real, dimension(SZDIB_(G),SZDJB_(G)), &
2368  intent(in) :: u !< The zonal ice shelf velocity at vertices [m year-1]
2369  real, dimension(SZDIB_(G),SZDJB_(G)), &
2370  intent(in) :: v !< The meridional ice shelf velocity at vertices [m year-1]
2371  real, dimension(SZDIB_(G),SZDJB_(G)), &
2372  intent(in) :: umask !< A coded mask indicating the nature of the
2373  !! zonal flow at the corner point
2374  real, dimension(SZDIB_(G),SZDJB_(G)), &
2375  intent(in) :: vmask !< A coded mask indicating the nature of the
2376  !! meridional flow at the corner point
2377  real, dimension(SZDIB_(G),SZDJB_(G)), &
2378  intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
2379  !! points [Z ~> m].
2380  real, dimension(SZDI_(G),SZDJ_(G)), &
2381  intent(in) :: hmask !< A mask indicating which tracer points are
2382  !! partly or fully covered by an ice-shelf
2383  real, dimension(SZDIB_(G),SZDJB_(G)), &
2384  intent(in) :: nu !< A field related to the ice viscosity from Glen's
2385  !! flow law. The exact form and units depend on the
2386  !! basal law exponent.
2387  real, dimension(SZDI_(G),SZDJ_(G)), &
2388  intent(in) :: float_cond !< An array indicating where the ice
2389  !! shelf is floating: 0 if floating, 1 if not.
2390  real, dimension(SZDI_(G),SZDJ_(G)), &
2391  intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
2392  real, dimension(SZDIB_(G),SZDJB_(G)), &
2393  intent(in) :: beta !< A field related to the nonlinear part of the
2394  !! "linearized" basal stress. The exact form and
2395  !! units depend on the basal law exponent.
2396  ! and/or whether flow is "hybridized"
2397  real, dimension(SZDI_(G),SZDJ_(G)), &
2398  intent(in) :: dxdyh !< The tracer cell area [m2]
2399  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2400  !! of seawater, nondimensional
2401  integer, intent(in) :: is !< The starting i-index to work on
2402  integer, intent(in) :: ie !< The ending i-index to work on
2403  integer, intent(in) :: js !< The starting j-index to work on
2404  integer, intent(in) :: je !< The ending j-index to work on
2405 
2406 ! the linear action of the matrix on (u,v) with bilinear finite elements
2407 ! as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
2408 ! but this may change pursuant to conversations with others
2409 !
2410 ! is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
2411 ! in order to make less frequent halo updates
2412 
2413 ! the linear action of the matrix on (u,v) with bilinear finite elements
2414 ! Phi has the form
2415 ! Phi(i,j,k,q) - applies to cell i,j
2416 
2417  ! 3 - 4
2418  ! | |
2419  ! 1 - 2
2420 
2421 ! Phi(i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
2422 ! Phi(i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
2423 ! Phi_k is equal to 1 at vertex k, and 0 at vertex l /= k, and bilinear
2424 
2425  real :: ux, vx, uy, vy, uq, vq, area, basel
2426  integer :: iq, jq, iphi, jphi, i, j, ilq, jlq
2427  real, dimension(2) :: xquad
2428  real, dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr,Ucontr
2429 
2430  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2431 
2432  do j=js,je
2433  do i=is,ie ; if (hmask(i,j) == 1) then
2434 ! dxh = G%dxh(i,j)
2435 ! dyh = G%dyh(i,j)
2436 !
2437 ! X(:,:) = G%geoLonBu(i-1:i,j-1:j)
2438 ! Y(:,:) = G%geoLatBu(i-1:i,j-1:j)
2439 !
2440 ! call bilinear_shape_functions (X, Y, Phi, area)
2441 
2442  ! X and Y must be passed in the form
2443  ! 3 - 4
2444  ! | |
2445  ! 1 - 2
2446  ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2447  ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
2448 
2449  area = dxdyh(i,j)
2450 
2451  ucontr=0
2452  do iq=1,2 ; do jq=1,2
2453 
2454 
2455  if (iq == 2) then
2456  ilq = 2
2457  else
2458  ilq = 1
2459  endif
2460 
2461  if (jq == 2) then
2462  jlq = 2
2463  else
2464  jlq = 1
2465  endif
2466 
2467  uq = u(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2468  u(i,j-1) * xquad(iq) * xquad(3-jq) + &
2469  u(i-1,j) * xquad(3-iq) * xquad(jq) + &
2470  u(i,j) * xquad(iq) * xquad(jq)
2471 
2472  vq = v(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2473  v(i,j-1) * xquad(iq) * xquad(3-jq) + &
2474  v(i-1,j) * xquad(3-iq) * xquad(jq) + &
2475  v(i,j) * xquad(iq) * xquad(jq)
2476 
2477  ux = u(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
2478  u(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
2479  u(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
2480  u(i,j) * phi(i,j,7,2*(jq-1)+iq)
2481 
2482  vx = v(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
2483  v(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
2484  v(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
2485  v(i,j) * phi(i,j,7,2*(jq-1)+iq)
2486 
2487  uy = u(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
2488  u(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
2489  u(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
2490  u(i,j) * phi(i,j,8,2*(jq-1)+iq)
2491 
2492  vy = v(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
2493  v(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
2494  v(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
2495  v(i,j) * phi(i,j,8,2*(jq-1)+iq)
2496 
2497  do iphi=1,2 ; do jphi=1,2
2498  if (umask(i-2+iphi,j-2+jphi) == 1) then
2499 
2500  uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
2501  .25 * area * nu(i,j) * ((4*ux+2*vy) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2502  (uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2503  endif
2504  if (vmask(i-2+iphi,j-2+jphi) == 1) then
2505 
2506  vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
2507  .25 * area * nu(i,j) * ((uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2508  (4*vy+2*ux) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2509  endif
2510 
2511  if (iq == iphi) then
2512  ilq = 2
2513  else
2514  ilq = 1
2515  endif
2516 
2517  if (jq == jphi) then
2518  jlq = 2
2519  else
2520  jlq = 1
2521  endif
2522 
2523  if (float_cond(i,j) == 0) then
2524 
2525  if (umask(i-2+iphi,j-2+jphi) == 1) then
2526 
2527  uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
2528  .25 * beta(i,j) * area * uq * xquad(ilq) * xquad(jlq)
2529 
2530  endif
2531 
2532  if (vmask(i-2+iphi,j-2+jphi) == 1) then
2533 
2534  vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
2535  .25 * beta(i,j) * area * vq * xquad(ilq) * xquad(jlq)
2536 
2537  endif
2538 
2539  endif
2540  ucontr(iphi,jphi) = ucontr(iphi,jphi) + .25 * area * uq * xquad(ilq) * xquad(jlq) * beta(i,j)
2541  enddo ; enddo
2542  enddo ; enddo
2543 
2544  if (float_cond(i,j) == 1) then
2545  usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = bathyt(i,j)
2546  ucell(:,:) = u(i-1:i,j-1:j) ; vcell(:,:) = v(i-1:i,j-1:j) ; hcell(:,:) = h_node(i-1:i,j-1:j)
2547  call cg_action_subgrid_basal(phisub, hcell, ucell, vcell, area, basel, &
2548  dens_ratio, usubcontr, vsubcontr)
2549  do iphi=1,2 ; do jphi=1,2
2550  if (umask(i-2+iphi,j-2+jphi) == 1) then
2551  uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
2552  endif
2553  if (vmask(i-2+iphi,j-2+jphi) == 1) then
2554  vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
2555  endif
2556  enddo ; enddo
2557  endif
2558 
2559  endif
2560  enddo ; enddo
2561 

References cg_action_subgrid_basal().

Referenced by ice_shelf_solve_inner(), and ice_shelf_solve_outer().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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)  DXDYH,
real, intent(in)  bathyT,
real, intent(in)  dens_ratio,
real, dimension(2,2), intent(inout)  Ucontr,
real, dimension(2,2), intent(inout)  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 [m year-1]
[in]vThe meridional ice shelf velocity at vertices [m year-1]
[in]dxdyhThe tracer cell area [m2]
[in]bathytThe depth of ocean bathymetry at tracer points [Z ~> m].
[in]dens_ratioThe density of ice divided by the density of seawater, nondimensional
[in,out]ucontrA field related to the subgridscale contributions to the u-direction basal stress.
[in,out]vcontrA field related to the subgridscale contributions to the v-direction basal stress.

Definition at line 2565 of file MOM_ice_shelf_dynamics.F90.

2565  real, dimension(:,:,:,:,:,:), &
2566  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2567  !! locations for finite element calculations
2568  real, dimension(2,2), intent(in) :: H !< The ice shelf thickness at nodal (corner) points [Z ~> m].
2569  real, dimension(2,2), intent(in) :: U !< The zonal ice shelf velocity at vertices [m year-1]
2570  real, dimension(2,2), intent(in) :: V !< The meridional ice shelf velocity at vertices [m year-1]
2571  real, intent(in) :: DXDYH !< The tracer cell area [m2]
2572  real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
2573  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2574  !! of seawater, nondimensional
2575  real, dimension(2,2), intent(inout) :: Ucontr !< A field related to the subgridscale contributions to
2576  !! the u-direction basal stress.
2577  real, dimension(2,2), intent(inout) :: Vcontr !< A field related to the subgridscale contributions to
2578  !! the v-direction basal stress.
2579 
2580  integer :: nsub, i, j, k, l, qx, qy, m, n
2581  real :: subarea, hloc, uq, vq
2582 
2583  nsub = size(phisub,1)
2584  subarea = dxdyh / (nsub**2)
2585 
2586  do m=1,2
2587  do n=1,2
2588  do j=1,nsub
2589  do i=1,nsub
2590  do qx=1,2
2591  do qy = 1,2
2592 
2593  hloc = phisub(i,j,1,1,qx,qy)*h(1,1) + phisub(i,j,1,2,qx,qy)*h(1,2) + &
2594  phisub(i,j,2,1,qx,qy)*h(2,1) + phisub(i,j,2,2,qx,qy)*h(2,2)
2595 
2596  if (dens_ratio * hloc - bathyt > 0) then
2597  !if (.true.) then
2598  uq = 0 ; vq = 0
2599  do k=1,2
2600  do l=1,2
2601  !Ucontr(m,n) = Ucontr(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * Phisub(i,j,k,l,qx,qy) * U(k,l)
2602  !Vcontr(m,n) = Vcontr(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * Phisub(i,j,k,l,qx,qy) * V(k,l)
2603  uq = uq + phisub(i,j,k,l,qx,qy) * u(k,l) ; vq = vq + phisub(i,j,k,l,qx,qy) * v(k,l)
2604  enddo
2605  enddo
2606 
2607  ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * uq
2608  vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * vq
2609 
2610  endif
2611 
2612  enddo
2613  enddo
2614  enddo
2615  enddo
2616  enddo
2617  enddo
2618 

Referenced by apply_boundary_values(), and cg_action().

Here is the caller graph for this function:

◆ 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)  DXDYH,
real, intent(in)  bathyT,
real, intent(in)  dens_ratio,
real, dimension(2,2), intent(inout)  Ucontr,
real, dimension(2,2), intent(inout)  Vcontr 
)
private
Parameters
[in]phisubQuadrature structure weights at subgridscale
[in]h_nodeThe ice shelf thickness at nodal (corner) points [Z ~> m].
[in]dxdyhThe tracer cell area [m2]
[in]bathytThe depth of ocean bathymetry at tracer points [Z ~> m].
[in]dens_ratioThe density of ice divided by the density of seawater, nondimensional
[in,out]ucontrA field related to the subgridscale contributions to the u-direction diagonal elements from basal stress.
[in,out]vcontrA field related to the subgridscale contributions to the v-direction diagonal elements from basal stress.

Definition at line 2770 of file MOM_ice_shelf_dynamics.F90.

2770  real, dimension(:,:,:,:,:,:), &
2771  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2772  !! locations for finite element calculations
2773  real, dimension(2,2), intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
2774  !! points [Z ~> m].
2775  real, intent(in) :: DXDYH !< The tracer cell area [m2]
2776  real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m].
2777  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2778  !! of seawater, nondimensional
2779  real, dimension(2,2), intent(inout) :: Ucontr !< A field related to the subgridscale contributions to
2780  !! the u-direction diagonal elements from basal stress.
2781  real, dimension(2,2), intent(inout) :: Vcontr !< A field related to the subgridscale contributions to
2782  !! the v-direction diagonal elements from basal stress.
2783 
2784  ! bathyT = cellwise-constant bed elevation
2785 
2786  integer :: nsub, i, j, k, l, qx, qy, m, n
2787  real :: subarea, hloc
2788 
2789  nsub = size(phisub,1)
2790  subarea = dxdyh / (nsub**2)
2791 
2792  do m=1,2 ; do n=1,2 ; do j=1,nsub ; do i=1,nsub ; do qx=1,2 ; do qy = 1,2
2793 
2794  hloc = phisub(i,j,1,1,qx,qy)*h_node(1,1) + phisub(i,j,1,2,qx,qy)*h_node(1,2) + &
2795  phisub(i,j,2,1,qx,qy)*h_node(2,1) + phisub(i,j,2,2,qx,qy)*h_node(2,2)
2796 
2797  if (dens_ratio * hloc - bathyt > 0) then
2798  ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
2799  vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
2800  endif
2801 
2802  enddo ; enddo ; enddo ; enddo ; enddo ; enddo
2803 

Referenced by matrix_diagonal().

Here is the caller graph for this function:

◆ 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 [s]
[in]timeThe current model time

Definition at line 672 of file MOM_ice_shelf_dynamics.F90.

672  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
673  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
674  !! the ice-shelf state
675  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
676  real, intent(in) :: time_step !< time step [s]
677  type(time_type), intent(in) :: Time !< The current model time
678 
679 
680 ! 3/8/11 DNG
681 !
682 ! This subroutine takes the velocity (on the Bgrid) and timesteps h_t = - div (uh) once.
683 ! ADDITIONALLY, it will update the volume of ice in partially-filled cells, and update
684 ! hmask accordingly
685 !
686 ! The flux overflows are included here. That is because they will be used to advect 3D scalars
687 ! into partial cells
688 
689  !
690  ! flux_enter: this is to capture flow into partially covered cells; it gives the mass flux into a given
691  ! cell across its boundaries.
692  ! ###Perhaps flux_enter should be changed into u-face and v-face
693  ! ###fluxes, which can then be used in halo updates, etc.
694  !
695  ! from left neighbor: flux_enter(:,:,1)
696  ! from right neighbor: flux_enter(:,:,2)
697  ! from bottom neighbor: flux_enter(:,:,3)
698  ! from top neighbor: flux_enter(:,:,4)
699  !
700  ! THESE ARE NOT CONSISTENT ==> FIND OUT WHAT YOU IMPLEMENTED
701 
702  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
703  !
704  ! o--- (4) ---o
705  ! | |
706  ! (1) (2)
707  ! | |
708  ! o--- (3) ---o
709  !
710 
711  real, dimension(SZDI_(G),SZDJ_(G)) :: h_after_uflux, h_after_vflux ! Ice thicknesses [Z ~> m].
712  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter
713  integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
714  real :: rho, spy
715 
716  rho = cs%density_ice
717  spy = 365 * 86400 ! seconds per year; is there a global constant for this? No - it is dependent upon a calendar.
718 
719  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
720  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
721  flux_enter(:,:,:) = 0.0
722 
723  h_after_uflux(:,:) = 0.0
724  h_after_vflux(:,:) = 0.0
725  ! call MOM_mesg("MOM_ice_shelf.F90: ice_shelf_advect called")
726 
727  do j=jsd,jed ; do i=isd,ied ; if (cs%thickness_bdry_val(i,j) /= 0.0) then
728  iss%h_shelf(i,j) = cs%thickness_bdry_val(i,j)
729  endif ; enddo ; enddo
730 
731  call ice_shelf_advect_thickness_x(cs, g, time_step/spy, iss%hmask, iss%h_shelf, h_after_uflux, flux_enter)
732 
733 ! call enable_averaging(time_step,Time,CS%diag)
734  ! call pass_var(h_after_uflux, G%domain)
735 ! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag)
736 ! call disable_averaging(CS%diag)
737 
738  call ice_shelf_advect_thickness_y(cs, g, time_step/spy, iss%hmask, h_after_uflux, h_after_vflux, flux_enter)
739 
740 ! call enable_averaging(time_step,Time,CS%diag)
741 ! call pass_var(h_after_vflux, G%domain)
742 ! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag)
743 ! call disable_averaging(CS%diag)
744 
745  do j=jsd,jed
746  do i=isd,ied
747  if (iss%hmask(i,j) == 1) iss%h_shelf(i,j) = h_after_vflux(i,j)
748  enddo
749  enddo
750 
751  if (cs%moving_shelf_front) then
752  call shelf_advance_front(cs, iss, g, flux_enter)
753  if (cs%min_thickness_simple_calve > 0.0) then
754  call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
755  cs%min_thickness_simple_calve)
756  endif
757  if (cs%calve_to_mask) then
758  call calve_to_mask(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, cs%calve_mask)
759  endif
760  endif
761 
762  !call enable_averaging(time_step,Time,CS%diag)
763  !if (CS%id_h_after_adv > 0) call post_data(CS%id_h_after_adv, ISS%h_shelf, CS%diag)
764  !call disable_averaging(CS%diag)
765 
766  !call change_thickness_using_melt(ISS, G, time_step, fluxes, CS%density_ice)
767 
768  call update_velocity_masks(cs, g, iss%hmask, cs%umask, cs%vmask, cs%u_face_mask, cs%v_face_mask)
769 

References calve_to_mask(), ice_shelf_advect_thickness_x(), ice_shelf_advect_thickness_y(), ice_shelf_min_thickness_calve(), shelf_advance_front(), and update_velocity_masks().

Referenced by update_ice_shelf().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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( g %isd: g %ied, g %jsd: g %jed), intent(in)  hmask,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h0,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  h_after_uflux,
real, dimension( g %isd: g %ied, g %jsd: g %jed,4), intent(inout)  flux_enter 
)
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 [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
[in,out]flux_enterThe integrated temperature flux into

Definition at line 3613 of file MOM_ice_shelf_dynamics.F90.

3613  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
3614  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
3615  real, intent(in) :: time_step !< The time step for this update [s].
3616  real, dimension(SZDI_(G),SZDJ_(G)), &
3617  intent(in) :: hmask !< A mask indicating which tracer points are
3618  !! partly or fully covered by an ice-shelf
3619  real, dimension(SZDI_(G),SZDJ_(G)), &
3620  intent(in) :: h0 !< The initial ice shelf thicknesses [Z ~> m].
3621  real, dimension(SZDI_(G),SZDJ_(G)), &
3622  intent(inout) :: h_after_uflux !< The ice shelf thicknesses after
3623  !! the zonal mass fluxes [Z ~> m].
3624  real, dimension(SZDI_(G),SZDJ_(G),4), &
3625  intent(inout) :: flux_enter !< The integrated temperature flux into
3626  !! the cell through the 4 cell boundaries [degC Z m2 ~> degC m3]
3627 
3628  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
3629 
3630  ! if there is an input bdry condition, the thickness there will be set in initialization
3631 
3632  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
3633  !
3634  ! from left neighbor: flux_enter(:,:,1)
3635  ! from right neighbor: flux_enter(:,:,2)
3636  ! from bottom neighbor: flux_enter(:,:,3)
3637  ! from top neighbor: flux_enter(:,:,4)
3638  !
3639  ! o--- (4) ---o
3640  ! | |
3641  ! (1) (2)
3642  ! | |
3643  ! o--- (3) ---o
3644  !
3645 
3646  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3647  integer :: i_off, j_off
3648  logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
3649  real, dimension(-2:2) :: stencil
3650  real :: u_face, & ! positive if out
3651  flux_diff_cell, phi, dxh, dyh, dxdyh
3652 
3653  character (len=1) :: debug_str
3654 
3655 
3656  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
3657  i_off = g%idg_offset ; j_off = g%jdg_offset
3658 
3659  do j=jsd+1,jed-1
3660  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
3661  ((j+j_off) >= g%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries
3662 
3663  stencil(:) = -1
3664 ! if (i+i_off == G%domain%nihalo+G%domain%nihalo)
3665  do i=is,ie
3666 
3667  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
3668  ((i+i_off) >= g%domain%nihalo+1)) then
3669 
3670  if (i+i_off == g%domain%nihalo+1) then
3671  at_west_bdry=.true.
3672  else
3673  at_west_bdry=.false.
3674  endif
3675 
3676  if (i+i_off == g%domain%niglobal+g%domain%nihalo) then
3677  at_east_bdry=.true.
3678  else
3679  at_east_bdry=.false.
3680  endif
3681 
3682  if (hmask(i,j) == 1) then
3683 
3684  dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
3685 
3686  h_after_uflux(i,j) = h0(i,j)
3687 
3688  stencil(:) = h0(i-2:i+2,j) ! fine as long has nx_halo >= 2
3689 
3690  flux_diff_cell = 0
3691 
3692  ! 1ST DO LEFT FACE
3693 
3694  if (cs%u_face_mask(i-1,j) == 4.) then
3695 
3696  flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i-1,j) * &
3697  cs%t_bdry_val(i-1,j) / dxdyh
3698  else
3699 
3700  ! get u-velocity at center of left face
3701  u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3702 
3703  if (u_face > 0) then !flux is into cell - we need info from h(i-2), h(i-1) if available
3704 
3705  ! i may not cover all the cases.. but i cover the realistic ones
3706 
3707  if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a
3708  ! thickness bdry condition, and the stencil contains it
3709  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
3710 
3711  elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then ! h(i-2) and h(i-1) are valid
3712  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3713  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
3714  (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3715 
3716  else ! h(i-1) is valid
3717  ! (o.w. flux would most likely be out of cell)
3718  ! but h(i-2) is not
3719 
3720  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(-1)
3721 
3722  endif
3723 
3724  elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
3725  if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid
3726  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3727  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
3728  (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3729 
3730  else
3731  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
3732 
3733  if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
3734  flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
3735  endif
3736  endif
3737  endif
3738  endif
3739 
3740  ! NEXT DO RIGHT FACE
3741 
3742  ! get u-velocity at center of right face
3743 
3744  if (cs%u_face_mask(i+1,j) == 4.) then
3745 
3746  flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i+1,j) *&
3747  cs%t_bdry_val(i+1,j)/ dxdyh
3748  else
3749 
3750  u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3751 
3752  if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available
3753 
3754  if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a
3755  ! thickness bdry condition, and the stencil contains it
3756 
3757  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
3758 
3759  elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then ! h(i+2) and h(i+1) are valid
3760 
3761  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3762  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
3763  (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3764 
3765  else ! h(i+1) is valid
3766  ! (o.w. flux would most likely be out of cell)
3767  ! but h(i+2) is not
3768 
3769  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(1)
3770 
3771  endif
3772 
3773  elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
3774 
3775  if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid
3776 
3777  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3778  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
3779  (stencil(0) - phi * (stencil(0)-stencil(1))/2)
3780 
3781  else ! h(i+1) is valid
3782  ! (o.w. flux would most likely be out of cell)
3783  ! but h(i+2) is not
3784 
3785  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
3786 
3787  if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then
3788 
3789  flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step * stencil(0)
3790  endif
3791 
3792  endif
3793 
3794  endif
3795 
3796  h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
3797 
3798  endif
3799 
3800  elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
3801 
3802  if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then
3803  u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3804  flux_enter(i,j,1) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%t_bdry_val(i-1,j)* &
3805  cs%thickness_bdry_val(i+1,j)
3806  elseif (cs%u_face_mask(i-1,j) == 4.) then
3807  flux_enter(i,j,1) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i-1,j)*cs%t_bdry_val(i-1,j)
3808  endif
3809 
3810  if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then
3811  u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3812  flux_enter(i,j,2) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%t_bdry_val(i+1,j)* &
3813  cs%thickness_bdry_val(i+1,j)
3814  elseif (cs%u_face_mask(i+1,j) == 4.) then
3815  flux_enter(i,j,2) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i+1,j) * cs%t_bdry_val(i+1,j)
3816  endif
3817 
3818 ! if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1)) then
3819  ! this is solely for the purposes of keeping the mask consistent while advancing
3820  ! the front without having to call pass_var - if cell is empty and cell to left
3821  ! is ice-covered then this cell will become partly covered
3822 ! hmask(i,j) = 2
3823 ! elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1)) then
3824  ! this is solely for the purposes of keeping the mask consistent while advancing
3825  ! the front without having to call pass_var - if cell is empty and cell to left
3826  ! is ice-covered then this cell will become partly covered
3827 ! hmask(i,j) = 2
3828 
3829 ! endif
3830 
3831  endif
3832 
3833  endif
3834 
3835  enddo ! i loop
3836 
3837  endif
3838 
3839  enddo ! j loop
3840 

References slope_limiter().

Referenced by ice_shelf_temp().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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,
real, dimension(szdi_(g),szdj_(g),4), intent(inout)  flux_enter 
)
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 [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
[in,out]flux_enterThe integrated temperature flux into

Definition at line 3844 of file MOM_ice_shelf_dynamics.F90.

3844  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
3845  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
3846  real, intent(in) :: time_step !< The time step for this update [s].
3847  real, dimension(SZDI_(G),SZDJ_(G)), &
3848  intent(in) :: hmask !< A mask indicating which tracer points are
3849  !! partly or fully covered by an ice-shelf
3850  real, dimension(SZDI_(G),SZDJ_(G)), &
3851  intent(in) :: h_after_uflux !< The ice shelf thicknesses after
3852  !! the zonal mass fluxes [Z ~> m].
3853  real, dimension(SZDI_(G),SZDJ_(G)), &
3854  intent(inout) :: h_after_vflux !< The ice shelf thicknesses after
3855  !! the meridional mass fluxes [Z ~> m].
3856  real, dimension(SZDI_(G),SZDJ_(G),4), &
3857  intent(inout) :: flux_enter !< The integrated temperature flux into
3858  !! the cell through the 4 cell boundaries [degC Z m2 ~> degC m3]
3859 
3860  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
3861 
3862  ! if there is an input bdry condition, the thickness there will be set in initialization
3863 
3864  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
3865  !
3866  ! from left neighbor: flux_enter(:,:,1)
3867  ! from right neighbor: flux_enter(:,:,2)
3868  ! from bottom neighbor: flux_enter(:,:,3)
3869  ! from top neighbor: flux_enter(:,:,4)
3870  !
3871  ! o--- (4) ---o
3872  ! | |
3873  ! (1) (2)
3874  ! | |
3875  ! o--- (3) ---o
3876  !
3877 
3878  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3879  integer :: i_off, j_off
3880  logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
3881  real, dimension(-2:2) :: stencil
3882  real :: v_face, & ! positive if out
3883  flux_diff_cell, phi, dxh, dyh, dxdyh
3884  character(len=1) :: debug_str
3885 
3886  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
3887  i_off = g%idg_offset ; j_off = g%jdg_offset
3888 
3889  do i=isd+2,ied-2
3890  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
3891  ((i+i_off) >= g%domain%nihalo+1)) then ! based on mehmet's code - only if btw east & west boundaries
3892 
3893  stencil(:) = -1
3894 
3895  do j=js,je
3896 
3897  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
3898  ((j+j_off) >= g%domain%njhalo+1)) then
3899 
3900  if (j+j_off == g%domain%njhalo+1) then
3901  at_south_bdry=.true.
3902  else
3903  at_south_bdry=.false.
3904  endif
3905  if (j+j_off == g%domain%njglobal+g%domain%njhalo) then
3906  at_north_bdry=.true.
3907  else
3908  at_north_bdry=.false.
3909  endif
3910 
3911  if (hmask(i,j) == 1) then
3912  dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
3913  h_after_vflux(i,j) = h_after_uflux(i,j)
3914 
3915  stencil(:) = h_after_uflux(i,j-2:j+2) ! fine as long has ny_halo >= 2
3916  flux_diff_cell = 0
3917 
3918  ! 1ST DO south FACE
3919 
3920  if (cs%v_face_mask(i,j-1) == 4.) then
3921 
3922  flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j-1) * &
3923  cs%t_bdry_val(i,j-1)/ dxdyh
3924  else
3925 
3926  ! get u-velocity at center of left face
3927  v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
3928 
3929  if (v_face > 0) then !flux is into cell - we need info from h(j-2), h(j-1) if available
3930 
3931  ! i may not cover all the cases.. but i cover the realistic ones
3932 
3933  if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a
3934  ! thickness bdry condition, and the stencil contains it
3935  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
3936 
3937  elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then ! h(j-2) and h(j-1) are valid
3938 
3939  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3940  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
3941  (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3942 
3943  else ! h(j-1) is valid
3944  ! (o.w. flux would most likely be out of cell)
3945  ! but h(j-2) is not
3946  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(-1)
3947  endif
3948 
3949  elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
3950 
3951  if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid
3952  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3953  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
3954  (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3955  else
3956  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
3957 
3958  if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
3959  flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
3960  endif
3961 
3962  endif
3963 
3964  endif
3965 
3966  endif
3967 
3968  ! NEXT DO north FACE
3969 
3970  if (cs%v_face_mask(i,j+1) == 4.) then
3971 
3972  flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j+1) *&
3973  cs%t_bdry_val(i,j+1)/ dxdyh
3974  else
3975 
3976  ! get u-velocity at center of right face
3977  v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
3978 
3979  if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available
3980 
3981  if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a
3982  ! thickness bdry condition, and the stencil contains it
3983  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
3984  elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then ! h(j+2) and h(j+1) are valid
3985  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3986  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
3987  (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3988  else ! h(j+1) is valid
3989  ! (o.w. flux would most likely be out of cell)
3990  ! but h(j+2) is not
3991  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
3992  endif
3993 
3994  elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
3995 
3996  if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid
3997  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3998  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
3999  (stencil(0) - phi * (stencil(0)-stencil(1))/2)
4000  else ! h(j+1) is valid
4001  ! (o.w. flux would most likely be out of cell)
4002  ! but h(j+2) is not
4003  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
4004  if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
4005  flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
4006  endif
4007  endif
4008 
4009  endif
4010 
4011  endif
4012 
4013  h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
4014 
4015  elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
4016 
4017  if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then
4018  v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
4019  flux_enter(i,j,3) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%t_bdry_val(i,j-1)* &
4020  cs%thickness_bdry_val(i,j-1)
4021  elseif (cs%v_face_mask(i,j-1) == 4.) then
4022  flux_enter(i,j,3) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j-1)*cs%t_bdry_val(i,j-1)
4023  endif
4024 
4025  if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then
4026  v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
4027  flux_enter(i,j,4) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%t_bdry_val(i,j+1)* &
4028  cs%thickness_bdry_val(i,j+1)
4029  elseif (cs%v_face_mask(i,j+1) == 4.) then
4030  flux_enter(i,j,4) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j+1)*cs%t_bdry_val(i,j+1)
4031  endif
4032 
4033 ! if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1)) then
4034  ! this is solely for the purposes of keeping the mask consistent while advancing
4035  ! the front without having to call pass_var - if cell is empty and cell to left
4036  ! is ice-covered then this cell will become partly covered
4037  ! hmask(i,j) = 2
4038  ! elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1)) then
4039  ! this is solely for the purposes of keeping the mask consistent while advancing the
4040  ! front without having to call pass_var - if cell is empty and cell to left is
4041  ! ice-covered then this cell will become partly covered
4042 ! hmask(i,j) = 2
4043 ! endif
4044 
4045  endif
4046  endif
4047  enddo ! j loop
4048  endif
4049  enddo ! i loop
4050 

References slope_limiter().

Referenced by ice_shelf_temp().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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,
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(szdi_(g),szdj_(g),4), intent(inout)  flux_enter 
)
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 [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]flux_enterThe ice volume flux into the cell

Definition at line 1421 of file MOM_ice_shelf_dynamics.F90.

1421  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1422  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1423  real, intent(in) :: time_step !< The time step for this update [s].
1424  real, dimension(SZDI_(G),SZDJ_(G)), &
1425  intent(inout) :: hmask !< A mask indicating which tracer points are
1426  !! partly or fully covered by an ice-shelf
1427  real, dimension(SZDI_(G),SZDJ_(G)), &
1428  intent(in) :: h0 !< The initial ice shelf thicknesses [Z ~> m].
1429  real, dimension(SZDI_(G),SZDJ_(G)), &
1430  intent(inout) :: h_after_uflux !< The ice shelf thicknesses after
1431  !! the zonal mass fluxes [Z ~> m].
1432  real, dimension(SZDI_(G),SZDJ_(G),4), &
1433  intent(inout) :: flux_enter !< The ice volume flux into the cell
1434  !! through the 4 cell boundaries [Z m2 ~> m3].
1435 
1436  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
1437 
1438  ! if there is an input bdry condition, the thickness there will be set in initialization
1439 
1440  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
1441  !
1442  ! from left neighbor: flux_enter(:,:,1)
1443  ! from right neighbor: flux_enter(:,:,2)
1444  ! from bottom neighbor: flux_enter(:,:,3)
1445  ! from top neighbor: flux_enter(:,:,4)
1446  !
1447  ! o--- (4) ---o
1448  ! | |
1449  ! (1) (2)
1450  ! | |
1451  ! o--- (3) ---o
1452  !
1453 
1454  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
1455  integer :: i_off, j_off
1456  logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
1457  real, dimension(-2:2) :: stencil ! Thicknesses [Z ~> m].
1458  real :: u_face, & ! positive if out
1459  flux_diff_cell, phi, dxh, dyh, dxdyh
1460  character (len=1) :: debug_str
1461 
1462  is = g%isc-2 ; ie = g%iec+2 ; js = g%jsc ; je = g%jec
1463  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1464  i_off = g%idg_offset ; j_off = g%jdg_offset
1465 
1466  do j=jsd+1,jed-1
1467  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1468  ((j+j_off) >= g%domain%njhalo+1)) then ! based on mehmet's code - only if btw north & south boundaries
1469 
1470  stencil(:) = -1.
1471 ! if (i+i_off == G%domain%nihalo+G%domain%nihalo)
1472  do i=is,ie
1473 
1474  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1475  ((i+i_off) >= g%domain%nihalo+1)) then
1476 
1477  if (i+i_off == g%domain%nihalo+1) then
1478  at_west_bdry=.true.
1479  else
1480  at_west_bdry=.false.
1481  endif
1482 
1483  if (i+i_off == g%domain%niglobal+g%domain%nihalo) then
1484  at_east_bdry=.true.
1485  else
1486  at_east_bdry=.false.
1487  endif
1488 
1489  if (hmask(i,j) == 1) then
1490 
1491  dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
1492 
1493  h_after_uflux(i,j) = h0(i,j)
1494 
1495  stencil(:) = h0(i-2:i+2,j) ! fine as long has nx_halo >= 2
1496 
1497  flux_diff_cell = 0
1498 
1499  ! 1ST DO LEFT FACE
1500 
1501  if (cs%u_face_mask(i-1,j) == 4.) then
1502 
1503  flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i-1,j) / dxdyh
1504 
1505  else
1506 
1507  ! get u-velocity at center of left face
1508  u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
1509 
1510  if (u_face > 0) then !flux is into cell - we need info from h(i-2), h(i-1) if available
1511 
1512  ! i may not cover all the cases.. but i cover the realistic ones
1513 
1514  if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then ! at western bdry but there is a
1515  ! thickness bdry condition, and the stencil contains it
1516  stencil(-1) = cs%thickness_bdry_val(i-1,j)
1517  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
1518 
1519  elseif (hmask(i-1,j) * hmask(i-2,j) == 1) then ! h(i-2) and h(i-1) are valid
1520  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
1521  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
1522  (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
1523 
1524  else ! h(i-1) is valid
1525  ! (o.w. flux would most likely be out of cell)
1526  ! but h(i-2) is not
1527 
1528  flux_diff_cell = flux_diff_cell + abs(u_face) * (dyh * time_step / dxdyh) * stencil(-1)
1529 
1530  endif
1531 
1532  elseif (u_face < 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
1533  if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid
1534  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
1535  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
1536  (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
1537 
1538  else
1539  flux_diff_cell = flux_diff_cell - abs(u_face) * (dyh * time_step / dxdyh) * stencil(0)
1540 
1541  if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
1542  flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
1543  endif
1544  endif
1545  endif
1546  endif
1547 
1548  ! NEXT DO RIGHT FACE
1549 
1550  ! get u-velocity at center of right face
1551 
1552  if (cs%u_face_mask(i+1,j) == 4.) then
1553 
1554  flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i+1,j) / dxdyh
1555 
1556  else
1557 
1558  u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
1559 
1560  if (u_face < 0) then !flux is into cell - we need info from h(i+2), h(i+1) if available
1561 
1562  if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then ! at eastern bdry but there is a
1563  ! thickness bdry condition, and the stencil contains it
1564 
1565  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
1566 
1567  elseif (hmask(i+1,j) * hmask(i+2,j) == 1) then ! h(i+2) and h(i+1) are valid
1568 
1569  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
1570  flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
1571  (stencil(1) - phi * (stencil(1)-stencil(0))/2)
1572 
1573  else ! h(i+1) is valid
1574  ! (o.w. flux would most likely be out of cell)
1575  ! but h(i+2) is not
1576 
1577  flux_diff_cell = flux_diff_cell + abs(u_face) * (dyh * time_step / dxdyh) * stencil(1)
1578 
1579  endif
1580 
1581  elseif (u_face > 0) then !flux is out of cell - we need info from h(i-1), h(i+1) if available
1582 
1583  if (hmask(i-1,j) * hmask(i+1,j) == 1) then ! h(i-1) and h(i+1) are both valid
1584 
1585  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
1586  flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
1587  (stencil(0) - phi * (stencil(0)-stencil(1))/2)
1588 
1589  else ! h(i+1) is valid
1590  ! (o.w. flux would most likely be out of cell)
1591  ! but h(i+2) is not
1592 
1593  flux_diff_cell = flux_diff_cell - abs(u_face) * (dyh * time_step / dxdyh) * stencil(0)
1594 
1595  if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then
1596  flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step * stencil(0)
1597  endif
1598 
1599  endif
1600 
1601  endif
1602 
1603  h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
1604 
1605  endif
1606 
1607  elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
1608 
1609  if (at_west_bdry .AND. (hmask(i-1,j) == 3)) then
1610  u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
1611  flux_enter(i,j,1) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%thickness_bdry_val(i-1,j)
1612  elseif (cs%u_face_mask(i-1,j) == 4.) then
1613  flux_enter(i,j,1) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i-1,j)
1614  endif
1615 
1616  if (at_east_bdry .AND. (hmask(i+1,j) == 3)) then
1617  u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
1618  flux_enter(i,j,2) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%thickness_bdry_val(i+1,j)
1619  elseif (cs%u_face_mask(i+1,j) == 4.) then
1620  flux_enter(i,j,2) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i+1,j)
1621  endif
1622 
1623  if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1)) then
1624  ! this is solely for the purposes of keeping the mask consistent while advancing
1625  ! the front without having to call pass_var - if cell is empty and cell to left
1626  ! is ice-covered then this cell will become partly covered
1627 
1628  hmask(i,j) = 2
1629  elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1)) then
1630  ! this is solely for the purposes of keeping the mask consistent while advancing
1631  ! the front without having to call pass_var - if cell is empty and cell to left
1632  ! is ice-covered then this cell will become partly covered
1633 
1634  hmask(i,j) = 2
1635 
1636  endif
1637 
1638  endif
1639 
1640  endif
1641 
1642  enddo ! i loop
1643 
1644  endif
1645 
1646  enddo ! j loop
1647 

References slope_limiter().

Referenced by ice_shelf_advect().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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,
real, intent(in)  time_step,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  hmask,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_after_uflux,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  h_after_vflux,
real, dimension( g %isd: g %ied, g %jsd: g %jed,4), intent(inout)  flux_enter 
)
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 [s].
[in,out]hmaskA mask indicating which tracer points are
[in]h_after_ufluxThe ice shelf thicknesses after
[in,out]h_after_vfluxThe ice shelf thicknesses after
[in,out]flux_enterThe ice volume flux into the cell

Definition at line 1651 of file MOM_ice_shelf_dynamics.F90.

1651  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1652  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1653  real, intent(in) :: time_step !< The time step for this update [s].
1654  real, dimension(SZDI_(G),SZDJ_(G)), &
1655  intent(inout) :: hmask !< A mask indicating which tracer points are
1656  !! partly or fully covered by an ice-shelf
1657  real, dimension(SZDI_(G),SZDJ_(G)), &
1658  intent(in) :: h_after_uflux !< The ice shelf thicknesses after
1659  !! the zonal mass fluxes [Z ~> m].
1660  real, dimension(SZDI_(G),SZDJ_(G)), &
1661  intent(inout) :: h_after_vflux !< The ice shelf thicknesses after
1662  !! the meridional mass fluxes [Z ~> m].
1663  real, dimension(SZDI_(G),SZDJ_(G),4), &
1664  intent(inout) :: flux_enter !< The ice volume flux into the cell
1665  !! through the 4 cell boundaries [Z m2 ~> m3].
1666 
1667  ! use will be made of ISS%hmask here - its value at the boundary will be zero, just like uncovered cells
1668 
1669  ! if there is an input bdry condition, the thickness there will be set in initialization
1670 
1671  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
1672  !
1673  ! from left neighbor: flux_enter(:,:,1)
1674  ! from right neighbor: flux_enter(:,:,2)
1675  ! from bottom neighbor: flux_enter(:,:,3)
1676  ! from top neighbor: flux_enter(:,:,4)
1677  !
1678  ! o--- (4) ---o
1679  ! | |
1680  ! (1) (2)
1681  ! | |
1682  ! o--- (3) ---o
1683  !
1684 
1685  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
1686  integer :: i_off, j_off
1687  logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
1688  real, dimension(-2:2) :: stencil ! Thicknesses [Z ~> m].
1689  real :: v_face, & ! positive if out
1690  flux_diff_cell, phi, dxh, dyh, dxdyh
1691  character(len=1) :: debug_str
1692 
1693  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
1694  i_off = g%idg_offset ; j_off = g%jdg_offset
1695 
1696  do i=isd+2,ied-2
1697  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1698  ((i+i_off) >= g%domain%nihalo+1)) then ! based on Mehmet's code - only if btw east & west boundaries
1699 
1700  stencil(:) = -1
1701 
1702  do j=js,je
1703 
1704  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1705  ((j+j_off) >= g%domain%njhalo+1)) then
1706 
1707  if (j+j_off == g%domain%njhalo+1) then
1708  at_south_bdry=.true.
1709  else
1710  at_south_bdry=.false.
1711  endif
1712 
1713  if (j+j_off == g%domain%njglobal+g%domain%njhalo) then
1714  at_north_bdry=.true.
1715  else
1716  at_north_bdry=.false.
1717  endif
1718 
1719  if (hmask(i,j) == 1) then
1720  dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
1721  h_after_vflux(i,j) = h_after_uflux(i,j)
1722 
1723  stencil(:) = h_after_uflux(i,j-2:j+2) ! fine as long has ny_halo >= 2
1724  flux_diff_cell = 0
1725 
1726  ! 1ST DO south FACE
1727 
1728  if (cs%v_face_mask(i,j-1) == 4.) then
1729 
1730  flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j-1) / dxdyh
1731 
1732  else
1733 
1734  ! get u-velocity at center of left face
1735  v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
1736 
1737  if (v_face > 0) then !flux is into cell - we need info from h(j-2), h(j-1) if available
1738 
1739  ! i may not cover all the cases.. but i cover the realistic ones
1740 
1741  if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then ! at western bdry but there is a
1742  ! thickness bdry condition, and the stencil contains it
1743  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
1744 
1745  elseif (hmask(i,j-1) * hmask(i,j-2) == 1) then ! h(j-2) and h(j-1) are valid
1746 
1747  phi = slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
1748  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
1749  (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
1750 
1751  else ! h(j-1) is valid
1752  ! (o.w. flux would most likely be out of cell)
1753  ! but h(j-2) is not
1754  flux_diff_cell = flux_diff_cell + abs(v_face) * (dxh * time_step / dxdyh) * stencil(-1)
1755  endif
1756 
1757  elseif (v_face < 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
1758 
1759  if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid
1760  phi = slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
1761  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
1762  (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
1763  else
1764  flux_diff_cell = flux_diff_cell - abs(v_face) * (dxh * time_step / dxdyh) * stencil(0)
1765 
1766  if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
1767  flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
1768  endif
1769 
1770  endif
1771 
1772  endif
1773 
1774  endif
1775 
1776  ! NEXT DO north FACE
1777 
1778  if (cs%v_face_mask(i,j+1) == 4.) then
1779 
1780  flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j+1) / dxdyh
1781 
1782  else
1783 
1784  ! get u-velocity at center of right face
1785  v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
1786 
1787  if (v_face < 0) then !flux is into cell - we need info from h(j+2), h(j+1) if available
1788 
1789  if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then ! at eastern bdry but there is a
1790  ! thickness bdry condition, and the stencil contains it
1791  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
1792  elseif (hmask(i,j+1) * hmask(i,j+2) == 1) then ! h(j+2) and h(j+1) are valid
1793  phi = slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
1794  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
1795  (stencil(1) - phi * (stencil(1)-stencil(0))/2)
1796  else ! h(j+1) is valid
1797  ! (o.w. flux would most likely be out of cell)
1798  ! but h(j+2) is not
1799  flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
1800  endif
1801 
1802  elseif (v_face > 0) then !flux is out of cell - we need info from h(j-1), h(j+1) if available
1803 
1804  if (hmask(i,j-1) * hmask(i,j+1) == 1) then ! h(j-1) and h(j+1) are both valid
1805  phi = slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
1806  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
1807  (stencil(0) - phi * (stencil(0)-stencil(1))/2)
1808  else ! h(j+1) is valid
1809  ! (o.w. flux would most likely be out of cell)
1810  ! but h(j+2) is not
1811  flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
1812  if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
1813  flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
1814  endif
1815  endif
1816 
1817  endif
1818 
1819  endif
1820 
1821  h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
1822 
1823  elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2)) then
1824 
1825  if (at_south_bdry .AND. (hmask(i,j-1) == 3)) then
1826  v_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i,j-1))
1827  flux_enter(i,j,3) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%thickness_bdry_val(i,j-1)
1828  elseif (cs%v_face_mask(i,j-1) == 4.) then
1829  flux_enter(i,j,3) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j-1)
1830  endif
1831 
1832  if (at_north_bdry .AND. (hmask(i,j+1) == 3)) then
1833  v_face = 0.5 * (cs%u_shelf(i-1,j) + cs%u_shelf(i,j))
1834  flux_enter(i,j,4) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%thickness_bdry_val(i,j+1)
1835  elseif (cs%v_face_mask(i,j+1) == 4.) then
1836  flux_enter(i,j,4) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j+1)
1837  endif
1838 
1839  if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1)) then
1840  ! this is solely for the purposes of keeping the mask consistent while advancing
1841  ! the front without having to call pass_var - if cell is empty and cell to left
1842  ! is ice-covered then this cell will become partly covered
1843  hmask(i,j) = 2
1844  elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1)) then
1845  ! this is solely for the purposes of keeping the mask consistent while advancing
1846  ! the front without having to call pass_var - if cell is empty and cell to left
1847  ! is ice-covered then this cell will become partly covered
1848  hmask(i,j) = 2
1849  endif
1850 
1851  endif
1852  endif
1853  enddo ! j loop
1854  endif
1855  enddo ! i loop
1856 

References slope_limiter().

Referenced by ice_shelf_advect().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 3454 of file MOM_ice_shelf_dynamics.F90.

3454  type(ice_shelf_dyn_CS), pointer :: CS !< A pointer to the ice shelf dynamics control structure
3455 
3456  if (.not.associated(cs)) return
3457 
3458  deallocate(cs%u_shelf, cs%v_shelf)
3459  deallocate(cs%t_shelf, cs%tmask)
3460  deallocate(cs%u_bdry_val, cs%v_bdry_val, cs%t_bdry_val)
3461  deallocate(cs%u_face_mask, cs%v_face_mask)
3462  deallocate(cs%umask, cs%vmask)
3463 
3464  deallocate(cs%ice_visc, cs%taub_beta_eff)
3465  deallocate(cs%OD_rt, cs%OD_av)
3466  deallocate(cs%float_frac, cs%float_frac_rt)
3467 
3468  deallocate(cs)
3469 

Referenced by mom_ice_shelf::ice_shelf_end().

Here is the caller graph for this function:

◆ 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 
)

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 [m2].
[in,out]hmaskA mask indicating which tracer points are
[in]thickness_calveThe thickness at which to trigger calving [Z ~> m].

Definition at line 2039 of file MOM_ice_shelf_dynamics.F90.

2039  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2040  real, dimension(SZDI_(G),SZDJ_(G)), &
2041  intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m].
2042  real, dimension(SZDI_(G),SZDJ_(G)), &
2043  intent(inout) :: area_shelf_h !< The area per cell covered by the ice shelf [m2].
2044  real, dimension(SZDI_(G),SZDJ_(G)), &
2045  intent(inout) :: hmask !< A mask indicating which tracer points are
2046  !! partly or fully covered by an ice-shelf
2047  real, intent(in) :: thickness_calve !< The thickness at which to trigger calving [Z ~> m].
2048 
2049  integer :: i,j
2050 
2051  do j=g%jsd,g%jed
2052  do i=g%isd,g%ied
2053 ! if ((h_shelf(i,j) < CS%thickness_calve) .and. (hmask(i,j) == 1) .and. &
2054 ! (CS%float_frac(i,j) == 0.0)) then
2055  if ((h_shelf(i,j) < thickness_calve) .and. (area_shelf_h(i,j) > 0.)) then
2056  h_shelf(i,j) = 0.0
2057  area_shelf_h(i,j) = 0.0
2058  hmask(i,j) = 0.0
2059  endif
2060  enddo
2061  enddo
2062 

Referenced by ice_shelf_advect().

Here is the caller graph for this function:

◆ 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,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  u,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  v,
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(szdi_(g),szdj_(g),8,4), 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,out]uThe zonal ice shelf velocity at vertices [m year-1]
[in,out]vThe meridional ice shelf velocity at vertices [m year-1]
[in]taudxThe x-direction driving stress, in ???
[in]taudyThe y-direction driving stress, in ???
[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 tolerence
[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 1036 of file MOM_ice_shelf_dynamics.F90.

1036  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1037  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
1038  !! the ice-shelf state
1039  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
1040  real, dimension(SZDIB_(G),SZDJB_(G)), &
1041  intent(inout) :: u !< The zonal ice shelf velocity at vertices [m year-1]
1042  real, dimension(SZDIB_(G),SZDJB_(G)), &
1043  intent(inout) :: v !< The meridional ice shelf velocity at vertices [m year-1]
1044  real, dimension(SZDIB_(G),SZDJB_(G)), &
1045  intent(in) :: taudx !< The x-direction driving stress, in ???
1046  real, dimension(SZDIB_(G),SZDJB_(G)), &
1047  intent(in) :: taudy !< The y-direction driving stress, in ???
1048  real, dimension(SZDIB_(G),SZDJB_(G)), &
1049  intent(in) :: H_node !< The ice shelf thickness at nodal (corner)
1050  !! points [Z ~> m].
1051  real, dimension(SZDI_(G),SZDJ_(G)), &
1052  intent(in) :: float_cond !< An array indicating where the ice
1053  !! shelf is floating: 0 if floating, 1 if not.
1054  real, dimension(SZDI_(G),SZDJ_(G)), &
1055  intent(in) :: hmask !< A mask indicating which tracer points are
1056  !! partly or fully covered by an ice-shelf
1057  integer, intent(out) :: conv_flag !< A flag indicating whether (1) or not (0) the
1058  !! iterations have converged to the specified tolerence
1059  integer, intent(out) :: iters !< The number of iterations used in the solver.
1060  type(time_type), intent(in) :: Time !< The current model time
1061  real, dimension(SZDI_(G),SZDJ_(G),8,4), &
1062  intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian
1063  !! quadrature points surrounding the cell verticies.
1064  real, dimension(:,:,:,:,:,:), &
1065  intent(in) :: Phisub !< Quadrature structure weights at subgridscale
1066  !! locations for finite element calculations
1067 ! one linear solve (nonlinear iteration) of the solution for velocity
1068 
1069 ! in this subroutine:
1070 ! boundary contributions are added to taud to get the RHS
1071 ! diagonal of matrix is found (for Jacobi precondition)
1072 ! CG iteration is carried out for max. iterations or until convergence
1073 
1074 ! assumed - u, v, taud, visc, beta_eff are valid on the halo
1075 
1076  real, dimension(SZDIB_(G),SZDJB_(G)) :: &
1077  Ru, Rv, Zu, Zv, DIAGu, DIAGv, RHSu, RHSv, &
1078  ubd, vbd, Au, Av, Du, Dv, &
1079  Zu_old, Zv_old, Ru_old, Rv_old, &
1080  sum_vec, sum_vec_2
1081  integer :: iter, i, j, isd, ied, jsd, jed, &
1082  isc, iec, jsc, jec, is, js, ie, je, isumstart, jsumstart, &
1083  isdq, iedq, jsdq, jedq, iscq, iecq, jscq, jecq, nx_halo, ny_halo
1084  real :: tol, beta_k, alpha_k, area, dot_p1, dot_p2, resid0, cg_halo, dot_p1a, dot_p2a
1085  character(2) :: gridsize
1086 
1087  real, dimension(8,4) :: Phi_temp
1088  real, dimension(2,2) :: X,Y
1089 
1090  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
1091  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
1092  ny_halo = g%domain%njhalo ; nx_halo = g%domain%nihalo
1093  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1094  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1095 
1096  zu(:,:) = 0 ; zv(:,:) = 0 ; diagu(:,:) = 0 ; diagv(:,:) = 0
1097  ru(:,:) = 0 ; rv(:,:) = 0 ; au(:,:) = 0 ; av(:,:) = 0
1098  du(:,:) = 0 ; dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0
1099  dot_p1 = 0 ; dot_p2 = 0
1100 
1101  isumstart = g%isc
1102  ! Include the edge if tile is at the western bdry; Should add a test to avoid this if reentrant.
1103  if (g%isc+g%idg_offset==g%isg) isumstart = g%iscB
1104 
1105  jsumstart = g%jsc
1106  ! Include the edge if tile is at the southern bdry; Should add a test to avoid this if reentrant.
1107  if (g%jsc+g%jdg_offset==g%jsg) jsumstart = g%jscB
1108 
1109  call apply_boundary_values(cs, iss, g, time, phisub, h_node, cs%ice_visc, &
1110  cs%taub_beta_eff, float_cond, &
1111  cs%density_ice/cs%density_ocean_avg, ubd, vbd)
1112 
1113  rhsu(:,:) = taudx(:,:) - ubd(:,:)
1114  rhsv(:,:) = taudy(:,:) - vbd(:,:)
1115 
1116 
1117  call pass_vector(rhsu, rhsv, g%domain, to_all, bgrid_ne)
1118 
1119  call matrix_diagonal(cs, g, float_cond, h_node, cs%ice_visc, &
1120  cs%taub_beta_eff, hmask, &
1121  cs%density_ice/cs%density_ocean_avg, phisub, diagu, diagv)
1122 ! DIAGu(:,:) = 1 ; DIAGv(:,:) = 1
1123 
1124  call pass_vector(diagu, diagv, g%domain, to_all, bgrid_ne)
1125 
1126  call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, hmask, &
1127  h_node, cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, &
1128  g%US%L_to_m**2*g%areaT, g, isc-1, iec+1, jsc-1, jec+1, cs%density_ice/cs%density_ocean_avg)
1129 
1130  call pass_vector(au, av, g%domain, to_all, bgrid_ne)
1131 
1132  ru(:,:) = rhsu(:,:) - au(:,:) ; rv(:,:) = rhsv(:,:) - av(:,:)
1133 
1134  if (.not. cs%use_reproducing_sums) then
1135 
1136  do j=jsumstart,jecq
1137  do i=isumstart,iecq
1138  if (cs%umask(i,j) == 1) dot_p1 = dot_p1 + ru(i,j)**2
1139  if (cs%vmask(i,j) == 1) dot_p1 = dot_p1 + rv(i,j)**2
1140  enddo
1141  enddo
1142 
1143  call sum_across_pes(dot_p1)
1144 
1145  else
1146 
1147  sum_vec(:,:) = 0.0
1148 
1149  do j=jsumstart,jecq
1150  do i=isumstart,iecq
1151  if (cs%umask(i,j) == 1) sum_vec(i,j) = ru(i,j)**2
1152  if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
1153  enddo
1154  enddo
1155 
1156  dot_p1 = reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1157  jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1158 
1159  endif
1160 
1161  resid0 = sqrt(dot_p1)
1162 
1163  do j=jsdq,jedq
1164  do i=isdq,iedq
1165  if (cs%umask(i,j) == 1) zu(i,j) = ru(i,j) / diagu(i,j)
1166  if (cs%vmask(i,j) == 1) zv(i,j) = rv(i,j) / diagv(i,j)
1167  enddo
1168  enddo
1169 
1170  du(:,:) = zu(:,:) ; dv(:,:) = zv(:,:)
1171 
1172  cg_halo = 3
1173  conv_flag = 0
1174 
1175  !!!!!!!!!!!!!!!!!!
1176  !! !!
1177  !! MAIN CG LOOP !!
1178  !! !!
1179  !!!!!!!!!!!!!!!!!!
1180 
1181 
1182 
1183  ! initially, c-grid data is valid up to 3 halo nodes out
1184 
1185  do iter = 1,cs%cg_max_iterations
1186 
1187  ! assume asymmetry
1188  ! thus we can never assume that any arrays are legit more than 3 vertices past
1189  ! the computational domain - this is their state in the initial iteration
1190 
1191 
1192  is = isc - cg_halo ; ie = iecq + cg_halo
1193  js = jscq - cg_halo ; je = jecq + cg_halo
1194 
1195  au(:,:) = 0 ; av(:,:) = 0
1196 
1197  call cg_action(au, av, du, dv, phi, phisub, cs%umask, cs%vmask, hmask, &
1198  h_node, cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, &
1199  g%US%L_to_m**2*g%areaT, g, is, ie, js, je, cs%density_ice/cs%density_ocean_avg)
1200 
1201  ! Au, Av valid region moves in by 1
1202 
1203  if ( .not. cs%use_reproducing_sums) then
1204 
1205 
1206  ! alpha_k = (Z \dot R) / (D \dot AD}
1207  dot_p1 = 0 ; dot_p2 = 0
1208  do j=jsumstart,jecq
1209  do i=isumstart,iecq
1210  if (cs%umask(i,j) == 1) then
1211  dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
1212  dot_p2 = dot_p2 + du(i,j)*au(i,j)
1213  endif
1214  if (cs%vmask(i,j) == 1) then
1215  dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
1216  dot_p2 = dot_p2 + dv(i,j)*av(i,j)
1217  endif
1218  enddo
1219  enddo
1220  call sum_across_pes(dot_p1) ; call sum_across_pes(dot_p2)
1221  else
1222 
1223  sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
1224 
1225  do j=jscq,jecq
1226  do i=iscq,iecq
1227  if (cs%umask(i,j) == 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
1228  if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + zv(i,j) * rv(i,j)
1229 
1230  if (cs%umask(i,j) == 1) sum_vec_2(i,j) = du(i,j) * au(i,j)
1231  if (cs%vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + dv(i,j) * av(i,j)
1232  enddo
1233  enddo
1234 
1235  dot_p1 = reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1236  jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1237 
1238  dot_p2 = reproducing_sum( sum_vec_2, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1239  jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1240  endif
1241 
1242  alpha_k = dot_p1/dot_p2
1243 
1244  do j=jsd,jed
1245  do i=isd,ied
1246  if (cs%umask(i,j) == 1) u(i,j) = u(i,j) + alpha_k * du(i,j)
1247  if (cs%vmask(i,j) == 1) v(i,j) = v(i,j) + alpha_k * dv(i,j)
1248  enddo
1249  enddo
1250 
1251  do j=jsd,jed
1252  do i=isd,ied
1253  if (cs%umask(i,j) == 1) then
1254  ru_old(i,j) = ru(i,j) ; zu_old(i,j) = zu(i,j)
1255  endif
1256  if (cs%vmask(i,j) == 1) then
1257  rv_old(i,j) = rv(i,j) ; zv_old(i,j) = zv(i,j)
1258  endif
1259  enddo
1260  enddo
1261 
1262 ! Ru(:,:) = Ru(:,:) - alpha_k * Au(:,:)
1263 ! Rv(:,:) = Rv(:,:) - alpha_k * Av(:,:)
1264 
1265  do j=jsd,jed
1266  do i=isd,ied
1267  if (cs%umask(i,j) == 1) ru(i,j) = ru(i,j) - alpha_k * au(i,j)
1268  if (cs%vmask(i,j) == 1) rv(i,j) = rv(i,j) - alpha_k * av(i,j)
1269  enddo
1270  enddo
1271 
1272 
1273  do j=jsdq,jedq
1274  do i=isdq,iedq
1275  if (cs%umask(i,j) == 1) then
1276  zu(i,j) = ru(i,j) / diagu(i,j)
1277  endif
1278  if (cs%vmask(i,j) == 1) then
1279  zv(i,j) = rv(i,j) / diagv(i,j)
1280  endif
1281  enddo
1282  enddo
1283 
1284  ! R,u,v,Z valid region moves in by 1
1285 
1286  if (.not. cs%use_reproducing_sums) then
1287 
1288  ! beta_k = (Z \dot R) / (Zold \dot Rold}
1289  dot_p1 = 0 ; dot_p2 = 0
1290  do j=jsumstart,jecq
1291  do i=isumstart,iecq
1292  if (cs%umask(i,j) == 1) then
1293  dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
1294  dot_p2 = dot_p2 + zu_old(i,j)*ru_old(i,j)
1295  endif
1296  if (cs%vmask(i,j) == 1) then
1297  dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
1298  dot_p2 = dot_p2 + zv_old(i,j)*rv_old(i,j)
1299  endif
1300  enddo
1301  enddo
1302  call sum_across_pes(dot_p1) ; call sum_across_pes(dot_p2)
1303 
1304 
1305  else
1306 
1307  sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
1308 
1309  do j=jsumstart,jecq
1310  do i=isumstart,iecq
1311  if (cs%umask(i,j) == 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
1312  if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + &
1313  zv(i,j) * rv(i,j)
1314 
1315  if (cs%umask(i,j) == 1) sum_vec_2(i,j) = zu_old(i,j) * ru_old(i,j)
1316  if (cs%vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + &
1317  zv_old(i,j) * rv_old(i,j)
1318  enddo
1319  enddo
1320 
1321 
1322  dot_p1 = reproducing_sum(sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1323  jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1324 
1325  dot_p2 = reproducing_sum(sum_vec_2, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1326  jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1327 
1328  endif
1329 
1330  beta_k = dot_p1/dot_p2
1331 
1332 
1333 ! Du(:,:) = Zu(:,:) + beta_k * Du(:,:)
1334 ! Dv(:,:) = Zv(:,:) + beta_k * Dv(:,:)
1335 
1336  do j=jsd,jed
1337  do i=isd,ied
1338  if (cs%umask(i,j) == 1) du(i,j) = zu(i,j) + beta_k * du(i,j)
1339  if (cs%vmask(i,j) == 1) dv(i,j) = zv(i,j) + beta_k * dv(i,j)
1340  enddo
1341  enddo
1342 
1343  ! D valid region moves in by 1
1344 
1345  dot_p1 = 0
1346 
1347  if (.not. cs%use_reproducing_sums) then
1348 
1349  do j=jsumstart,jecq
1350  do i=isumstart,iecq
1351  if (cs%umask(i,j) == 1) then
1352  dot_p1 = dot_p1 + ru(i,j)**2
1353  endif
1354  if (cs%vmask(i,j) == 1) then
1355  dot_p1 = dot_p1 + rv(i,j)**2
1356  endif
1357  enddo
1358  enddo
1359  call sum_across_pes(dot_p1)
1360 
1361  else
1362 
1363  sum_vec(:,:) = 0.0
1364 
1365  do j=jsumstart,jecq
1366  do i=isumstart,iecq
1367  if (cs%umask(i,j) == 1) sum_vec(i,j) = ru(i,j)**2
1368  if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
1369  enddo
1370  enddo
1371 
1372  dot_p1 = reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1373  jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1374  endif
1375 
1376  dot_p1 = sqrt(dot_p1)
1377 
1378  if (dot_p1 <= cs%cg_tolerance * resid0) then
1379  iters = iter
1380  conv_flag = 1
1381  exit
1382  endif
1383 
1384  cg_halo = cg_halo - 1
1385 
1386  if (cg_halo == 0) then
1387  ! pass vectors
1388  call pass_vector(du, dv, g%domain, to_all, bgrid_ne)
1389  call pass_vector(u, v, g%domain, to_all, bgrid_ne)
1390  call pass_vector(ru, rv, g%domain, to_all, bgrid_ne)
1391  cg_halo = 3
1392  endif
1393 
1394  enddo ! end of CG loop
1395 
1396  do j=jsdq,jedq
1397  do i=isdq,iedq
1398  if (cs%umask(i,j) == 3) then
1399  u(i,j) = cs%u_bdry_val(i,j)
1400  elseif (cs%umask(i,j) == 0) then
1401  u(i,j) = 0
1402  endif
1403 
1404  if (cs%vmask(i,j) == 3) then
1405  v(i,j) = cs%v_bdry_val(i,j)
1406  elseif (cs%vmask(i,j) == 0) then
1407  v(i,j) = 0
1408  endif
1409  enddo
1410  enddo
1411 
1412  call pass_vector(u,v, g%domain, to_all, bgrid_ne)
1413 
1414  if (conv_flag == 0) then
1415  iters = cs%cg_max_iterations
1416  endif
1417 

References apply_boundary_values(), cg_action(), matrix_diagonal(), and mom_domains::to_all.

Referenced by ice_shelf_solve_outer().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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,
real, dimension(szdib_(g),szdjb_(g)), intent(inout)  v,
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]usPointer to a structure containing unit conversion factors
[in,out]uThe zonal ice shelf velocity at vertices [m year-1]
[in,out]vThe meridional ice shelf velocity at vertices [m year-1]
[out]itersThe number of iterations used in the solver.
[in]timeThe current model time

Definition at line 773 of file MOM_ice_shelf_dynamics.F90.

773  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
774  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
775  !! the ice-shelf state
776  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
777  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
778  real, dimension(SZDIB_(G),SZDJB_(G)), &
779  intent(inout) :: u !< The zonal ice shelf velocity at vertices [m year-1]
780  real, dimension(SZDIB_(G),SZDJB_(G)), &
781  intent(inout) :: v !< The meridional ice shelf velocity at vertices [m year-1]
782  integer, intent(out) :: iters !< The number of iterations used in the solver.
783  type(time_type), intent(in) :: Time !< The current model time
784 
785  real, dimension(SZDIB_(G),SZDJB_(G)) :: TAUDX, TAUDY, u_prev_iterate, v_prev_iterate, &
786  u_bdry_cont, v_bdry_cont, Au, Av, err_u, err_v, &
787  u_last, v_last
788  real, dimension(SZDIB_(G),SZDJB_(G)) :: H_node ! Ice shelf thickness at corners [Z ~> m].
789  real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! An array indicating where the ice
790  ! shelf is floating: 0 if floating, 1 if not.
791  character(len=160) :: mesg ! The text of an error message
792  integer :: conv_flag, i, j, k,l, iter
793  integer :: isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub
794  real :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi_rhow
795  real, pointer, dimension(:,:,:,:) :: Phi => null()
796  real, pointer, dimension(:,:,:,:,:,:) :: Phisub => null()
797  real, dimension(8,4) :: Phi_temp
798  real, dimension(2,2) :: X,Y
799  character(2) :: iternum
800  character(2) :: numproc
801 
802  ! for GL interpolation - need to make this a readable parameter
803  nsub = cs%n_sub_regularize
804 
805  isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
806  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
807  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
808 
809  taudx(:,:) = 0.0 ; taudy(:,:) = 0.0
810  u_bdry_cont(:,:) = 0.0 ; v_bdry_cont(:,:) = 0.0
811  au(:,:) = 0.0 ; av(:,:) = 0.0
812 
813  ! need to make these conditional on GL interpolation
814  float_cond(:,:) = 0.0 ; h_node(:,:)=0
815  allocate(phisub(nsub,nsub,2,2,2,2)) ; phisub = 0.0
816 
817  isumstart = g%isc
818  ! Include the edge if tile is at the western bdry; Should add a test to avoid this if reentrant.
819  if (g%isc+g%idg_offset==g%isg) isumstart = g%iscB
820 
821  jsumstart = g%jsc
822  ! Include the edge if tile is at the southern bdry; Should add a test to avoid this if reentrant.
823  if (g%jsc+g%jdg_offset==g%jsg) jsumstart = g%jscB
824 
825  call calc_shelf_driving_stress(cs, iss, g, us, taudx, taudy, cs%OD_av)
826 
827  ! this is to determine which cells contain the grounding line,
828  ! the criterion being that the cell is ice-covered, with some nodes
829  ! floating and some grounded
830  ! floatation condition is estimated by assuming topography is cellwise constant
831  ! and H is bilinear in a cell; floating where rho_i/rho_w * H_node + D is nonpositive
832 
833  ! need to make this conditional on GL interp
834 
835  if (cs%GL_regularize) then
836 
837  call interpolate_h_to_b(g, iss%h_shelf, iss%hmask, h_node)
838 
839  do j=g%jsc,g%jec
840  do i=g%isc,g%iec
841  nodefloat = 0
842  do k=0,1
843  do l=0,1
844  if ((iss%hmask(i,j) == 1) .and. &
845  (rhoi_rhow * h_node(i-1+k,j-1+l) - g%bathyT(i,j) <= 0)) then
846  nodefloat = nodefloat + 1
847  endif
848  enddo
849  enddo
850  if ((nodefloat > 0) .and. (nodefloat < 4)) then
851  float_cond(i,j) = 1.0
852  cs%float_frac(i,j) = 1.0
853  endif
854  enddo
855  enddo
856 
857  call pass_var(float_cond, g%Domain)
858 
859  call bilinear_shape_functions_subgrid(phisub, nsub)
860 
861  endif
862 
863  ! make above conditional
864 
865  u_prev_iterate(:,:) = u(:,:)
866  v_prev_iterate(:,:) = v(:,:)
867 
868  ! must prepare phi
869  allocate(phi(isd:ied,jsd:jed,1:8,1:4)) ; phi(:,:,:,:) = 0.0
870 
871  do j=jsd,jed ; do i=isd,ied
872  if (((i > isd) .and. (j > jsd))) then
873  x(:,:) = g%geoLonBu(i-1:i,j-1:j)*1000
874  y(:,:) = g%geoLatBu(i-1:i,j-1:j)*1000
875  else
876  x(2,:) = g%geoLonBu(i,j)*1000
877  x(1,:) = g%geoLonBu(i,j)*1000 - us%L_to_m*g%dxT(i,j)
878  y(:,2) = g%geoLatBu(i,j)*1000
879  y(:,1) = g%geoLatBu(i,j)*1000 - us%L_to_m*g%dyT(i,j)
880  endif
881 
882  call bilinear_shape_functions(x, y, phi_temp, area)
883  phi(i,j,:,:) = phi_temp
884  enddo ; enddo
885 
886  call calc_shelf_visc(cs, iss, g, us, u, v)
887 
888  call pass_var(cs%ice_visc, g%domain)
889  call pass_var(cs%taub_beta_eff, g%domain)
890 
891  ! makes sure basal stress is only applied when it is supposed to be
892 
893  do j=g%jsd,g%jed ; do i=g%isd,g%ied
894  cs%taub_beta_eff(i,j) = cs%taub_beta_eff(i,j) * cs%float_frac(i,j)
895  enddo ; enddo
896 
897  call apply_boundary_values(cs, iss, g, time, phisub, h_node, cs%ice_visc, &
898  cs%taub_beta_eff, float_cond, &
899  rhoi_rhow, u_bdry_cont, v_bdry_cont)
900 
901  au(:,:) = 0.0 ; av(:,:) = 0.0
902 
903  call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
904  cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, g%US%L_to_m**2*g%areaT, &
905  g, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
906 
907  err_init = 0 ; err_tempu = 0; err_tempv = 0
908  do j=jsumstart,g%jecB
909  do i=isumstart,g%iecB
910  if (cs%umask(i,j) == 1) then
911  err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
912  endif
913  if (cs%vmask(i,j) == 1) then
914  err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
915  endif
916  if (err_tempv >= err_init) then
917  err_init = err_tempv
918  endif
919  enddo
920  enddo
921 
922  call max_across_pes(err_init)
923 
924  write(mesg,*) "ice_shelf_solve_outer: INITIAL nonlinear residual = ",err_init
925  call mom_mesg(mesg, 5)
926 
927  u_last(:,:) = u(:,:) ; v_last(:,:) = v(:,:)
928 
929  !! begin loop
930 
931  do iter=1,100
932 
933  call ice_shelf_solve_inner(cs, iss, g, u, v, taudx, taudy, h_node, float_cond, &
934  iss%hmask, conv_flag, iters, time, phi, phisub)
935 
936  if (cs%debug) then
937  call qchksum(u, "u shelf", g%HI, haloshift=2)
938  call qchksum(v, "v shelf", g%HI, haloshift=2)
939  endif
940 
941  write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations"
942  call mom_mesg(mesg, 5)
943 
944  call calc_shelf_visc(cs, iss, g, us, u, v)
945  call pass_var(cs%ice_visc, g%domain)
946  call pass_var(cs%taub_beta_eff, g%domain)
947 
948  ! makes sure basal stress is only applied when it is supposed to be
949 
950  do j=g%jsd,g%jed ; do i=g%isd,g%ied
951  cs%taub_beta_eff(i,j) = cs%taub_beta_eff(i,j) * cs%float_frac(i,j)
952  enddo ; enddo
953 
954  u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0
955 
956  call apply_boundary_values(cs, iss, g, time, phisub, h_node, cs%ice_visc, &
957  cs%taub_beta_eff, float_cond, &
958  rhoi_rhow, u_bdry_cont, v_bdry_cont)
959 
960  au(:,:) = 0 ; av(:,:) = 0
961 
962  call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
963  cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, g%US%L_to_m**2*g%areaT, &
964  g, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
965 
966  err_max = 0
967 
968  if (cs%nonlin_solve_err_mode == 1) then
969 
970  do j=jsumstart,g%jecB
971  do i=isumstart,g%iecB
972  if (cs%umask(i,j) == 1) then
973  err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
974  endif
975  if (cs%vmask(i,j) == 1) then
976  err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
977  endif
978  if (err_tempv >= err_max) then
979  err_max = err_tempv
980  endif
981  enddo
982  enddo
983 
984  call max_across_pes(err_max)
985 
986  elseif (cs%nonlin_solve_err_mode == 2) then
987 
988  max_vel = 0 ; tempu = 0 ; tempv = 0
989 
990  do j=jsumstart,g%jecB
991  do i=isumstart,g%iecB
992  if (cs%umask(i,j) == 1) then
993  err_tempu = abs(u_last(i,j)-u(i,j))
994  tempu = u(i,j)
995  endif
996  if (cs%vmask(i,j) == 1) then
997  err_tempv = max(abs(v_last(i,j)- v(i,j)), err_tempu)
998  tempv = sqrt(v(i,j)**2+tempu**2)
999  endif
1000  if (err_tempv >= err_max) then
1001  err_max = err_tempv
1002  endif
1003  if (tempv >= max_vel) then
1004  max_vel = tempv
1005  endif
1006  enddo
1007  enddo
1008 
1009  u_last(:,:) = u(:,:)
1010  v_last(:,:) = v(:,:)
1011 
1012  call max_across_pes(max_vel)
1013  call max_across_pes(err_max)
1014  err_init = max_vel
1015 
1016  endif
1017 
1018  write(mesg,*) "ice_shelf_solve_outer: nonlinear residual = ",err_max/err_init
1019  call mom_mesg(mesg, 5)
1020 
1021  if (err_max <= cs%nonlinear_tolerance * err_init) then
1022  write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations"
1023  call mom_mesg(mesg, 5)
1024  exit
1025  endif
1026 
1027  enddo
1028 
1029  deallocate(phi)
1030  deallocate(phisub)
1031 

References apply_boundary_values(), bilinear_shape_functions(), bilinear_shape_functions_subgrid(), calc_shelf_driving_stress(), calc_shelf_visc(), cg_action(), ice_shelf_solve_inner(), interpolate_h_to_b(), and mom_error_handler::mom_mesg().

Referenced by initialize_diagnostic_fields(), initialize_ice_shelf_dyn(), and update_ice_shelf().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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( g %isd: g %ied, g %jsd: g %jed), 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]usPointer to a structure containing unit conversion factors
[in]time_stepThe time step for this update [s].
[in]melt_ratebasal melt rate [kg m-2 s-1]
[in]timeThe current model time

Definition at line 3475 of file MOM_ice_shelf_dynamics.F90.

3475  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
3476  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
3477  !! the ice-shelf state
3478  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
3479  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
3480  real, intent(in) :: time_step !< The time step for this update [s].
3481  real, dimension(SZDI_(G),SZDJ_(G)), &
3482  intent(in) :: melt_rate !< basal melt rate [kg m-2 s-1]
3483  type(time_type), intent(in) :: Time !< The current model time
3484 
3485 ! 5/23/12 OVS
3486 ! This subroutine takes the velocity (on the Bgrid) and timesteps
3487 ! (HT)_t = - div (uHT) + (adot Tsurf -bdot Tbot) once and then calculates T=HT/H
3488 !
3489 ! The flux overflows are included here. That is because they will be used to advect 3D scalars
3490 ! into partial cells
3491 
3492  !
3493  ! flux_enter: this is to capture flow into partially covered cells; it gives the mass flux into a given
3494  ! cell across its boundaries.
3495  ! ###Perhaps flux_enter should be changed into u-face and v-face
3496  ! ###fluxes, which can then be used in halo updates, etc.
3497  !
3498  ! from left neighbor: flux_enter(:,:,1)
3499  ! from right neighbor: flux_enter(:,:,2)
3500  ! from bottom neighbor: flux_enter(:,:,3)
3501  ! from top neighbor: flux_enter(:,:,4)
3502  !
3503  ! THESE ARE NOT CONSISTENT ==> FIND OUT WHAT YOU IMPLEMENTED
3504 
3505  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
3506  !
3507  ! o--- (4) ---o
3508  ! | |
3509  ! (1) (2)
3510  ! | |
3511  ! o--- (3) ---o
3512  !
3513 
3514  real, dimension(SZDI_(G),SZDJ_(G)) :: th_after_uflux, th_after_vflux, TH
3515  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter
3516  integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
3517  real :: rho, spy, t_bd, Tsurf, adot
3518 
3519  rho = cs%density_ice
3520  spy = 365 * 86400 ! seconds per year; is there a global constant for this? No - it is dependent upon a calendar.
3521 
3522  adot = 0.1*us%m_to_Z/spy ! for now adot and Tsurf are defined here adot=surf acc 0.1m/yr, Tsurf=-20oC, vary them later
3523  tsurf = -20.0
3524 
3525  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3526  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
3527  flux_enter(:,:,:) = 0.0
3528 
3529  th_after_uflux(:,:) = 0.0
3530  th_after_vflux(:,:) = 0.0
3531 
3532  do j=jsd,jed
3533  do i=isd,ied
3534  t_bd = cs%t_bdry_val(i,j)
3535 ! if (ISS%hmask(i,j) > 1) then
3536  if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2)) then
3537  cs%t_shelf(i,j) = cs%t_bdry_val(i,j)
3538  endif
3539  enddo
3540  enddo
3541 
3542  do j=jsd,jed
3543  do i=isd,ied
3544  th(i,j) = cs%t_shelf(i,j)*iss%h_shelf(i,j)
3545  enddo
3546  enddo
3547 
3548 
3549 ! call enable_averaging(time_step,Time,CS%diag)
3550 ! call pass_var(h_after_uflux, G%domain)
3551 ! call pass_var(h_after_vflux, G%domain)
3552 ! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag)
3553 ! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag)
3554 ! call disable_averaging(CS%diag)
3555 
3556  call ice_shelf_advect_temp_x(cs, g, time_step/spy, iss%hmask, th, th_after_uflux, flux_enter)
3557  call ice_shelf_advect_temp_y(cs, g, time_step/spy, iss%hmask, th_after_uflux, th_after_vflux, flux_enter)
3558 
3559  do j=jsd,jed
3560  do i=isd,ied
3561 ! if (ISS%hmask(i,j) == 1) then
3562  if (iss%h_shelf(i,j) > 0.0) then
3563  cs%t_shelf(i,j) = th_after_vflux(i,j)/(iss%h_shelf(i,j))
3564  else
3565  cs%t_shelf(i,j) = -10.0
3566  endif
3567  enddo
3568  enddo
3569 
3570  do j=jsd,jed
3571  do i=isd,ied
3572  t_bd = cs%t_bdry_val(i,j)
3573 ! if (ISS%hmask(i,j) > 1) then
3574  if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2)) then
3575  cs%t_shelf(i,j) = t_bd
3576 ! CS%t_shelf(i,j) = -15.0
3577  endif
3578  enddo
3579  enddo
3580 
3581  do j=jsc,jec
3582  do i=isc,iec
3583  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
3584  if (iss%h_shelf(i,j) > 0.0) then
3585 ! CS%t_shelf(i,j) = CS%t_shelf(i,j) + &
3586 ! time_step*(adot*Tsurf - US%m_to_Z*melt_rate(i,j)*ISS%tfreeze(i,j))/(ISS%h_shelf(i,j))
3587  cs%t_shelf(i,j) = cs%t_shelf(i,j) + &
3588  time_step*(adot*tsurf - (3.0*us%m_to_Z/spy)*iss%tfreeze(i,j)) / iss%h_shelf(i,j)
3589  else
3590  ! the ice is about to melt away
3591  ! in this case set thickness, area, and mask to zero
3592  ! NOTE: not mass conservative
3593  ! should maybe scale salt & heat flux for this cell
3594 
3595  cs%t_shelf(i,j) = -10.0
3596  cs%tmask(i,j) = 0.0
3597  endif
3598  endif
3599  enddo
3600  enddo
3601 
3602  call pass_var(cs%t_shelf, g%domain)
3603  call pass_var(cs%tmask, g%domain)
3604 
3605  if (cs%debug) then
3606  call hchksum(cs%t_shelf, "temp after front", g%HI, haloshift=3)
3607  endif
3608 

References ice_shelf_advect_temp_x(), and ice_shelf_advect_temp_y().

Referenced by update_ice_shelf().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 timestep that can be taken based on the current ice velocities. Because it involves finding a global minimum, it can be suprisingly 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 [s].

Definition at line 576 of file MOM_ice_shelf_dynamics.F90.

576  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
577  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
578  !! the ice-shelf state
579  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
580  real :: ice_time_step_CFL !< The maximum permitted timestep based on the ice velocities [s].
581 
582  real :: ratio, min_ratio
583  real :: local_u_max, local_v_max
584  integer :: i, j
585 
586  min_ratio = 1.0e16 ! This is just an arbitrary large nondiensional value.
587  do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (iss%hmask(i,j) == 1.0) then
588  local_u_max = max(abs(cs%u_shelf(i,j)), abs(cs%u_shelf(i+1,j+1)), &
589  abs(cs%u_shelf(i+1,j)), abs(cs%u_shelf(i,j+1)))
590  local_v_max = max(abs(cs%v_shelf(i,j)), abs(cs%v_shelf(i+1,j+1)), &
591  abs(cs%v_shelf(i+1,j)), abs(cs%v_shelf(i,j+1)))
592 
593  ! Here the hard-coded 1e-12 has units of m s-1. Consider revising.
594  ratio = g%US%L_to_m**2*min(g%areaT(i,j) / (local_u_max + 1.0e-12), &
595  g%areaT(i,j) / (local_v_max + 1.0e-12))
596  min_ratio = min(min_ratio, ratio)
597  endif ; enddo ; enddo ! i- and j- loops
598 
599  call min_across_pes(min_ratio)
600 
601  ! solved velocities are in m/yr; we want time_step_int in seconds
602  ice_time_step_cfl = cs%CFL_factor * min_ratio * (365*86400)
603 

Referenced by mom_ice_shelf::solo_time_step().

Here is the caller graph for this function:

◆ 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 [Z m2 s-1 ~> m3 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 2284 of file MOM_ice_shelf_dynamics.F90.

2284  type(ice_shelf_dyn_CS),intent(inout) :: CS !< A pointer to the ice shelf control structure
2285  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
2286  type(time_type), intent(in) :: Time !< The current model time
2287  real, dimension(SZDI_(G),SZDJ_(G)), &
2288  intent(in) :: hmask !< A mask indicating which tracer points are
2289  !! partly or fully covered by an ice-shelf
2290  real, intent(in) :: input_flux !< The integrated inward ice thickness flux [Z m2 s-1 ~> m3 s-1]
2291  real, intent(in) :: input_thick !< The ice thickness at boundaries [Z ~> m].
2292  logical, optional, intent(in) :: new_sim !< If present and false, this run is being restarted
2293 
2294 ! this will be a per-setup function. the boundary values of thickness and velocity
2295 ! (and possibly other variables) will be updated in this function
2296 
2297 ! FOR RESTARTING PURPOSES: if grid is not symmetric and the model is restarted, we will
2298 ! need to update those velocity points not *technically* in any
2299 ! computational domain -- if this function gets moves to another module,
2300 ! DO NOT TAKE THE RESTARTING BIT WITH IT
2301  integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
2302  integer :: gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
2303  integer :: i_off, j_off
2304  real :: A, n, ux, uy, vx, vy, eps_min, domain_width
2305 
2306 
2307  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2308 ! iscq = G%iscq ; iecq = G%iecq ; jscq = G%jscq ; jecq = G%jecq
2309  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
2310 ! iegq = G%iegq ; jegq = G%jegq
2311  i_off = g%idg_offset ; j_off = g%jdg_offset
2312 
2313  domain_width = g%len_lat
2314 
2315  ! this loop results in some values being set twice but... eh.
2316 
2317  do j=jsd,jed
2318  do i=isd,ied
2319 
2320  if (hmask(i,j) == 3) then
2321  cs%thickness_bdry_val(i,j) = input_thick
2322  endif
2323 
2324  if ((hmask(i,j) == 0) .or. (hmask(i,j) == 1) .or. (hmask(i,j) == 2)) then
2325  if ((i <= iec).and.(i >= isc)) then
2326  if (cs%u_face_mask(i-1,j) == 3) then
2327  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) * &
2328  1.5 * input_flux / input_thick
2329  cs%u_bdry_val(i-1,j) = (1 - ((g%geoLatBu(i-1,j) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
2330  1.5 * input_flux / input_thick
2331  endif
2332  endif
2333  endif
2334 
2335  if (.not.(new_sim)) then
2336  if (.not. g%symmetric) then
2337  if (((i+i_off) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3)) then
2338  cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
2339  cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
2340  endif
2341  if (((j+j_off) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3)) then
2342  cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
2343  cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
2344  endif
2345  endif
2346  endif
2347  enddo
2348  enddo
2349 

◆ 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]usPointer to a structure containing unit conversion factors
[in]timeThe current model time

Definition at line 540 of file MOM_ice_shelf_dynamics.F90.

540  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
541  type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe
542  !! the ice-shelf state
543  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
544  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
545  type(time_type), intent(in) :: Time !< The current model time
546 
547  integer :: i, j, iters, isd, ied, jsd, jed
548  real :: rhoi_rhow, OD
549  type(time_type) :: dummy_time
550 
551  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
552  dummy_time = set_time(0,0)
553  isd=g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
554 
555  do j=jsd,jed
556  do i=isd,ied
557  od = g%bathyT(i,j) - rhoi_rhow * iss%h_shelf(i,j)
558  if (od >= 0) then
559  ! ice thickness does not take up whole ocean column -> floating
560  cs%OD_av(i,j) = od
561  cs%float_frac(i,j) = 0.
562  else
563  cs%OD_av(i,j) = 0.
564  cs%float_frac(i,j) = 1.
565  endif
566  enddo
567  enddo
568 
569  call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, dummy_time)
570 

References ice_shelf_solve_outer().

Here is the call graph for this function:

◆ 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]usPointer to a 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 263 of file MOM_ice_shelf_dynamics.F90.

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

References ice_shelf_solve_outer(), mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), mom_domains::to_all, update_od_ffrac_uncoupled(), and update_velocity_masks().

Here is the call graph for this function:

◆ interpolate_h_to_b()

subroutine mom_ice_shelf_dynamics::interpolate_h_to_b ( type(ocean_grid_type), intent(inout)  G,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_shelf,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  hmask,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), 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 3410 of file MOM_ice_shelf_dynamics.F90.

3410  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
3411  real, dimension(SZDI_(G),SZDJ_(G)), &
3412  intent(in) :: h_shelf !< The ice shelf thickness at tracer points [Z ~> m].
3413  real, dimension(SZDI_(G),SZDJ_(G)), &
3414  intent(in) :: hmask !< A mask indicating which tracer points are
3415  !! partly or fully covered by an ice-shelf
3416  real, dimension(SZDIB_(G),SZDJB_(G)), &
3417  intent(inout) :: H_node !< The ice shelf thickness at nodal (corner)
3418  !! points [Z ~> m].
3419 
3420  integer :: i, j, isc, iec, jsc, jec, num_h, k, l
3421  real :: summ
3422 
3423  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3424 
3425  h_node(:,:) = 0.0
3426 
3427  ! H_node is node-centered; average over all cells that share that node
3428  ! if no (active) cells share the node then its value there is irrelevant
3429 
3430  do j=jsc-1,jec
3431  do i=isc-1,iec
3432  summ = 0.0
3433  num_h = 0
3434  do k=0,1
3435  do l=0,1
3436  if (hmask(i+k,j+l) == 1.0) then
3437  summ = summ + h_shelf(i+k,j+l)
3438  num_h = num_h + 1
3439  endif
3440  enddo
3441  enddo
3442  if (num_h > 0) then
3443  h_node(i,j) = summ / num_h
3444  endif
3445  enddo
3446  enddo
3447 
3448  call pass_var(h_node, g%domain, position=corner)
3449 

Referenced by ice_shelf_solve_outer().

Here is the caller graph for this function:

◆ matrix_diagonal()

subroutine mom_ice_shelf_dynamics::matrix_diagonal ( type(ice_shelf_dyn_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  float_cond,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(in)  H_node,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(in)  nu,
real, dimension(szdib_(g),szdjb_(g)), intent(in)  beta,
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]float_condAn array indicating where the ice
[in]h_nodeThe ice shelf thickness at nodal
[in]nuA field related to the ice viscosity from Glen's
[in]betaA 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, nondimensional
[in]phisubQuadrature structure weights at subgridscale locations for finite element calculations
[in,out]u_diagonalThe diagonal elements of the u-velocity
[in,out]v_diagonalThe diagonal elements of the v-velocity

Definition at line 2624 of file MOM_ice_shelf_dynamics.F90.

2624 
2625  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
2626  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
2627  real, dimension(SZDI_(G),SZDJ_(G)), &
2628  intent(in) :: float_cond !< An array indicating where the ice
2629  !! shelf is floating: 0 if floating, 1 if not.
2630  real, dimension(SZDIB_(G),SZDJB_(G)), &
2631  intent(in) :: H_node !< The ice shelf thickness at nodal
2632  !! (corner) points [Z ~> m].
2633  real, dimension(SZDIB_(G),SZDJB_(G)), &
2634  intent(in) :: nu !< A field related to the ice viscosity from Glen's
2635  !! flow law. The exact form and units depend on the
2636  !! basal law exponent.
2637  real, dimension(SZDIB_(G),SZDJB_(G)), &
2638  intent(in) :: beta !< A field related to the nonlinear part of the
2639  !! "linearized" basal stress. The exact form and
2640  !! units depend on the basal law exponent
2641  real, dimension(SZDI_(G),SZDJ_(G)), &
2642  intent(in) :: hmask !< A mask indicating which tracer points are
2643  !! partly or fully covered by an ice-shelf
2644  real, intent(in) :: dens_ratio !< The density of ice divided by the density
2645  !! of seawater, nondimensional
2646  real, dimension(:,:,:,:,:,:), intent(in) :: Phisub !< Quadrature structure weights at subgridscale
2647  !! locations for finite element calculations
2648  real, dimension(SZDIB_(G),SZDJB_(G)), &
2649  intent(inout) :: u_diagonal !< The diagonal elements of the u-velocity
2650  !! matrix from the left-hand side of the solver.
2651  real, dimension(SZDIB_(G),SZDJB_(G)), &
2652  intent(inout) :: v_diagonal !< The diagonal elements of the v-velocity
2653  !! matrix from the left-hand side of the solver.
2654 
2655 
2656 ! returns the diagonal entries of the matrix for a Jacobi preconditioning
2657 
2658  integer :: i, j, is, js, cnt, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq
2659  real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, area, uq, vq, basel
2660  real, dimension(8,4) :: Phi
2661  real, dimension(4) :: X, Y
2662  real, dimension(2) :: xquad
2663  real, dimension(2,2) :: Hcell,Usubcontr,Vsubcontr
2664 
2665 
2666  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2667 
2668  xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2669 
2670 ! X and Y must be passed in the form
2671  ! 3 - 4
2672  ! | |
2673  ! 1 - 2
2674 ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2675 ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
2676 
2677  do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1) then
2678 
2679  dxh = g%US%L_to_m*g%dxT(i,j)
2680  dyh = g%US%L_to_m*g%dyT(i,j)
2681  dxdyh = g%US%L_to_m**2*g%areaT(i,j)
2682 
2683  x(1:2) = g%geoLonBu(i-1:i,j-1)*1000
2684  x(3:4) = g%geoLonBu(i-1:i,j) *1000
2685  y(1:2) = g%geoLatBu(i-1:i,j-1) *1000
2686  y(3:4) = g%geoLatBu(i-1:i,j)*1000
2687 
2688  call bilinear_shape_functions(x, y, phi, area)
2689 
2690  ! X and Y must be passed in the form
2691  ! 3 - 4
2692  ! | |
2693  ! 1 - 2
2694  ! Phi (2*i-1,j) gives d(Phi_i)/dx at quadrature point j
2695  ! Phi (2*i,j) gives d(Phi_i)/dy at quadrature point j
2696 
2697  do iq=1,2 ; do jq=1,2
2698 
2699  do iphi=1,2 ; do jphi=1,2
2700 
2701  if (iq == iphi) then
2702  ilq = 2
2703  else
2704  ilq = 1
2705  endif
2706 
2707  if (jq == jphi) then
2708  jlq = 2
2709  else
2710  jlq = 1
2711  endif
2712 
2713  if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
2714 
2715  ux = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
2716  uy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
2717  vx = 0.
2718  vy = 0.
2719 
2720  u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
2721  .25 * dxdyh * nu(i,j) * ((4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2722  (uy+vy) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2723 
2724  uq = xquad(ilq) * xquad(jlq)
2725 
2726  if (float_cond(i,j) == 0) then
2727  u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
2728  .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
2729  endif
2730 
2731  endif
2732 
2733  if (cs%vmask(i-2+iphi,j-2+jphi) == 1) then
2734 
2735  vx = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
2736  vy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
2737  ux = 0.
2738  uy = 0.
2739 
2740  v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
2741  .25 * dxdyh * nu(i,j) * ((uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2742  (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2743 
2744  vq = xquad(ilq) * xquad(jlq)
2745 
2746  if (float_cond(i,j) == 0) then
2747  v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
2748  .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
2749  endif
2750 
2751  endif
2752  enddo ; enddo
2753  enddo ; enddo
2754  if (float_cond(i,j) == 1) then
2755  usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
2756  hcell(:,:) = h_node(i-1:i,j-1:j)
2757  call cg_diagonal_subgrid_basal(phisub, hcell, dxdyh, basel, dens_ratio, usubcontr, vsubcontr)
2758  do iphi=1,2 ; do jphi=1,2
2759  if (cs%umask(i-2+iphi,j-2+jphi) == 1) then
2760  u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
2761  v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
2762  endif
2763  enddo ; enddo
2764  endif
2765  endif ; enddo ; enddo
2766 

References bilinear_shape_functions(), and cg_diagonal_subgrid_basal().

Referenced by ice_shelf_solve_inner().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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.
[in]yThe y-positions of the vertices of the quadrilateral.

Definition at line 184 of file MOM_ice_shelf_dynamics.F90.

184  real, dimension(4), intent(in) :: X !< The x-positions of the vertices of the quadrilateral.
185  real, dimension(4), intent(in) :: Y !< The y-positions of the vertices of the quadrilateral.
186  real :: quad_area, p2, q2, a2, c2, b2, d2
187 
188 ! X and Y must be passed in the form
189  ! 3 - 4
190  ! | |
191  ! 1 - 2
192 
193  p2 = (x(4)-x(1))**2 + (y(4)-y(1))**2 ; q2 = (x(3)-x(2))**2 + (y(3)-y(2))**2
194  a2 = (x(3)-x(4))**2 + (y(3)-y(4))**2 ; c2 = (x(1)-x(2))**2 + (y(1)-y(2))**2
195  b2 = (x(2)-x(4))**2 + (y(2)-y(4))**2 ; d2 = (x(3)-x(1))**2 + (y(3)-y(1))**2
196  quad_area = .25 * sqrt(4*p2*q2-(b2+d2-a2-c2)**2)
197 

Referenced by bilinear_shape_functions().

Here is the caller graph for this function:

◆ 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 203 of file MOM_ice_shelf_dynamics.F90.

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

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ 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),4), intent(inout)  flux_enter 
)
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]flux_enterThe ice volume flux into the cell

Definition at line 1860 of file MOM_ice_shelf_dynamics.F90.

1860  type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1861  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
1862  !! the ice-shelf state
1863  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
1864  real, dimension(SZDI_(G),SZDJ_(G),4), &
1865  intent(inout) :: flux_enter !< The ice volume flux into the cell
1866  !! through the 4 cell boundaries [Z m2 ~> m3].
1867 
1868  ! in this subroutine we go through the computational cells only and, if they are empty or partial cells,
1869  ! we find the reference thickness and update the shelf mass and partial area fraction and the hmask if necessary
1870 
1871  ! if any cells go from partial to complete, we then must set the thickness, update hmask accordingly,
1872  ! and divide the overflow across the adjacent EMPTY (not partly-covered) cells.
1873  ! (it is highly unlikely there will not be any; in which case this will need to be rethought.)
1874 
1875  ! most likely there will only be one "overflow". if not, though, a pass_var of all relevant variables
1876  ! is done; there will therefore be a loop which, in practice, will hopefully not have to go through
1877  ! many iterations
1878 
1879  ! when 3d advected scalars are introduced, they will be impacted by what is done here
1880 
1881  ! flux_enter(isd:ied,jsd:jed,1:4): if cell is not ice-covered, gives flux of ice into cell from kth boundary
1882  !
1883  ! from left neighbor: flux_enter(:,:,1)
1884  ! from right neighbor: flux_enter(:,:,2)
1885  ! from bottom neighbor: flux_enter(:,:,3)
1886  ! from top neighbor: flux_enter(:,:,4)
1887  !
1888  ! o--- (4) ---o
1889  ! | |
1890  ! (1) (2)
1891  ! | |
1892  ! o--- (3) ---o
1893  !
1894 
1895  integer :: i, j, isc, iec, jsc, jec, n_flux, k, l, iter_count
1896  integer :: i_off, j_off
1897  integer :: iter_flag
1898 
1899  real :: h_reference, dxh, dyh, dxdyh, rho, partial_vol, tot_flux
1900  character(len=160) :: mesg ! The text of an error message
1901  integer, dimension(4) :: mapi, mapj, new_partial
1902 ! real, dimension(size(flux_enter,1),size(flux_enter,2),size(flux_enter,2)) :: flux_enter_replace
1903  real, dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter_replace
1904 
1905  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1906  i_off = g%idg_offset ; j_off = g%jdg_offset
1907  rho = cs%density_ice
1908  iter_count = 0 ; iter_flag = 1
1909 
1910 
1911  mapi(1) = -1 ; mapi(2) = 1 ; mapi(3:4) = 0
1912  mapj(3) = -1 ; mapj(4) = 1 ; mapj(1:2) = 0
1913 
1914  do while (iter_flag == 1)
1915 
1916  iter_flag = 0
1917 
1918  if (iter_count > 0) then
1919  flux_enter(:,:,:) = flux_enter_replace(:,:,:)
1920  endif
1921  flux_enter_replace(:,:,:) = 0.0
1922 
1923  iter_count = iter_count + 1
1924 
1925  ! if iter_count >= 3 then some halo updates need to be done...
1926 
1927  do j=jsc-1,jec+1
1928 
1929  if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1930  ((j+j_off) >= g%domain%njhalo+1)) then
1931 
1932  do i=isc-1,iec+1
1933 
1934  if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1935  ((i+i_off) >= g%domain%nihalo+1)) then
1936  ! first get reference thickness by averaging over cells that are fluxing into this cell
1937  n_flux = 0
1938  h_reference = 0.0
1939  tot_flux = 0.0
1940 
1941  do k=1,2
1942  if (flux_enter(i,j,k) > 0) then
1943  n_flux = n_flux + 1
1944  h_reference = h_reference + iss%h_shelf(i+2*k-3,j)
1945  tot_flux = tot_flux + flux_enter(i,j,k)
1946  flux_enter(i,j,k) = 0.0
1947  endif
1948  enddo
1949 
1950  do k=1,2
1951  if (flux_enter(i,j,k+2) > 0) then
1952  n_flux = n_flux + 1
1953  h_reference = h_reference + iss%h_shelf(i,j+2*k-3)
1954  tot_flux = tot_flux + flux_enter(i,j,k+2)
1955  flux_enter(i,j,k+2) = 0.0
1956  endif
1957  enddo
1958 
1959  if (n_flux > 0) then
1960  dxdyh = g%US%L_to_m**2*g%areaT(i,j)
1961  h_reference = h_reference / real(n_flux)
1962  partial_vol = iss%h_shelf(i,j) * iss%area_shelf_h(i,j) + tot_flux
1963 
1964  if ((partial_vol / dxdyh) == h_reference) then ! cell is exactly covered, no overflow
1965  iss%hmask(i,j) = 1
1966  iss%h_shelf(i,j) = h_reference
1967  iss%area_shelf_h(i,j) = dxdyh
1968  elseif ((partial_vol / dxdyh) < h_reference) then
1969  iss%hmask(i,j) = 2
1970  ! ISS%mass_shelf(i,j) = partial_vol * rho
1971  iss%area_shelf_h(i,j) = partial_vol / h_reference
1972  iss%h_shelf(i,j) = h_reference
1973  else
1974 
1975  iss%hmask(i,j) = 1
1976  iss%area_shelf_h(i,j) = dxdyh
1977  !h_temp(i,j) = h_reference
1978  partial_vol = partial_vol - h_reference * dxdyh
1979 
1980  iter_flag = 1
1981 
1982  n_flux = 0 ; new_partial(:) = 0
1983 
1984  do k=1,2
1985  if (cs%u_face_mask(i-2+k,j) == 2) then
1986  n_flux = n_flux + 1
1987  elseif (iss%hmask(i+2*k-3,j) == 0) then
1988  n_flux = n_flux + 1
1989  new_partial(k) = 1
1990  endif
1991  enddo
1992  do k=1,2
1993  if (cs%v_face_mask(i,j-2+k) == 2) then
1994  n_flux = n_flux + 1
1995  elseif (iss%hmask(i,j+2*k-3) == 0) then
1996  n_flux = n_flux + 1
1997  new_partial(k+2) = 1
1998  endif
1999  enddo
2000 
2001  if (n_flux == 0) then ! there is nowhere to put the extra ice!
2002  iss%h_shelf(i,j) = h_reference + partial_vol / dxdyh
2003  else
2004  iss%h_shelf(i,j) = h_reference
2005 
2006  do k=1,2
2007  if (new_partial(k) == 1) &
2008  flux_enter_replace(i+2*k-3,j,3-k) = partial_vol / real(n_flux)
2009  enddo
2010  do k=1,2 ! ### Combine these two loops?
2011  if (new_partial(k+2) == 1) &
2012  flux_enter_replace(i,j+2*k-3,5-k) = partial_vol / real(n_flux)
2013  enddo
2014  endif
2015 
2016  endif ! Parital_vol test.
2017  endif ! n_flux gt 0 test.
2018 
2019  endif
2020  enddo ! j-loop
2021  endif
2022  enddo
2023 
2024  ! call max_across_PEs(iter_flag)
2025 
2026  enddo ! End of do while(iter_flag) loop
2027 
2028  call max_across_pes(iter_count)
2029 
2030  if (is_root_pe() .and. (iter_count > 1)) then
2031  write(mesg,*) "shelf_advance_front: ", iter_count, " max iterations"
2032  call mom_mesg(mesg, 5)
2033  endif
2034 

References mom_error_handler::is_root_pe(), and mom_error_handler::mom_mesg().

Referenced by ice_shelf_advect().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ slope_limiter()

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

used for flux limiting in advective subroutines Van Leer limiter (source: Wikipedia)

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 166 of file MOM_ice_shelf_dynamics.F90.

166  real, intent(in) :: num !< The numerator of the ratio used in the Van Leer slope limiter
167  real, intent(in) :: denom !< The denominator of the ratio used in the Van Leer slope limiter
168  real :: slope_limiter
169  real :: r
170 
171  if (denom == 0) then
172  slope_limiter = 0
173  elseif (num*denom <= 0) then
174  slope_limiter = 0
175  else
176  r = num/denom
177  slope_limiter = (r+abs(r))/(1+abs(r))
178  endif
179 

Referenced by ice_shelf_advect_temp_x(), ice_shelf_advect_temp_y(), ice_shelf_advect_thickness_x(), and ice_shelf_advect_thickness_y().

Here is the caller graph for this function:

◆ 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]usPointer to a structure containing unit conversion factors
[in]time_steptime step [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 609 of file MOM_ice_shelf_dynamics.F90.

609  type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure
610  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
611  !! the ice-shelf state
612  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
613  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
614  real, intent(in) :: time_step !< time step [s]
615  type(time_type), intent(in) :: Time !< The current model time
616  real, dimension(SZDI_(G),SZDJ_(G)), &
617  optional, intent(in) :: ocean_mass !< If present this is the mass per unit area
618  !! of the ocean [kg m-2].
619  logical, optional, intent(in) :: coupled_grounding !< If true, the grounding line is
620  !! determined by coupled ice-ocean dynamics
621  logical, optional, intent(in) :: must_update_vel !< Always update the ice velocities if true.
622 
623  integer :: iters
624  logical :: update_ice_vel, coupled_GL
625 
626  update_ice_vel = .false.
627  if (present(must_update_vel)) update_ice_vel = must_update_vel
628 
629  coupled_gl = .false.
630  if (present(ocean_mass) .and. present(coupled_grounding)) coupled_gl = coupled_grounding
631 
632  call ice_shelf_advect(cs, iss, g, time_step, time)
633  cs%elapsed_velocity_time = cs%elapsed_velocity_time + time_step
634  if (cs%elapsed_velocity_time >= cs%velocity_update_time_step) update_ice_vel = .true.
635 
636  if (coupled_gl) then
637  call update_od_ffrac(cs, g, us, ocean_mass, update_ice_vel)
638  elseif (update_ice_vel) then
639  call update_od_ffrac_uncoupled(cs, g, iss%h_shelf(:,:))
640  endif
641 
642  if (update_ice_vel) then
643  call ice_shelf_solve_outer(cs, iss, g, us, cs%u_shelf, cs%v_shelf, iters, time)
644  endif
645 
646  call ice_shelf_temp(cs, iss, g, us, time_step, iss%water_flux, time)
647 
648  if (update_ice_vel) then
649  call enable_averaging(cs%elapsed_velocity_time, time, cs%diag)
650  if (cs%id_col_thick > 0) call post_data(cs%id_col_thick, cs%OD_av, cs%diag)
651  if (cs%id_u_shelf > 0) call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
652  if (cs%id_v_shelf > 0) call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
653  if (cs%id_t_shelf > 0) call post_data(cs%id_t_shelf,cs%t_shelf,cs%diag)
654  if (cs%id_float_frac > 0) call post_data(cs%id_float_frac,cs%float_frac,cs%diag)
655  if (cs%id_OD_av >0) call post_data(cs%id_OD_av, cs%OD_av,cs%diag)
656 
657  if (cs%id_u_mask > 0) call post_data(cs%id_u_mask,cs%umask,cs%diag)
658  if (cs%id_v_mask > 0) call post_data(cs%id_v_mask,cs%vmask,cs%diag)
659  if (cs%id_t_mask > 0) call post_data(cs%id_t_mask,cs%tmask,cs%diag)
660 
661  call disable_averaging(cs%diag)
662 
663  cs%elapsed_velocity_time = 0.0
664  endif
665 

References mom_diag_mediator::disable_averaging(), mom_diag_mediator::enable_averaging(), ice_shelf_advect(), ice_shelf_solve_outer(), ice_shelf_temp(), update_od_ffrac(), and update_od_ffrac_uncoupled().

Here is the call graph for this function:

◆ 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( g %isd: g %ied, g %jsd: g %jed), 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]usPointer to a 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 3050 of file MOM_ice_shelf_dynamics.F90.

3050  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
3051  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
3052  type(unit_scale_type), intent(in) :: US !< Pointer to a structure containing unit conversion factors
3053  real, dimension(SZDI_(G),SZDJ_(G)), &
3054  intent(in) :: ocean_mass !< The mass per unit area of the ocean [kg m-2].
3055  logical, intent(in) :: find_avg !< If true, find the average of OD and ffrac, and
3056  !! reset the underlying running sums to 0.
3057 
3058  integer :: isc, iec, jsc, jec, i, j
3059  real :: I_rho_ocean
3060  real :: I_counter
3061 
3062  i_rho_ocean = 1.0 / (us%Z_to_m*cs%density_ocean_avg)
3063 
3064  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3065 
3066  do j=jsc,jec ; do i=isc,iec
3067  cs%OD_rt(i,j) = cs%OD_rt(i,j) + ocean_mass(i,j)*i_rho_ocean
3068  if (ocean_mass(i,j)*i_rho_ocean > cs%thresh_float_col_depth) then
3069  cs%float_frac_rt(i,j) = cs%float_frac_rt(i,j) + 1.0
3070  endif
3071  enddo ; enddo
3072  cs%OD_rt_counter = cs%OD_rt_counter + 1
3073 
3074  if (find_avg) then
3075  i_counter = 1.0 / real(cs%OD_rt_counter)
3076  do j=jsc,jec ; do i=isc,iec
3077  cs%float_frac(i,j) = 1.0 - (cs%float_frac_rt(i,j) * i_counter)
3078  cs%OD_av(i,j) = cs%OD_rt(i,j) * i_counter
3079 
3080  cs%OD_rt(i,j) = 0.0 ; cs%float_frac_rt(i,j) = 0.0
3081  enddo ; enddo
3082 
3083  call pass_var(cs%float_frac, g%domain)
3084  call pass_var(cs%OD_av, g%domain)
3085  endif
3086 

Referenced by update_ice_shelf().

Here is the caller graph for this function:

◆ 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( g %isd: g %ied, g %jsd: g %jed), 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 3090 of file MOM_ice_shelf_dynamics.F90.

3090  type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure
3091  type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf.
3092  real, dimension(SZDI_(G),SZDJ_(G)), &
3093  intent(in) :: h_shelf !< the thickness of the ice shelf [Z ~> m].
3094 
3095  integer :: i, j, iters, isd, ied, jsd, jed
3096  real :: rhoi_rhow, OD
3097 
3098  rhoi_rhow = cs%density_ice / cs%density_ocean_avg
3099  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3100 
3101  do j=jsd,jed
3102  do i=isd,ied
3103  od = g%bathyT(i,j) - rhoi_rhow * h_shelf(i,j)
3104  if (od >= 0) then
3105  ! ice thickness does not take up whole ocean column -> floating
3106  cs%OD_av(i,j) = od
3107  cs%float_frac(i,j) = 0.
3108  else
3109  cs%OD_av(i,j) = 0.
3110  cs%float_frac(i,j) = 1.
3111  endif
3112  enddo
3113  enddo
3114 

Referenced by initialize_ice_shelf_dyn(), and update_ice_shelf().

Here is the caller graph for this function:

◆ 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( g %isd: g %ied, g %jsd: g %jed), intent(in)  hmask,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(out)  umask,
real, dimension( g %isdb: g %iedb, g %jsdb: g %jedb), intent(out)  vmask,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(out)  u_face_mask,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), 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 3253 of file MOM_ice_shelf_dynamics.F90.

3253  type(ice_shelf_dyn_CS),intent(in) :: CS !< A pointer to the ice shelf dynamics control structure
3254  type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf.
3255  real, dimension(SZDI_(G),SZDJ_(G)), &
3256  intent(in) :: hmask !< A mask indicating which tracer points are
3257  !! partly or fully covered by an ice-shelf
3258  real, dimension(SZDIB_(G),SZDJB_(G)), &
3259  intent(out) :: umask !< A coded mask indicating the nature of the
3260  !! zonal flow at the corner point
3261  real, dimension(SZDIB_(G),SZDJB_(G)), &
3262  intent(out) :: vmask !< A coded mask indicating the nature of the
3263  !! meridional flow at the corner point
3264  real, dimension(SZDIB_(G),SZDJ_(G)), &
3265  intent(out) :: u_face_mask !< A coded mask for velocities at the C-grid u-face
3266  real, dimension(SZDI_(G),SZDJB_(G)), &
3267  intent(out) :: v_face_mask !< A coded mask for velocities at the C-grid v-face
3268  ! sets masks for velocity solve
3269  ! ignores the fact that their might be ice-free cells - this only considers the computational boundary
3270 
3271  ! !!!IMPORTANT!!! relies on thickness mask - assumed that this is called after hmask has been updated & halo-updated
3272 
3273  integer :: i, j, k, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
3274  integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec
3275  integer :: i_off, j_off
3276 
3277  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3278  iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
3279  i_off = g%idg_offset ; j_off = g%jdg_offset
3280  isd = g%isd ; jsd = g%jsd
3281  iegq = g%iegB ; jegq = g%jegB
3282  gisc = g%Domain%nihalo ; gjsc = g%Domain%njhalo
3283  giec = g%Domain%niglobal+gisc ; gjec = g%Domain%njglobal+gjsc
3284 
3285  umask(:,:) = 0 ; vmask(:,:) = 0
3286  u_face_mask(:,:) = 0 ; v_face_mask(:,:) = 0
3287 
3288  if (g%symmetric) then
3289  is = isd ; js = jsd
3290  else
3291  is = isd+1 ; js = jsd+1
3292  endif
3293 
3294  do j=js,g%jed
3295  do i=is,g%ied
3296 
3297  if (hmask(i,j) == 1) then
3298 
3299  umask(i-1:i,j-1:j) = 1.
3300  vmask(i-1:i,j-1:j) = 1.
3301 
3302  do k=0,1
3303 
3304  select case (int(cs%u_face_mask_bdry(i-1+k,j)))
3305  case (3)
3306  umask(i-1+k,j-1:j)=3.
3307  vmask(i-1+k,j-1:j)=0.
3308  u_face_mask(i-1+k,j)=3.
3309  case (2)
3310  u_face_mask(i-1+k,j)=2.
3311  case (4)
3312  umask(i-1+k,j-1:j)=0.
3313  vmask(i-1+k,j-1:j)=0.
3314  u_face_mask(i-1+k,j)=4.
3315  case (0)
3316  umask(i-1+k,j-1:j)=0.
3317  vmask(i-1+k,j-1:j)=0.
3318  u_face_mask(i-1+k,j)=0.
3319  case (1) ! stress free x-boundary
3320  umask(i-1+k,j-1:j)=0.
3321  case default
3322  end select
3323  enddo
3324 
3325  do k=0,1
3326 
3327  select case (int(cs%v_face_mask_bdry(i,j-1+k)))
3328  case (3)
3329  vmask(i-1:i,j-1+k)=3.
3330  umask(i-1:i,j-1+k)=0.
3331  v_face_mask(i,j-1+k)=3.
3332  case (2)
3333  v_face_mask(i,j-1+k)=2.
3334  case (4)
3335  umask(i-1:i,j-1+k)=0.
3336  vmask(i-1:i,j-1+k)=0.
3337  v_face_mask(i,j-1+k)=4.
3338  case (0)
3339  umask(i-1:i,j-1+k)=0.
3340  vmask(i-1:i,j-1+k)=0.
3341  u_face_mask(i,j-1+k)=0.
3342  case (1) ! stress free y-boundary
3343  vmask(i-1:i,j-1+k)=0.
3344  case default
3345  end select
3346  enddo
3347 
3348  !if (CS%u_face_mask_bdry(i-1,j) >= 0) then !left boundary
3349  ! u_face_mask(i-1,j) = CS%u_face_mask_bdry(i-1,j)
3350  ! umask(i-1,j-1:j) = 3.
3351  ! vmask(i-1,j-1:j) = 0.
3352  !endif
3353 
3354  !if (j_off+j == gjsc+1) then !bot boundary
3355  ! v_face_mask(i,j-1) = 0.
3356  ! umask (i-1:i,j-1) = 0.
3357  ! vmask (i-1:i,j-1) = 0.
3358  !elseif (j_off+j == gjec) then !top boundary
3359  ! v_face_mask(i,j) = 0.
3360  ! umask (i-1:i,j) = 0.
3361  ! vmask (i-1:i,j) = 0.
3362  !endif
3363 
3364  if (i < g%ied) then
3365  if ((hmask(i+1,j) == 0) &
3366  .OR. (hmask(i+1,j) == 2)) then
3367  !right boundary or adjacent to unfilled cell
3368  u_face_mask(i,j) = 2.
3369  endif
3370  endif
3371 
3372  if (i > g%isd) then
3373  if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2)) then
3374  !adjacent to unfilled cell
3375  u_face_mask(i-1,j) = 2.
3376  endif
3377  endif
3378 
3379  if (j > g%jsd) then
3380  if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2)) then
3381  !adjacent to unfilled cell
3382  v_face_mask(i,j-1) = 2.
3383  endif
3384  endif
3385 
3386  if (j < g%jed) then
3387  if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2)) then
3388  !adjacent to unfilled cell
3389  v_face_mask(i,j) = 2.
3390  endif
3391  endif
3392 
3393 
3394  endif
3395 
3396  enddo
3397  enddo
3398 
3399  ! note: if the grid is nonsymmetric, there is a part that will not be transferred with a halo update
3400  ! so this subroutine must update its own symmetric part of the halo
3401 
3402  call pass_vector(u_face_mask, v_face_mask, g%domain, to_all, cgrid_ne)
3403  call pass_vector(umask, vmask, g%domain, to_all, bgrid_ne)
3404 

References mom_domains::to_all.

Referenced by ice_shelf_advect(), and initialize_ice_shelf_dyn().

Here is the caller graph for this function:
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54