*Use tv%SpV_avg in non-Boussinesq regridding

  Use SpV_avg in both hybgen_unmix and regridding_main to estimate the nominal
ocean bottom depth in thickness units when in fully non-Boussinesq mode.  This
change eliminates certain dependencies on the Boussinesq reference density via
GV%Z_to_H.  Also added a call to calc_derived_thermo in ALE_regrid_accelerated
that is necessary for the specific volume to be updated, along with tests for
the validity of the specific volume with the 1-point halos that are currently
used with ALE regridding and remapping.

  Also, use the total thickness place of the nominal depth in build_grid_rho and
build_grid_sigma when in fully non-Boussinesq mode.  This is mathematically
equivalent but changes answers at roundoff.

  All answers are bitwise identical in Boussinesq mode, but ALE-based solutions
change in fully non-Boussinesq mode with this commit.

src/ALE/MOM_ALE.F90

@@ -20,7 +20,7 @@ module MOM_ALE
 use MOM_hybgen_unmix,     only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_CS
 use MOM_hybgen_regrid,    only : hybgen_regrid_CS
 use MOM_file_parser,      only : get_param, param_file_type, log_param
-use MOM_interface_heights,only : find_eta
+use MOM_interface_heights,only : find_eta, calc_derived_thermo
 use MOM_open_boundary,    only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
 use MOM_open_boundary,    only : OBC_DIRECTION_N, OBC_DIRECTION_S
 use MOM_regridding,       only : initialize_regridding, regridding_main, end_regridding
@@ -659,6 +659,9 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d
     ! generate new grid
     if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local)
 
+    ! Update the layer specific volumes if necessary
+    if (allocated(tv_local%SpV_avg)) call calc_derived_thermo(tv_local, h, G, GV, US, halo=1)
+
     call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface)
     dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:)
 

src/ALE/MOM_hybgen_unmix.F90

@@ -9,6 +9,7 @@ module MOM_hybgen_unmix
 use MOM_file_parser,     only : get_param, param_file_type, log_param
 use MOM_hybgen_regrid,   only : hybgen_column_init
 use MOM_hybgen_regrid,   only : hybgen_regrid_CS, get_hybgen_regrid_params
+use MOM_interface_heights, only : calc_derived_thermo
 use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chkinv
 use MOM_unit_scaling,    only : unit_scale_type
 use MOM_variables,       only : ocean_grid_type, thermo_var_ptrs
@@ -146,7 +147,8 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
   real :: p_col(GV%ke)      ! A column of reference pressures [R L2 T-2 ~> Pa]
   real :: tracer(GV%ke,max(ntr,1)) ! Columns of each tracer [Conc]
   real :: h_tot             ! Total thickness of the water column [H ~> m or kg m-2]
-  real :: nominalDepth      ! Depth of ocean bottom (positive downward) [H ~> m or kg m-2]
+  real :: dz_tot            ! Vertical distance between the top and bottom of the water column [Z ~> m]
+  real :: nominalDepth      ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2]
   real :: h_thin            ! A negligibly small thickness to identify essentially
                             ! vanished layers [H ~> m or kg m-2]
   real :: dilate            ! A factor by which to dilate the target positions from z to z* [nondim]
@@ -169,6 +171,15 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
   h_thin = 1e-6*GV%m_to_H
   debug_conservation = .false. !  Set this to true for debugging
 
+  if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
+    if (tv%valid_SpV_halo < 0) then
+      mesg = "invalid values of SpV_avg."
+    else
+      mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
+    endif
+    call MOM_error(FATAL, "hybgen_unmix called in fully non-Boussinesq mode with "//trim(mesg))
+  endif
+
   p_col(:) = CS%ref_pressure
 
   do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then
@@ -203,13 +214,27 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
     endif
 
     ! The following block of code is used to trigger z* stretching of the targets heights.
-    nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H
-    if (h_tot <= CS%min_dilate*nominalDepth) then
-      dilate = CS%min_dilate
-    elseif (h_tot >= CS%max_dilate*nominalDepth) then
-      dilate = CS%max_dilate
+    if (allocated(tv%SpV_avg)) then  ! This is the fully non-Boussiesq version
+      dz_tot = 0.0
+      do k=1,nk
+        dz_tot = dz_tot + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h_col(k)
+      enddo
+      if (dz_tot <= CS%min_dilate*(G%bathyT(i,j)+G%Z_ref)) then
+        dilate = CS%min_dilate
+      elseif (dz_tot >= CS%max_dilate*(G%bathyT(i,j)+G%Z_ref)) then
+        dilate = CS%max_dilate
+      else
+        dilate = dz_tot / (G%bathyT(i,j)+G%Z_ref)
+      endif
     else
-      dilate = h_tot / nominalDepth
+      nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H
+      if (h_tot <= CS%min_dilate*nominalDepth) then
+        dilate = CS%min_dilate
+      elseif (h_tot >= CS%max_dilate*nominalDepth) then
+        dilate = CS%max_dilate
+      else
+        dilate = h_tot / nominalDepth
+      endif
     endif
 
     terrain_following = (h_tot < dilate*CS%dpns) .and. (CS%dpns >= CS%dsns)
@@ -268,6 +293,9 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
     endif
   endif ; enddo ; enddo !i & j.
 
+  ! Update the layer properties
+  if (allocated(tv%SpV_avg)) call calc_derived_thermo(tv, h, G, GV, US, halo=1)
+
 end subroutine hybgen_unmix
 
 

src/ALE/MOM_regridding.F90

@@ -814,20 +814,47 @@ subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, &
 
   ! Local variables
   real :: nom_depth_H(SZI_(G),SZJ_(G))  !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2]
