+Rescaled density units in MOM_T_S_init_from_Z

  Rescaled density units in MOM_temp_salt_initialize_from_Z for dimensional
consistency testing, and added a new optional argument, eps_rho, to
find_interfaces.  All answers are bitwise identical, but there is a new optional
argument to a public interface.

src/initialization/MOM_state_initialization.F90

@@ -1969,6 +1969,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
   integer :: kd, inconsistent
   integer :: nkd      ! number of levels to use for regridding input arrays
   real    :: eps_Z    ! A negligibly thin layer thickness [Z ~> m].
+  real    :: eps_rho  ! A negligibly small density difference [R ~> kg m-3].
   real    :: PI_180             ! for conversion from degrees to radians
 
   real, dimension(:,:), pointer :: shelf_area => NULL()
@@ -1988,9 +1989,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
   logical :: debug = .false.  ! manually set this to true for verbose output
 
   ! data arrays
-  real, dimension(:), allocatable :: z_edges_in, z_in, Rb
+  real, dimension(:), allocatable :: z_edges_in, z_in
+  real, dimension(:), allocatable :: Rb  ! Interface densities [R ~> kg m-3]
   real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z
-  real, dimension(:,:,:), allocatable :: rho_z
+  real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi   ! Interface heights [Z ~> m].
   real, dimension(SZI_(G),SZJ_(G))  :: nlevs
   real, dimension(SZI_(G))   :: press  ! Pressures [Pa].
@@ -2115,6 +2117,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
 
   !### Change this to GV%Angstrom_Z
   eps_z = 1.0e-10*US%m_to_Z
+  eps_rho = 1.0e-10*US%kg_m3_to_R
 
   ! Read input grid coordinates for temperature and salinity field
   ! in z-coordinate dataset. The file is REQUIRED to contain the
@@ -2154,7 +2157,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
   call convert_temp_salt_for_TEOS10(temp_z, salt_z, press, G, kd, mask_z, eos)
 
   do k=1,kd ; do j=js,je
-    call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos)
+    call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), is, ie, &
+                           eos, scale=US%kg_m3_to_R)
   enddo ; enddo
 
   call pass_var(temp_z,G%Domain)
@@ -2286,11 +2290,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
 
     ! Rb contains the layer interface densities
     allocate(Rb(nz+1))
-    do k=2,nz ; Rb(k) = 0.5*US%R_to_kg_m3*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo
-    Rb(1) = 0.0 ;  Rb(nz+1) = US%R_to_kg_m3*( 2.0*GV%Rlay(nz) - GV%Rlay(nz-1) )
+    do k=2,nz ; Rb(k) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo
+    Rb(1) = 0.0 ;  Rb(nz+1) = 2.0*GV%Rlay(nz) - GV%Rlay(nz-1)
 
     zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, Rb, G%bathyT(is:ie,js:je), &
-                                        nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z)
+                                        nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z, &
+                                        eps_rho=eps_rho)
 
     if (correct_thickness) then
       call adjustEtaToFitBathymetry(G, GV, US, zi, h)

src/initialization/midas_vertmap.F90

@@ -559,12 +559,12 @@ function find_limited_slope(val, e, k) result(slope)
 end function find_limited_slope
 
 !> Find interface positions corresponding to density profile
-function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps_z) result(zi)
+function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps_z, eps_rho) result(zi)
   real, dimension(:,:,:), &
-                      intent(in) :: rho   !< potential density in z-space [R ~> kg m-3]
+                      intent(in) :: rho   !< potential density in z-space [kg m-3 or R ~> kg m-3]
   real, dimension(size(rho,3)), &
-                      intent(in) :: zin   !< Input data levels [Z ~> m or m].
-  real, dimension(:), intent(in) :: Rb    !< target interface densities [R ~> kg m-3]
+                      intent(in) :: zin   !< Input data levels [m or Z ~> m].
+  real, dimension(:), intent(in) :: Rb    !< target interface densities [kg m-3 or R ~> kg m-3]
   real, dimension(size(rho,1),size(rho,2)), &
                       intent(in) :: depth !< ocean depth [Z ~> m].
   real, dimension(size(rho,1),size(rho,2)), &
@@ -573,7 +573,8 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps
   integer,  optional, intent(in) :: nkml  !< number of mixed layer pieces
   integer,  optional, intent(in) :: nkbl  !< number of buffer layer pieces
   real,     optional, intent(in) :: hml   !< mixed layer depth [Z ~> m].
-  real,     optional, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m or m].
+  real,     optional, intent(in) :: eps_z !< A negligibly small layer thickness [m or Z ~> m].
+  real,     optional, intent(in) :: eps_rho !< A negligibly small density difference [kg m-3 or R ~> kg m-3].
   real, dimension(size(rho,1),size(rho,2),size(Rb,1)) :: zi !< The returned interface, in the same units az zin.
 
   ! Local variables
@@ -589,8 +590,8 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps
   integer :: n,i,j,k,l,nx,ny,nz,nt
   integer :: nlay,kk,nkml_,nkbl_
   logical :: debug_ = .false.
-  real    :: epsln_Z    ! A negligibly thin layer thickness [Z ~> m].
-  real    :: epsln_rho  ! A negligibly small density change [kg m-3].
+  real    :: epsln_Z    ! A negligibly thin layer thickness [m or Z ~> m].
+  real    :: epsln_rho  ! A negligibly small density change [kg m-3 or R ~> kg m-3].
   real, parameter :: zoff=0.999
 
   nlay=size(Rb)-1
@@ -606,7 +607,7 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps
   nkbl_ = 0 ;  if (PRESENT(nkbl)) nkbl_ = max(0, nkbl)
   hml_ = 0.0 ; if (PRESENT(hml)) hml_ = hml
   epsln_Z = 1.0e-10 ; if (PRESENT(eps_z)) epsln_Z = eps_z
-  epsln_rho = 1.0e-10
+  epsln_rho = 1.0e-10 ; if (PRESENT(eps_rho)) epsln_rho = eps_rho
 
   if (PRESENT(nlevs)) then
     nlevs_data(:,:) = nlevs(:,:)