+Add thickness_to_dz, calc_derived_thermo and avg_specific_vol

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_density_integrals.F90

@@ -4,12 +4,13 @@ module MOM_density_integrals
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_EOS,              only : EOS_type
-use MOM_EOS,              only : EOS_quadrature
+use MOM_EOS,              only : EOS_quadrature, EOS_domain
 use MOM_EOS,              only : analytic_int_density_dz
 use MOM_EOS,              only : analytic_int_specific_vol_dp
 use MOM_EOS,              only : calculate_density
 use MOM_EOS,              only : calculate_spec_vol
 use MOM_EOS,              only : calculate_specific_vol_derivs
+use MOM_EOS,              only : average_specific_vol
 use MOM_error_handler,    only : MOM_error, FATAL, WARNING, MOM_mesg
 use MOM_hor_index,        only : hor_index_type
 use MOM_string_functions, only : uppercase
@@ -28,6 +29,7 @@ module MOM_density_integrals
 public int_specific_vol_dp
 public int_spec_vol_dp_generic_pcm
 public int_spec_vol_dp_generic_plm
+public avg_specific_vol
 public find_depth_of_pressure_in_cell
 
 contains
@@ -1613,6 +1615,36 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
 
 end subroutine find_depth_of_pressure_in_cell
 
+!> Calculate the average in situ specific volume across layers
+subroutine avg_specific_vol(T, S, p_t, dp, HI, EOS, SpV_avg, halo_size)
+  type(hor_index_type), intent(in)  :: HI  !< The horizontal index structure
+  real, dimension(SZI_(HI),SZJ_(HI)), &
+                        intent(in)  :: T   !< Potential temperature of the layer [C ~> degC]
+  real, dimension(SZI_(HI),SZJ_(HI)), &
+                        intent(in)  :: S   !< Salinity of the layer [S ~> ppt]
+  real, dimension(SZI_(HI),SZJ_(HI)), &
+                        intent(in)  :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
+  real, dimension(SZI_(HI),SZJ_(HI)), &
+                        intent(in)  :: dp  !< Pressure change in the layer [R L2 T-2 ~> Pa]
+  type(EOS_type),       intent(in)  :: EOS !< Equation of state structure
+  real, dimension(SZI_(HI),SZJ_(HI)), &
+                        intent(inout) :: SpV_avg !< The vertical average specific volume
+                                           !! in the layer [R-1 ~> m3 kg-1]
+  integer,    optional, intent(in)  :: halo_size !< The number of halo points in which to work.
+
+  ! Local variables
+  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
+  integer :: jsh, jeh, j, halo
+
+  halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0)
+  jsh = HI%jsc-halo ; jeh = HI%jec+halo
+
+  EOSdom(:) = EOS_domain(HI, halo_size)
+  do j=jsh,jeh
+    call average_specific_vol(T(:,j), S(:,j), p_t(:,j), dp(:,j), SpV_avg(:,j), EOS, EOSdom)
+  enddo
+
+end subroutine avg_specific_vol
 
 !> Returns change in anomalous pressure change from top to non-dimensional
 !! position pos between z_t and z_b [R L2 T-2 ~> Pa]

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