+Add tv%valid_SpV_halo to debug non-Boussinesq mode

  Added the integer valid_SpV_halo to the thermo_var_ptrs type to indicate
whether the SpV_array has been updated and its valid halo size, to facilitate
error detection and debugging in non-Boussinesq mode.  Tv%valid_SpV_halo is set
to the halo size in calc_derived_thermo or after a halo update is done to
tv%SpV_avg, and it is set to a negative value right after calls that change
temperatures and salinities (such as by ALE remapping) unless there is a call to
calc_derived_thermo.  Tests for the validity of tv%SpV_avg are added to the
routines behind thickness_to_dz, with fatal errors issued if invalid arrays
would be used, but more tests could perhaps be used in any parameterization
routines where tv%SpV_avg is used directly.  Handling the updates to tv%SpV_avg
this way helps to avoid unnecessary calls to calc_derived_thermo, which in turn
has equation of state calls that can be expensive, while also providing
essential verification of new code related to the non-Boussinesq code.  These
tests can probably be commented out or removed for efficiency once there is a
full suite of regression tests for the fully non-Boussinesq mode of MOM6.  In
addition, a new optional debug argument was added to calc_derived_thermo which
can be used to triggers checksums for the variables used to calculate
tv%SpV_avg.  One call to calc_derived_thermo was also added just before the
initialization call to ALE_regrid that will be needed with the next commit, but
does not change answers yet.  All answers are bitwise identical, but there is a
new element in a transparent and widely used type and a new optional argument to
a public interface.

src/ALE/MOM_ALE.F90

@@ -532,6 +532,7 @@ subroutine ALE_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug,
 
   ! Remap all variables from old grid h onto new grid h_new
   call ALE_remap_tracers(CS, G, GV, h, h_new, Reg, debug=CS%show_call_tree)
+  if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
   if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_inputs)")
 
   ! Reintegrate mass transports from Zstar to the offline vertical coordinate
@@ -571,6 +572,8 @@ subroutine ALE_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug,
     h(i,j,k) = h_new(i,j,k)
   enddo ; enddo ; enddo
 
+  if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
+
   if (CS%show_call_tree) call callTree_leave("ALE_offline_inputs()")
 end subroutine ALE_offline_inputs
 
@@ -674,6 +677,9 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d
 
   ! save total dzregrid for diags if needed?
   if (present(dzRegrid)) dzRegrid(:,:,:) = dzIntTotal(:,:,:)
+
+  if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
+
 end subroutine ALE_regrid_accelerated
 
 !> This routine takes care of remapping all tracer variables between the old and the

src/core/MOM.F90

@@ -671,8 +671,9 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     dt_therm = dt ; ntstep = 1
     if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf
     CS%tv%p_surf => NULL()
-    if (associated(fluxes%p_surf)) then
-      if (CS%use_p_surf_in_EOS) CS%tv%p_surf => fluxes%p_surf
+    if (CS%use_p_surf_in_EOS .and. associated(fluxes%p_surf)) then
+      CS%tv%p_surf => fluxes%p_surf
+      if (allocated(CS%tv%SpV_avg)) call pass_var(fluxes%p_surf, G%Domain, clock=id_clock_pass)
     endif
     if (CS%UseWaves) call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass)
   endif
@@ -1110,6 +1111,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     endif
 
     if (CS%interface_filter) then
+      if (allocated(CS%tv%SpV_avg)) call pass_var(CS%tv%SpV_avg, G%Domain, clock=id_clock_pass)
+      CS%tv%valid_SpV_halo = min(G%Domain%nihalo, G%Domain%njhalo)
       call cpu_clock_begin(id_clock_int_filter)
       call interface_filter(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, &
                             CS%CDp, CS%interface_filter_CSp)
@@ -1245,6 +1248,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     endif
 
     if (CS%interface_filter) then
+      if (allocated(CS%tv%SpV_avg)) call pass_var(CS%tv%SpV_avg, G%Domain, clock=id_clock_pass)
+      CS%tv%valid_SpV_halo = min(G%Domain%nihalo, G%Domain%njhalo)
       call cpu_clock_begin(id_clock_int_filter)
       call interface_filter(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, &
                             CS%CDp, CS%interface_filter_CSp)
@@ -1392,6 +1397,8 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
 
   if (associated(CS%tv%T)) then
     call extract_diabatic_member(CS%diabatic_CSp, diabatic_halo=halo_sz)
+    ! The bottom boundary layer calculation may need halo values of SpV_avg, including the corners.
+    if (allocated(CS%tv%SpV_avg)) halo_sz = max(halo_sz, 1)
     if (halo_sz > 0) then
       call create_group_pass(pass_T_S, CS%tv%T, G%Domain, To_All, halo=halo_sz)
       call create_group_pass(pass_T_S, CS%tv%S, G%Domain, To_All, halo=halo_sz)
