+Add thickness_to_dz and calc_derived_thermo

  Added the new overloaded interface thickness_to_dz to convert the layer
thicknesses in thickness units [H ~> m or kg m-2] into vertical distances in
[Z ~> m], with variants that set full 3-d arrays or an i-/k- slice.  Also added
a field (SpV_avg) for the layer-averaged specific volume to the thermo_vars_ptr
type and the new subroutine calc_derived_thermo to set it.  This new subroutine
is being called after halo updates to the temperatures and salinities.  The new
runtime parameter SEMI_BOUSSINESQ was added to determine whether tv%SpV_avg is
allocated and used; it is stored in GV%semi_Boussinesq.  Also added the new
element GV%dZ_subroundoff to the verticalGrid_type as a counterpart to
GV%H_subroundoff but in height units.  All answers are bitwise identical, but
there is a new runtime parameter in some MOM_parameter_doc files, new elements
in a transparent type and a new public interface.

src/core/MOM.F90

@@ -91,7 +91,7 @@ module MOM
 use MOM_grid,                  only : set_first_direction, rescale_grid_bathymetry
 use MOM_hor_index,             only : hor_index_type, hor_index_init
 use MOM_hor_index,             only : rotate_hor_index
-use MOM_interface_heights,     only : find_eta
+use MOM_interface_heights,     only : find_eta, calc_derived_thermo
 use MOM_interface_filter,      only : interface_filter, interface_filter_init, interface_filter_end
 use MOM_interface_filter,      only : interface_filter_CS
 use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end
@@ -1400,6 +1400,12 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
       call create_group_pass(pass_T_S, CS%tv%T, G%Domain, To_All+Omit_Corners, halo=1)
       call create_group_pass(pass_T_S, CS%tv%S, G%Domain, To_All+Omit_Corners, halo=1)
       call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass)
+      halo_sz = 1
+    endif
+
+    ! Update derived thermodynamic quantities.
+    if (allocated(CS%tv%SpV_avg)) then
+      call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz)
     endif
   endif
 
@@ -1581,6 +1587,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
     call create_group_pass(pass_uv_T_S_h, h, G%Domain, halo=dynamics_stencil)
     call do_group_pass(pass_uv_T_S_h, G%Domain, clock=id_clock_pass)
 
+    ! Update derived thermodynamic quantities.
+    if (allocated(tv%SpV_avg)) then
+      call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil)
+    endif
+
     if (CS%debug .and. CS%use_ALE_algorithm) then
       call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US)
       call hchksum(tv%T, "Post-ALE T", G%HI, haloshift=1, scale=US%C_to_degC)
@@ -1623,13 +1634,19 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
     call cpu_clock_end(id_clock_adiabatic)
 
     if (associated(tv%T)) then
-      call create_group_pass(pass_T_S, tv%T, G%Domain, To_All+Omit_Corners, halo=1)
-      call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=1)
+      dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo)
+      call create_group_pass(pass_T_S, tv%T, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil)
+      call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil)
       call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass)
       if (CS%debug) then
         if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, scale=US%C_to_degC)
         if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, scale=US%S_to_ppt)
       endif
+
+      ! Update derived thermodynamic quantities.
+      if (allocated(tv%SpV_avg)) then
+        call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil)
+      endif
     endif
 
   endif   ! endif for the block "if (.not.CS%adiabatic)"
@@ -1676,6 +1693,8 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
 
   type(time_type), pointer :: accumulated_time => NULL()
   type(time_type), pointer :: vertical_time => NULL()
+  integer :: dynamics_stencil  ! The computational stencil for the calculations
+                               ! in the dynamic core.
   integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
 
   ! 3D pointers
@@ -1848,6 +1867,12 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
 
   fluxes%fluxes_used = .true.
 
+  ! Update derived thermodynamic quantities.
+  if (allocated(CS%tv%SpV_avg)) then
+    dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo)
+    call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil)
+  endif
+
   if (last_iter) then
     accumulated_time = real_to_time(0.0)
   endif
@@ -2817,6 +2842,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     endif
   endif
 
+  ! 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)
+  endif
+
   if (use_ice_shelf .and. CS%debug) then
     call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0)
     call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0,scale=US%RZ_to_kg_m2)
@@ -3103,6 +3133,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   call do_group_pass(pass_uv_T_S_h, G%Domain)
 
