Allow ALE boundary extrapolation behaviour to differ at initialisatio…

…n and in model run

Adds INIT_BOUNDARY_EXTRAP parameter and function, so that REMAP_BOUNDARY_EXTRAP can just be used in the initialisation and not in the model dynamics, as required for an ice shelf in ALE mode.

Background:
When an ice shelf is initialised in ALE mode, MOM_state_initialization first initialises the ocean thickness and T/S with topography, but without an ice shelf. It then calls the function trim_for_ice, which takes in input of the ice shelf pressure, and in turn calls cut_off_column_top which truncates the columns (and adds vanishing layers) to account for the pressure of the ice depressing the free surface. The output of this is a grid with the ice shelf boundary treated in a ALE z-coordinate fashion with vanishing layers at the surface, and @adcroft's upcoming changes will make T and S consistent to the grid. However, after adding the ice shelf, we still need to regrid and remap onto our desired ALE coordinate, for example sigma coordinates. This can be done with regrid_accelerate in MOM_state_initialization or by the block controlled by REMAP_AFTER_INITIALIZATION in MOM.F90

We want the initialisation to use boundary extrapolation in the remapping so that the boundary cells preserve the desired initialisation profile when remapped to their target coordinate. This is controlled by REMAP_BOUNDARY_EXTRAP, except this runtime parameter modifies the control structure for the dynamic part of the model run too, which is bad because it means extrema are not preserved in the ALE remapping and can lead to non-sensible values of tracers.

Changes:

* Added runtime parameter INIT_BOUNDARY_EXTRAP and set the ALE control structure to use that instead of REMAP_BOUNDARY_EXTRAP at first. The default of INIT_BOUNDARY_EXTRAP is REMAP_BOUNDARY_EXTRAP to preserve answers (though it is not recommended to use REMAP_BOUNDARY_EXTRAP during model dynamics...)
 * After ALE regrid/remap in MOM.F90 (final initialisation step using ALE), reset the boundary extrapolation flag in the ALE control structure to go back to REMAP_BOUNDARY_EXTRAP.

This PR combined with @adcroft's upcoming changes should make the initialisation of the ocean thickness and TS in an ALE ice shelf cavity perfect, when the new flag INIT_BOUNDARY_EXTRAP = True and old flag REMAP_BOUNDARY_EXTRAP = False.

src/ALE/MOM_ALE.F90

@@ -40,7 +40,7 @@ module MOM_ALE
 use MOM_remapping,        only : remapping_core_h, remapping_core_w
 use MOM_remapping,        only : remappingSchemesDoc, remappingDefaultScheme
 use MOM_remapping,        only : interpolate_column, reintegrate_column
-use MOM_remapping,        only : remapping_CS, dzFromH1H2
+use MOM_remapping,        only : remapping_CS, dzFromH1H2, remapping_set_param
 use MOM_string_functions, only : uppercase, extractWord, extract_integer
 use MOM_tracer_registry,  only : tracer_registry_type, tracer_type, MOM_tracer_chkinv
 use MOM_unit_scaling,     only : unit_scale_type
@@ -147,6 +147,7 @@ module MOM_ALE
 public pre_ALE_adjustments
 public ALE_remap_init_conds
 public ALE_register_diags
+public ALE_set_extrap_boundaries
 
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
@@ -176,6 +177,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   logical           :: force_bounds_in_subcell
   logical           :: local_logical
   logical           :: remap_boundary_extrap
+  logical           :: init_boundary_extrap
   type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
                                                          ! for sharing parameters.
 
@@ -225,6 +227,10 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
                  "If true, values at the interfaces of boundary cells are "//&
                  "extrapolated instead of piecewise constant", default=.false.)
+  call get_param(param_file, mdl, "INIT_BOUNDARY_EXTRAP", init_boundary_extrap, &
+                 "If true, values at the interfaces of boundary cells are "//&
+                 "extrapolated instead of piecewise constant during initialization."//&
+                 "Defaults to REMAP_BOUNDARY_EXTRAP.", default=remap_boundary_extrap)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)
@@ -237,13 +243,13 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
 
   call initialize_remapping( CS%remapCS, string, &
-                             boundary_extrapolation=remap_boundary_extrap, &
+                             boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
                              answer_date=CS%answer_date)
   call initialize_remapping( CS%vel_remapCS, vel_string, &
-                             boundary_extrapolation=remap_boundary_extrap, &
+                             boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
@@ -308,6 +314,18 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   if (CS%show_call_tree) call callTree_leave("ALE_init()")
 end subroutine ALE_init
 
+!> Sets the boundary extrapolation set for the remapping type.
+subroutine ALE_set_extrap_boundaries( param_file, CS)
+  type(param_file_type),   intent(in) :: param_file !< Parameter file
+  type(ALE_CS),            pointer    :: CS         !< Module control structure
+
+  logical :: remap_boundary_extrap
+  call get_param(param_file, "MOM_ALE", "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
+                 "If true, values at the interfaces of boundary cells are "//&
+                 "extrapolated instead of piecewise constant", default=.false.)
+  call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap)
+end subroutine ALE_set_extrap_boundaries
+
 !> Initialize diagnostics for the ALE module.
 subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   type(time_type),target,     intent(in)  :: Time  !< Time structure

src/core/MOM.F90

@@ -56,6 +56,7 @@ module MOM
 use MOM_ALE,                   only : ALE_remap_tracers, ALE_remap_velocities
 use MOM_ALE,                   only : ALE_remap_set_h_vel, ALE_remap_set_h_vel_via_dz
 use MOM_ALE,                   only : ALE_update_regrid_weights, pre_ALE_diagnostics, ALE_register_diags
+use MOM_ALE,                   only : ALE_set_extrap_boundaries
 use MOM_ALE_sponge,            only : rotate_ALE_sponge, update_ALE_sponge_field
 use MOM_barotropic,            only : Barotropic_CS
 use MOM_boundary_update,       only : call_OBC_register, OBC_register_end, update_OBC_CS
@@ -3120,8 +3121,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
       endif
     endif
   endif
-  if ( CS%use_ALE_algorithm ) call ALE_updateVerticalGridType( CS%ALE_CSp, GV )
-
+  if ( CS%use_ALE_algorithm ) then
+   call ALE_set_extrap_boundaries (param_file, CS%ALE_CSp)
+   call callTree_waypoint("returned from ALE_init() (initialize_MOM)")
+   call ALE_updateVerticalGridType( CS%ALE_CSp, GV )
+  endif
   ! The basic state variables have now been fully initialized, so update their halos and
   ! calculate any derived thermodynmics quantities.