+  real :: tot_h(SZI_(G),SZJ_(G))  !< The total thickness of the water column [H ~> m or kg m-2]
+  real :: tot_dz(SZI_(G),SZJ_(G)) !< The total distance between the top and bottom of the water column [Z ~> m]
   real :: Z_to_H  ! A conversion factor used by some routines to convert coordinate
                   ! parameters to depth units [H Z-1 ~> nondim or kg m-3]
   real :: trickGnuCompiler
-  integer :: i, j
+  character(len=128) :: mesg    ! A string for error messages
+  integer :: i, j, k
 
   if (present(PCM_cell)) PCM_cell(:,:,:) = .false.
 
   Z_to_H = US%Z_to_m * GV%m_to_H  ! Often this is equivalent to GV%Z_to_H.
-  do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
-    nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * Z_to_H
-    ! Consider using the following instead:
-    ! nom_depth_H(i,j) = max( (G%bathyT(i,j)+G%Z_ref) * Z_to_H , CS%min_nom_depth )
-    ! if (G%mask2dT(i,j)==0.) nom_depth_H(i,j) = 0.0
-  enddo ; enddo
+
+  if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
+    if (tv%valid_SpV_halo < 0) then
+      mesg = "invalid values of SpV_avg."
+    else
+      mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
+    endif
+    call MOM_error(FATAL, "Regridding_main called in fully non-Boussinesq mode with "//trim(mesg))
+  endif
+
+  if (allocated(tv%SpV_avg)) then  ! This is the fully non-Boussinesq case
+    do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      tot_h(i,j) = 0.0 ; tot_dz(i,j) = 0.0
+    enddo ; enddo
+    do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      tot_h(i,j) = tot_h(i,j) + h(i,j,k)
+      tot_dz(i,j) = tot_dz(i,j) + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h(i,j,k)
+    enddo ; enddo ; enddo
+    do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      if ((tot_dz(i,j) > 0.0) .and. (G%bathyT(i,j)+G%Z_ref > 0.0)) then
+        nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * (tot_h(i,j) / tot_dz(i,j))
+      else
+        nom_depth_H(i,j) = 0.0
+      endif
+    enddo ; enddo
+  else
+    do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      nom_depth_H(i,j) = max((G%bathyT(i,j)+G%Z_ref) * Z_to_H, 0.0)
+    enddo ; enddo
+  endif
 
   select case ( CS%regridding_scheme )
 
@@ -1308,12 +1335,12 @@ subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface )
       ! In sigma coordinates, the bathymetric depth is only used as an arbitrary offset that
       ! cancels out when determining coordinate motion, so referencing the column postions to
       ! the surface is perfectly acceptable, but for preservation of previous answers the
-      ! referencing is done relative to the bottom when in Boussinesq mode.
-      ! if (GV%Boussinesq) then
+      ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
+      if (GV%Boussinesq .or. GV%semi_Boussinesq) then
         nominalDepth = nom_depth_H(i,j)
-      ! else
-      !   nominalDepth = totalThickness
-      ! endif
+      else
+        nominalDepth = totalThickness
+      endif
 
       call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew)
 
@@ -1436,12 +1463,12 @@ subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS,
       ! In rho coordinates, the bathymetric depth is only used as an arbitrary offset that
       ! cancels out when determining coordinate motion, so referencing the column postions to
       ! the surface is perfectly acceptable, but for preservation of previous answers the
-      ! referencing is done relative to the bottom when in Boussinesq mode.
-      ! if (GV%Boussinesq) then
+      ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
+      if (GV%Boussinesq .or. GV%semi_Boussinesq) then
         nominalDepth = nom_depth_H(i,j)
-      ! else
-      !   nominalDepth = totalThickness
-      ! endif
+      else
+        nominalDepth = totalThickness
+      endif
 
       ! Determine absolute interface positions
       zOld(nz+1) = - nominalDepth