+  ! 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)
+  endif
+
   if (associated(CS%visc%Kv_shear)) &
     call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1)
 
@@ -3931,6 +3966,7 @@ subroutine MOM_end(CS)
   if (associated(CS%Hml)) deallocate(CS%Hml)
   if (associated(CS%tv%salt_deficit)) deallocate(CS%tv%salt_deficit)
   if (associated(CS%tv%frazil)) deallocate(CS%tv%frazil)
+  if (allocated(CS%tv%SpV_avg)) deallocate(CS%tv%SpV_avg)
 
   if (associated(CS%tv%T)) then
     DEALLOC_(CS%T) ; CS%tv%T => NULL() ; DEALLOC_(CS%S) ; CS%tv%S => NULL()

src/core/MOM_interface_heights.F90

@@ -3,7 +3,7 @@ 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
+use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol
 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
@@ -16,18 +16,26 @@ module MOM_interface_heights
 
 #include <MOM_memory.h>
 
-public find_eta, dz_to_thickness, dz_to_thickness_simple
+public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
+public calc_derived_thermo
 
 !> Calculates the heights of the free surface or all interfaces from layer thicknesses.
 interface find_eta
   module procedure find_eta_2d, find_eta_3d
 end interface find_eta
 
-!> Calculates layer thickness in thickness units from geometric thicknesses in height units.
+!> Calculates layer thickness in thickness units from geometric distance between the
+!! interfaces around that layer in height units.
 interface dz_to_thickness
   module procedure dz_to_thickness_tv, dz_to_thickness_EoS
 end interface dz_to_thickness
 
+!> Converts layer thickness in thickness units into the vertical distance between the
+!! interfaces around a layer in height units.
+interface thickness_to_dz
+  module procedure thickness_to_dz_3d, thickness_to_dz_jslice
+end interface thickness_to_dz
+
 contains
 
 !> Calculates the heights of all interfaces between layers, using the appropriate
@@ -253,6 +261,45 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
 end subroutine find_eta_2d
 
 
+!> Calculate derived thermodynamic quantities for re-use later.
+subroutine calc_derived_thermo(tv, h, G, GV, US, halo)
+  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
+  type(thermo_var_ptrs),   intent(inout) :: tv !< A structure pointing to various
+                                               !! thermodynamic variables, some of
+                                               !! 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
+                                               !! calculate thicknesses
+  ! 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]
+  integer :: i, j, k, is, ie, js, je, halos, nz
+
+  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
+
+  if (allocated(tv%Spv_avg) .and. associated(tv%eqn_of_state)) then
+    if (associated(tv%p_surf)) then
+      do j=js,je ; do i=is,ie ; p_t(i,j) = tv%p_surf(i,j) ; enddo ; enddo
+    else
+      do j=js,je ; do i=is,ie ; p_t(i,j) = 0.0 ; enddo ; enddo
+    endif
+    do k=1,nz
+      do j=js,je ; do i=is,ie
+        dp(i,j) = GV%g_Earth*GV%H_to_RZ*h(i,j,k)
+      enddo ; enddo
+      call avg_specific_vol(tv%T(:,:,k), tv%S(:,:,k), p_t, dp, G%HI, tv%eqn_of_state, tv%SpV_avg(:,:,k), halo)
+      if (k<nz) then ; do j=js,je ; do i=is,ie
+        p_t(i,j) = p_t(i,j) + dp(i,j)
+      enddo ; enddo ; endif
+    enddo
+  endif
+
+end subroutine calc_derived_thermo
+
 !> Converts thickness from geometric height units to thickness units, perhaps via an
 !! inversion of the integral of the density in pressure using variables stored in
 !! the thermo_var_ptrs type when in non-Boussinesq mode.
@@ -428,4 +475,75 @@ subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode)
 
 end subroutine dz_to_thickness_simple
 