@@ -1407,7 +1414,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
 
     ! Update derived thermodynamic quantities.
     if (allocated(CS%tv%SpV_avg)) then
-      call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz)
+      call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz, debug=CS%debug)
     endif
   endif
 
@@ -1559,6 +1566,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
       ! Remap all variables from the old grid h onto the new grid h_new
       call ALE_remap_tracers(CS%ALE_CSp, G, GV, h, h_new, CS%tracer_Reg, showCallTree, dtdia, PCM_cell)
       call ALE_remap_velocities(CS%ALE_CSp, G, GV, h, h_new, u, v, CS%OBC, dzRegrid, showCallTree, dtdia)
+      if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
 
       if (CS%remap_aux_vars) then
         if (CS%split) &
@@ -1591,7 +1599,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
 
     ! Update derived thermodynamic quantities.
     if (allocated(tv%SpV_avg)) then
-      call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil)
+      call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil, debug=CS%debug)
     endif
 
     if (CS%debug .and. CS%use_ALE_algorithm) then
@@ -1647,7 +1655,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
 
       ! Update derived thermodynamic quantities.
       if (allocated(tv%SpV_avg)) then
-        call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil)
+        call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil, debug=CS%debug)
       endif
     endif
 
@@ -1820,6 +1828,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
       ! are used are intended to ensure that in the case where transports don't quite conserve,
       ! the offline layer thicknesses do not drift too far away from the online model.
       call ALE_remap_tracers(CS%ALE_CSp, G, GV, CS%h, h_new, CS%tracer_Reg, debug=CS%debug)
+      if (allocated(CS%tv%SpV_avg)) CS%tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
 
       ! Update the tracer grid.
       do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
@@ -2847,6 +2856,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   ! Allocate any derived equation of state fields.
   if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
     allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0)
+    CS%tv%valid_SpV_halo = -1  ! This array does not yet have any valid data.
   endif
 
   if (use_ice_shelf .and. CS%debug) then
@@ -2895,6 +2905,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     endif
     call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
     call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h )
+    if (allocated(CS%tv%SpV_avg)) call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=1)
     call pre_ALE_adjustments(G, GV, US, CS%h, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%u, CS%v)
 
     call callTree_waypoint("Calling ALE_regrid() to remap initial conditions (initialize_MOM)")
@@ -2911,6 +2922,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     ! Remap all variables from the old grid h onto the new grid h_new
     call ALE_remap_tracers(CS%ALE_CSp, G, GV, CS%h, h_new, CS%tracer_Reg, CS%debug, PCM_cell=PCM_cell)
     call ALE_remap_velocities(CS%ALE_CSp, G, GV, CS%h, h_new, CS%u, CS%v, CS%OBC, dzRegrid, debug=CS%debug)
+    if (allocated(CS%tv%SpV_avg)) CS%tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
 
     ! Replace the old grid with new one.  All remapping must be done at this point.
     !$OMP parallel do default(shared)
@@ -3137,7 +3149,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   ! Update derived thermodynamic quantities.
   if (allocated(CS%tv%SpV_avg)) then
-    call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil)
+    call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil, debug=CS%debug)
   endif
 
   if (associated(CS%visc%Kv_shear)) &

src/core/MOM_interface_heights.F90

@@ -4,13 +4,14 @@ module MOM_interface_heights
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol
+use MOM_debugging,     only : hchksum
 use MOM_error_handler, only : MOM_error, FATAL
-use MOM_EOS, only : calculate_density, EOS_type, EOS_domain
-use MOM_file_parser, only : log_version
-use MOM_grid, only : ocean_grid_type
-use MOM_unit_scaling, only : unit_scale_type
-use MOM_variables, only : thermo_var_ptrs
-use MOM_verticalGrid, only : verticalGrid_type
+use MOM_EOS,           only : calculate_density, average_specific_vol, EOS_type, EOS_domain
+use MOM_file_parser,   only : log_version
+use MOM_grid,          only : ocean_grid_type
+use MOM_unit_scaling,  only : unit_scale_type
+use MOM_variables,     only : thermo_var_ptrs
+use MOM_verticalGrid,  only : verticalGrid_type
 
 implicit none ; private
 
@@ -262,7 +263,7 @@ end subroutine find_eta_2d
 
 
 !> Calculate derived thermodynamic quantities for re-use later.
-subroutine calc_derived_thermo(tv, h, G, GV, US, halo)
+subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
   type(ocean_grid_type),   intent(in)    :: G  !< The ocean's grid structure
   type(verticalGrid_type), intent(in)    :: GV !< The ocean's vertical grid structure
   type(unit_scale_type),   intent(in)    :: US !< A dimensional unit scaling type
