Ice sheet: consistency under dimensional rescaling and rotation

[alex-huth]

This PR guarantees consistency of the ice sheet code under rotated/flipped input and grids, and under dimensional rescaling to the +/-140 power. Bitwise identical results were achieved under rotation/rescaling using standard 2D test cases such as MISMIP+. Major changes include the following:

• Specified order of operations throughout, and modified the implementation of FEM integration, h-point-to-node operations, etc for consistency under rotation (see MOM_ice_shelf_dynamics.F90 subroutines calc_shelf_driving_stress, CG_action, CG_action_subgrid_basal, matrix_diagonal, shelf_advance_front, and interpolate_h_to_b).
• Added an ice-shelf version of 'first_direction' and 'alterate_first_direction', allowing users to control the order in which the x- and y-direction ice-shelf advection calls are made. Required for consistency under rotation.
• Dimensional rescaling fixes throughout MOM_ice_shelf_initialize.F90, and within MOM_ice_shelf_dynamics.F90 (see subroutines initialize_ice_shelf_dyn, ice_shelf_solve_outer, and calc_shelf_taub)
• Rotation/rescaling testing revealed two additional bugs that were subsequently fixed:
  (1) An index error in the conjugate gradient algorithm (subroutine CG_action)
  (2) Discretization error for the east/north fluxes in subroutines ice_shelf_advect_thickness_x and ice_shelf_advect_thickness_y.

The PR also includes several simple changes for code readability and computational efficiency:

• Saved Phi, PhiC, and Phisub in the ice shelf dynamics control structure at the start of each run, rather than reallocating and recalculating these static fields each time they are used. Reshaped the Phisub array for computational efficiency over loops.
• Modified the sub-cell grounding scheme so that it is only called for partially-grounded cells that contain the grounding line (i.e. where float_cond(i,j)==1). Added float_cond to the ice dynamics structure for output.
• Simplified the option to calculate ice viscosity at 4 quadrature points per cell rather than at cell centers, by adding parameter NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS and reshaping CS%ice_visc accordingly.
• Style changes and removal of unused code (e.g. removed subroutine apply_boundary_values and field thickness_bdry_val)

[codecov]

Codecov Report

Attention: 13 lines in your changes are missing coverage. Please review.

Comparison is base (5c14668) 37.46% compared to head (8e82f80) 37.44%.

Files	Patch %	Lines
src/ice_shelf/MOM_ice_shelf_initialize.F90	0.00%	13 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff              @@
##           dev/gfdl     #527      +/-   ##
============================================
- Coverage     37.46%   37.44%   -0.02%     
============================================
  Files           270      270              
  Lines         79636    79674      +38     
  Branches      14825    14828       +3     
============================================
+ Hits          29835    29836       +1     
- Misses        44252    44290      +38     
+ Partials       5549     5548       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Hallberg-NOAA]

I have carefully examined all of these changes, and they all look good to me, apart from one omitted change in the units in a comment on line 3074 (see the specific comment).

This will be ready to go once that very minor change is made.

[marshallward]

Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/21727 ✔️