+!> Converts layer thicknesses in thickness units to the vertical distance between edges in height
+!! units, perhaps by multiplication by the precomputed layer-mean specific volume stored in an
+!! array in the thermo_var_ptrs type when in non-Boussinesq mode.
+subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size)
+  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
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(in)    :: h  !< Input thicknesses in thickness units [H ~> m or kg m-2].
+  type(thermo_var_ptrs),   intent(in)    :: tv !< A structure pointing to various
+                                               !! thermodynamic variables
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
+                                               !! This is essentially intent out, but declared as intent
+                                               !! inout to preserve any initialized values in halo points.
+  integer,       optional, intent(in)    :: halo_size !< Width of halo within which to
+                                               !! calculate thicknesses
+  ! Local variables
+  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
+    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
+  else
+    do k=1,nz ; do j=js,je ; do i=is,ie
+      dz(i,j,k) = GV%H_to_Z * h(i,j,k)
+    enddo ; enddo ; enddo
+  endif
+
+end subroutine thickness_to_dz_3d
+
+
+!> Converts a vertical i- / k- slice of layer thicknesses in thickness units to the vertical
+!! distance between edges in height units, perhaps by multiplication by the precomputed layer-mean
+!! specific volume stored in an array in the thermo_var_ptrs type when in non-Boussinesq mode.
+subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size)
+  type(ocean_grid_type),   intent(in)    :: G  !< The ocean's grid structure
+  type(verticalGrid_type), intent(in)    :: GV !< The ocean's vertical grid structure
+   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(in)    :: h  !< Input thicknesses in thickness units [H ~> m or kg m-2].
+  type(thermo_var_ptrs),   intent(in)    :: tv !< A structure pointing to various
+                                               !! thermodynamic variables
+  real, dimension(SZI_(G),SZK_(GV)), &
+                           intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m]
+                                               !! This is essentially intent out, but declared as intent
+                                               !! inout to preserve any initialized values in halo points.
+  integer,                 intent(in)    :: j  !< The second (j-) index of the input thicknesses to work with
+  integer,       optional, intent(in)    :: halo_size !< Width of halo within which to
+                                               !! calculate thicknesses
+  ! Local variables
+  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
+    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
+  else
+    do k=1,nz ; do i=is,ie
+      dz(i,k) = GV%H_to_Z * h(i,j,k)
+    enddo ; enddo
+  endif
+
+end subroutine thickness_to_dz_jslice
+
 end module MOM_interface_heights

src/core/MOM_variables.F90

@@ -93,6 +93,9 @@ module MOM_variables
   logical :: S_is_absS = .false. !< If true, the salinity variable tv%S is
                          !! actually the absolute salinity in units of [gSalt kg-1].
   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].
+
   ! These arrays are accumulated fluxes for communication with other components.
   real, dimension(:,:), pointer :: frazil => NULL()
                          !< The energy needed to heat the ocean column to the

src/core/MOM_verticalGrid.F90

@@ -41,12 +41,18 @@ module MOM_verticalGrid
 
   ! The following variables give information about the vertical grid.
   logical :: Boussinesq !< If true, make the Boussinesq approximation.
+  logical :: semi_Boussinesq !< If true, do non-Boussinesq pressure force calculations and
+                        !! use mass-based "thicknesses, but use Rho0 to convert layer thicknesses
+                        !! into certain height changes.  This only applies if BOUSSINESQ is false.
   real :: Angstrom_H    !< A one-Angstrom thickness in the model thickness units [H ~> m or kg m-2].
   real :: Angstrom_Z    !< A one-Angstrom thickness in the model depth units [Z ~> m].
   real :: Angstrom_m    !< A one-Angstrom thickness [m].
   real :: H_subroundoff !< A thickness that is so small that it can be added to a thickness of
                         !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2].
                         !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m.
+  real :: dZ_subroundoff !< A thickness in height units that is so small that it can be added to a
+                        !! vertical distance of Angstrom_Z or 1e-17 m without changing it at the bit
+                        !! level [Z ~> m].  This is the height equivalent of H_subroundoff.
   real, allocatable, dimension(:) :: &
     g_prime, &          !< The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2].
     Rlay                !< The target coordinate value (potential density) in each layer [R ~> kg m-3].
@@ -114,6 +120,12 @@ subroutine verticalGridInit( param_file, GV, US )
                  units="kg m-3", default=1035.0, scale=US%kg_m3_to_R)
   call get_param(param_file, mdl, "BOUSSINESQ", GV%Boussinesq, &
                  "If true, make the Boussinesq approximation.", default=.true.)
