Merge remote-tracking branch 'gfdl/dev/gfdl' into esmg_work

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.
 

src/parameterizations/vertical/MOM_CVMix_KPP.F90

@@ -119,6 +119,10 @@ module MOM_CVMix_KPP
                                        !! enhancement with KPP [Z ~> m]
   logical :: STOKES_MIXING             !< Flag if model is mixing down Stokes gradient
                                        !! This is relevant for which current to use in RiB
+  integer :: answer_date               !< The vintage of the order of arithmetic in the CVMix KPP
+                                       !! calculations.  Values below 20240501 recover the answers
+                                       !! from early in 2024, while higher values use expressions
+                                       !! that have been refactored for rotational symmetry.
 
   !> CVMix parameters
   type(CVMix_kpp_params_type), pointer :: KPP_params => NULL()
@@ -199,6 +203,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
   character(len=20) :: langmuir_mixing_opt = 'NONE' !< Langmuir mixing option to be passed to CVMix, e.g., LWF16
   character(len=20) :: langmuir_entrainment_opt = 'NONE' !< Langmuir entrainment option to be
                                        !! passed to CVMix, e.g., LWF16
+  integer :: default_answer_date       ! The default setting for the various ANSWER_DATE flags.
   logical :: CS_IS_ONE=.false.         !< Logical for setting Cs based on Non-local
   logical :: lnoDGat1=.false.          !< True => G'(1) = 0 (shape function)
                                        !! False => compute G'(1) as in LMD94
@@ -217,6 +222,10 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
   if (.not. KPP_init) return
   allocate(CS)
 
+  call get_param(paramFile, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
+                 "This sets the default value for the various _ANSWER_DATE parameters.", &
+                 default=99991231, do_not_log=.true.)
+
   call openParameterBlock(paramFile,'KPP')
   call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode,           &
                  'If True, puts KPP into a passive-diagnostic mode.', &
@@ -471,6 +480,12 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
                    units="m", default=1.0, scale=US%m_to_Z)
   endif
 
+  call get_param(paramFile, mdl, "ANSWER_DATE", CS%answer_date, &
+                 "The vintage of the order of arithmetic in the CVMix KPP calculations.  Values "//&
+                 "below 20240501 recover the answers from early in 2024, while higher values "//&
+                 "use expressions that have been refactored for rotational symmetry.", &
+                 default=20240101) !### Change to: default=default_answer_date)
+
   call closeParameterBlock(paramFile)
 
   call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.)
@@ -1390,11 +1405,17 @@ subroutine KPP_smooth_BLD(CS, G, GV, US, dz)
         wn = 0.125 * G%mask2dT(i,j+1)
         wc = 1.0 - (ww+we+wn+ws)
 
-        CS%OBLdepth(i,j) =  wc * OBLdepth_prev(i,j)   &
-                          + ww * OBLdepth_prev(i-1,j) &
-                          + we * OBLdepth_prev(i+1,j) &
-                          + ws * OBLdepth_prev(i,j-1) &
-                          + wn * OBLdepth_prev(i,j+1)
+        if (CS%answer_date < 20240501) then
+          CS%OBLdepth(i,j) =  wc * OBLdepth_prev(i,j)   &
+                            + ww * OBLdepth_prev(i-1,j) &
+                            + we * OBLdepth_prev(i+1,j) &
+                            + ws * OBLdepth_prev(i,j-1) &
+                            + wn * OBLdepth_prev(i,j+1)
+        else
+          CS%OBLdepth(i,j) =  wc * OBLdepth_prev(i,j) &
+                            + ((ww * OBLdepth_prev(i-1,j) + we * OBLdepth_prev(i+1,j)) &
+                             + (ws * OBLdepth_prev(i,j-1) + wn * OBLdepth_prev(i,j+1)))
+        endif
 
         ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing.
         if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j), OBLdepth_prev(i,j))