@@ -271,13 +272,16 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo)
                                                !! which will be set here.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                            intent(in)    :: h  !< Layer thicknesses [H ~> m or kg m-2].
-  integer,         optional, intent(in)  :: halo !< Width of halo within which to
+  integer,       optional, intent(in)    :: halo !< Width of halo within which to
                                                !! calculate thicknesses
+  logical,       optional, intent(in)    :: debug !< If present and true, write debugging checksums
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G)) :: p_t  ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa]
   real, dimension(SZI_(G),SZJ_(G)) :: dp   ! Pressure change across a layer [R L2 T-2 ~> Pa]
+  logical :: do_debug  ! If true, write checksums for debugging.
   integer :: i, j, k, is, ie, js, je, halos, nz
 
+  do_debug = .false. ; if (present(debug)) do_debug = debug
   halos = 0 ; if (present(halo)) halos = max(0,halo)
   is = G%isc-halos ; ie = G%iec+halos ; js = G%jsc-halos ; je = G%jec+halos ; nz = GV%ke
 
@@ -296,6 +300,15 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo)
         p_t(i,j) = p_t(i,j) + dp(i,j)
       enddo ; enddo ; endif
     enddo
+    tv%valid_SpV_halo = halos
+
+    if (do_debug) then
+      call hchksum(h, "derived_thermo h", G%HI, haloshift=halos, scale=GV%H_to_MKS)
+      if (associated(tv%p_surf)) call hchksum(tv%p_surf, "derived_thermo p_surf", G%HI, &
+                                              haloshift=halos, scale=US%RL2_T2_to_Pa)
+      call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, scale=US%C_to_degC)
+      call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, scale=US%S_to_ppt)
+    endif
   endif
 
 end subroutine calc_derived_thermo
@@ -493,12 +506,23 @@ subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size)
   integer,       optional, intent(in)    :: halo_size !< Width of halo within which to
                                                !! calculate thicknesses
   ! Local variables
+  character(len=128) :: mesg    ! A string for error messages
   integer :: i, j, k, is, ie, js, je, halo, nz
 
   halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
   is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke
 
   if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg))  then
+    if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then
+      if (tv%valid_SpV_halo < 0) then
+        mesg = "invalid values of SpV_avg."
+      else
+        write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
+                     tv%valid_SpV_halo, halo
+      endif
+      call MOM_error(FATAL, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg))
+    endif
+
     do k=1,nz ; do j=js,je ; do i=is,ie
       dz(i,j,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
     enddo ; enddo ; enddo
@@ -529,12 +553,23 @@ subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size)
   integer,       optional, intent(in)    :: halo_size !< Width of halo within which to
                                                !! calculate thicknesses
   ! Local variables
+  character(len=128) :: mesg    ! A string for error messages
   integer :: i, k, is, ie, halo, nz
 
   halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
   is = G%isc-halo ; ie = G%iec+halo ; nz = GV%ke
 
   if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg))  then
+    if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then
+      if (tv%valid_SpV_halo < 0) then
+        mesg = "invalid values of SpV_avg."
+      else
+        write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
+                     tv%valid_SpV_halo, halo
+      endif
+      call MOM_error(FATAL, "thickness_to_dz called in fully non-Boussinesq mode with "//trim(mesg))
+    endif
+
     do k=1,nz ; do i=is,ie
       dz(i,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
     enddo ; enddo

src/core/MOM_variables.F90

@@ -95,6 +95,8 @@ module MOM_variables
   real :: min_salinity   !< The minimum value of salinity when BOUND_SALINITY=True [S ~> ppt].
   real, allocatable, dimension(:,:,:) :: SpV_avg
                          !< The layer averaged in situ specific volume [R-1 ~> m3 kg-1].
+  integer :: valid_SpV_halo = -1 !< If positive, the valid halo size for SpV_avg, or if negative
+                         !! SpV_avg is not currently set.
 
   ! These arrays are accumulated fluxes for communication with other components.
   real, dimension(:,:), pointer :: frazil => NULL()

src/tracer/MOM_offline_main.F90

@@ -364,6 +364,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C
       ! Remap all variables from the old grid h_new onto the new grid h_post_remap
       call ALE_remap_tracers(CS%ALE_CSp, G, GV, h_new, h_post_remap, CS%tracer_Reg, &
                              CS%debug, dt=CS%dt_offline)
+      if (allocated(CS%tv%SpV_avg)) CS%tv%valid_SpV_halo = -1   ! Record that SpV_avg is no longer valid.
 
       do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
         h_new(i,j,k) = h_post_remap(i,j,k)