+  call get_param(param_file, mdl, "SEMI_BOUSSINESQ", GV%semi_Boussinesq, &
+                 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
+                 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
+                 "height changes.  This only applies if BOUSSINESQ is false.", &
+                 default=.true., do_not_log=GV%Boussinesq)
+  if (GV%Boussinesq) GV%semi_Boussinesq = .true.
   call get_param(param_file, mdl, "ANGSTROM", GV%Angstrom_Z, &
                  "The minimum layer thickness, usually one-Angstrom.", &
                  units="m", default=1.0e-10, scale=US%m_to_Z)
@@ -165,7 +177,8 @@ subroutine verticalGridInit( param_file, GV, US )
     GV%Angstrom_H = US%Z_to_m*GV%Angstrom_Z * 1000.0*GV%kg_m2_to_H
     GV%H_to_MKS = GV%H_to_kg_m2
   endif
-  GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H,GV%m_to_H*1e-17)
+  GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H, GV%m_to_H*1e-17)
+  GV%dZ_subroundoff = 1e-20 * max(GV%Angstrom_Z, US%m_to_Z*1e-17)
   GV%H_to_Pa = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth * GV%H_to_kg_m2
 
   GV%H_to_Z = GV%H_to_m * US%m_to_Z

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -43,7 +43,7 @@ module MOM_diabatic_driver
 use MOM_grid,                only : ocean_grid_type
 use MOM_int_tide_input,      only : set_int_tide_input, int_tide_input_init
 use MOM_int_tide_input,      only : int_tide_input_end, int_tide_input_CS, int_tide_input_type
-use MOM_interface_heights,   only : find_eta
+use MOM_interface_heights,   only : find_eta, calc_derived_thermo
 use MOM_internal_tides,      only : propagate_int_tide
 use MOM_internal_tides,      only : internal_tides_init, internal_tides_end, int_tide_CS
 use MOM_kappa_shear,         only : kappa_shear_is_used
@@ -1828,9 +1828,15 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
   ! Also changes: visc%Kd_shear and visc%Kv_shear
   if ((CS%halo_TS_diff > 0) .and. (CS%ML_mix_first > 0.0)) then
     if (associated(tv%T)) call pass_var(tv%T, G%Domain, halo=CS%halo_TS_diff, complete=.false.)
-    if (associated(tv%T)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.)
+    if (associated(tv%S)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.)
     call pass_var(h, G%domain, halo=CS%halo_TS_diff, complete=.true.)
   endif
+
+  ! Update derived thermodynamic quantities.
+  if ((CS%ML_mix_first > 0.0) .and. allocated(tv%SpV_avg)) then
+    call calc_derived_thermo(tv, h, G, GV, US, halo=CS%halo_TS_diff)
+  endif
+
   if (CS%debug) &
     call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff)
   if (CS%double_diffuse) then

src/tracer/MOM_offline_main.F90

@@ -22,6 +22,7 @@ module MOM_offline_main
 use MOM_file_parser,          only : read_param, get_param, log_version, param_file_type
 use MOM_forcing_type,         only : forcing
 use MOM_grid,                 only : ocean_grid_type
+use MOM_interface_heights,    only : calc_derived_thermo
 use MOM_io,                   only : MOM_read_data, MOM_read_vector, CENTER
 use MOM_offline_aux,          only : update_offline_from_arrays, update_offline_from_files
 use MOM_offline_aux,          only : next_modulo_time, offline_add_diurnal_sw
@@ -1025,6 +1026,7 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
   type(forcing),              intent(inout) :: fluxes !< Pointers to forcing fields
   logical,                    intent(in   ) :: do_ale !< True if using ALE
   ! Local variables
+  integer :: stencil
   integer :: i, j, k, is, ie, js, je, nz
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_start ! Initial thicknesses [H ~> m or kg m-2]
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
@@ -1086,6 +1088,12 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale)
   call pass_var(CS%tv%T, G%Domain)
   call pass_var(CS%tv%S, G%Domain)
 
+  ! Update derived thermodynamic quantities.
+  if (allocated(CS%tv%SpV_avg)) then
+    stencil = min(3, G%Domain%nihalo, G%Domain%njhalo)
+    call calc_derived_thermo(CS%tv, CS%h_end, G, GV, US, halo=stencil)
+  endif
+
   ! Update the read indices
   CS%ridx_snap = next_modulo_time(CS%ridx_snap,CS%numtime)
   CS%ridx_sum = next_modulo_time(CS%ridx_sum,CS%numtime)