diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 62b06a0b57..1907a75c74 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -149,6 +149,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                      optional, intent(in)   :: frac_shelf_h    !< The fraction of the grid cell covered
                                                                !! by a floating ice shelf [nondim].
   ! Local variables
+  real :: depth_tot(SZI_(G),SZJ_(G))  ! The nominal total depth of the ocean [Z ~> m]
   character(len=200) :: filename   ! The name of an input file.
   character(len=200) :: filename2  ! The name of an input files.
   character(len=200) :: inputdir   ! The directory where NetCDF input files are.
@@ -179,8 +180,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   logical :: debug_layers = .false.
   logical :: use_ice_shelf
   character(len=80) :: mesg
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
 
@@ -227,6 +228,13 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
     !enddo
   endif
 
+  ! Set the nominal depth of the ocean, which might be different from the bathymetric
+  ! geopotential height, for use by the various initialization routines.  G%bathyT has
+  ! already been initialized in previous calls.
+  do j=jsd,jed ; do i=isd,ied
+    depth_tot(i,j) = G%bathyT(i,j)
+  enddo ; enddo
+
   ! The remaining initialization calls are done, regardless of whether the
   ! fields are actually initialized here (if just_read=.false.) or whether it
   ! is just to make sure that all valid parameters are read to enable the
@@ -241,8 +249,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
     if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//&
        "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")
 
-    call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params=just_read,&
-         frac_shelf_h=frac_shelf_h)
+    call MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, &
+                                         just_read_params=just_read, frac_shelf_h=frac_shelf_h)
   else
     ! Initialize thickness, h.
     call get_param(PF, mdl, "THICKNESS_CONFIG", config, &
@@ -275,9 +283,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
              default="uniform", do_not_log=just_read)
     select case (trim(config))
       case ("file")
-        call initialize_thickness_from_file(h, G, GV, US, PF, .false., just_read_params=just_read)
+        call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, .false., just_read_params=just_read)
       case ("thickness_file")
-        call initialize_thickness_from_file(h, G, GV, US, PF, .true., just_read_params=just_read)
+        call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, .true., just_read_params=just_read)
       case ("coord")
         if (new_sim .and. useALE) then
           call ALE_initThicknessToCoord( ALE_CSp, G, GV, h )
@@ -285,37 +293,37 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
           call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//&
                                 "for THICKNESS_CONFIG of 'coord'")
         endif
-      case ("uniform"); call initialize_thickness_uniform(h, G, GV, PF, &
+      case ("uniform"); call initialize_thickness_uniform(h, depth_tot, G, GV, PF, &
                                  just_read_params=just_read)
-      case ("list"); call initialize_thickness_list(h, G, GV, US, PF, &
+      case ("list"); call initialize_thickness_list(h, depth_tot, G, GV, US, PF, &
                                  just_read_params=just_read)
-      case ("DOME"); call DOME_initialize_thickness(h, G, GV, PF, &
+      case ("DOME"); call DOME_initialize_thickness(h, depth_tot, G, GV, PF, &
                               just_read_params=just_read)
-      case ("ISOMIP"); call ISOMIP_initialize_thickness(h, G, GV, US, PF, tv, &
+      case ("ISOMIP"); call ISOMIP_initialize_thickness(h, depth_tot, G, GV, US, PF, tv, &
                                 just_read_params=just_read)
-      case ("benchmark"); call benchmark_initialize_thickness(h, G, GV, US, PF, &
+      case ("benchmark"); call benchmark_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                    tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
-      case ("Neverwoorld","Neverland"); call Neverworld_initialize_thickness(h, G, GV, US, PF, &
+      case ("Neverworld","Neverland"); call Neverworld_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                 tv%eqn_of_state, tv%P_Ref)
-      case ("search"); call initialize_thickness_search
-      case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, GV, PF, &
+      case ("search"); call initialize_thickness_search()
+      case ("circle_obcs"); call circle_obcs_initialize_thickness(h, depth_tot, G, GV, PF, &
                                      just_read_params=just_read)
       case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, GV, US, &
                                        PF, just_read_params=just_read)
       case ("external_gwave"); call external_gwave_initialize_thickness(h, G, GV, US, &
                                         PF, just_read_params=just_read)
-      case ("DOME2D"); call DOME2d_initialize_thickness(h, G, GV, US, PF, &
+      case ("DOME2D"); call DOME2d_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                 just_read_params=just_read)
       case ("adjustment2d"); call adjustment_initialize_thickness(h, G, GV, US, &
                                       PF, just_read_params=just_read)
-      case ("sloshing"); call sloshing_initialize_thickness(h, G, GV, US, PF, &
+      case ("sloshing"); call sloshing_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                   just_read_params=just_read)
-      case ("seamount"); call seamount_initialize_thickness(h, G, GV, US, PF, &
+      case ("seamount"); call seamount_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                   just_read_params=just_read)
-      case ("dumbbell"); call dumbbell_initialize_thickness(h, G, GV, US, PF, &
+      case ("dumbbell"); call dumbbell_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                   just_read_params=just_read)
-      case ("soliton"); call soliton_initialize_thickness(h, G, GV, US)
-      case ("phillips"); call Phillips_initialize_thickness(h, G, GV, US, PF, &
+      case ("soliton"); call soliton_initialize_thickness(h, depth_tot, G, GV, US)
+      case ("phillips"); call Phillips_initialize_thickness(h, depth_tot, G, GV, US, PF, &
                                   just_read_params=just_read)
       case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, US, &
                                       PF, just_read_params=just_read)
@@ -363,11 +371,11 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
         case ("DOME2D"); call DOME2d_initialize_temperature_salinity ( tv%T, &
                                   tv%S, h, G, GV, PF, eos, just_read_params=just_read)
         case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity ( tv%T, &
-                                  tv%S, h, G, GV, US, PF, eos, just_read_params=just_read)
+                                  tv%S, h, depth_tot, G, GV, US, PF, eos, just_read_params=just_read)
         case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
-                                        tv%S, h, G, GV, PF, eos, just_read_params=just_read)
+                                        tv%S, h, depth_tot, G, GV, PF, eos, just_read_params=just_read)
         case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, &
-                                           tv%S, h, G, GV, US, PF, just_read_params=just_read)
+                                           tv%S, h, depth_tot, G, GV, US, PF, just_read_params=just_read)
         case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
                                     tv%S, h, G, GV, PF, eos, just_read_params=just_read)
         case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
@@ -547,22 +555,22 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                  " \t\t for buoyancy-forced basin case.\n"//&
                  " \t USER - call a user modified routine.", default="file")
     select case (trim(config))
-      case ("DOME"); call DOME_initialize_sponges(G, GV, US, tv, PF, sponge_CSp)
-      case ("DOME2D"); call DOME2d_initialize_sponges(G, GV, US, tv, PF, useALE, &
+      case ("DOME"); call DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, sponge_CSp)
+      case ("DOME2D"); call DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
                                                       sponge_CSp, ALE_sponge_CSp)
-      case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, US, tv, PF, useALE, &
+      case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
                                                       sponge_CSp, ALE_sponge_CSp)
-      case("RGC"); call RGC_initialize_sponges(G, GV, US, tv, u, v, PF, useALE, &
+      case("RGC"); call RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, useALE, &
                                                      sponge_CSp, ALE_sponge_CSp)
       case ("USER"); call user_initialize_sponges(G, GV, use_temperature, tv, PF, sponge_CSp, h)
-      case ("BFB"); call BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, PF, &
+      case ("BFB"); call BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, PF, &
                                                           sponge_CSp, h)
-      case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, PF, useALE, &
+      case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
                                                           sponge_CSp, ALE_sponge_CSp)
       case ("phillips"); call Phillips_initialize_sponges(G, GV, US, tv, PF, sponge_CSp, h)
-      case ("dense"); call dense_water_initialize_sponges(G, GV, US, tv, PF, useALE, &
+      case ("dense"); call dense_water_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
                                                           sponge_CSp, ALE_sponge_CSp)
-      case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, PF, &
+      case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, PF, &
                                                   sponge_CSp, ALE_sponge_CSp, Time)
       case default ; call MOM_error(FATAL,  "MOM_initialize_state: "//&
              "Unrecognized sponge configuration "//trim(config))
@@ -639,14 +647,16 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
 end subroutine MOM_initialize_state
 
 !> Reads the layer thicknesses or interface heights from a file.
-subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thickness, &
+subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, file_has_thickness, &
                                           just_read_params)
   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(out) :: h    !< The thickness that is being initialized [H ~> m or kg m-2].
-  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
+  type(param_file_type),   intent(in)  :: param_file !< A structure indicating the open file
                                                !! to parse for model parameter values.
   logical,                 intent(in)  :: file_has_thickness !< If true, this file contains layer
                                                !! thicknesses; otherwise it contains
@@ -696,7 +706,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi
     call MOM_read_data(filename, "eta", eta(:,:,:), G%Domain, scale=US%m_to_Z)
 
     if (correct_thickness) then
-      call adjustEtaToFitBathymetry(G, GV, US, eta, h,  dZ_ref_eta=0.0)
+      call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=0.0)
     else
       do k=nz,1,-1 ; do j=js,je ; do i=is,ie
         if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then
@@ -708,7 +718,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi
       enddo ; enddo ; enddo
 
       do j=js,je ; do i=is,ie
-        if (abs(eta(i,j,nz+1) + G%bathyT(i,j)) > 1.0*US%m_to_Z) &
+        if (abs(eta(i,j,nz+1) + depth_tot(i,j)) > 1.0*US%m_to_Z) &
           inconsistent = inconsistent + 1
       enddo ; enddo
       call sum_across_PEs(inconsistent)
@@ -811,11 +821,13 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta)
 end subroutine adjustEtaToFitBathymetry
 
 !> Initializes thickness to be uniform
-subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params)
+subroutine initialize_thickness_uniform(h, depth_tot, G, GV, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot   !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -850,7 +862,7 @@ subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params)
     ! Angstrom thick, and 2.  the interfaces are where they should be
     ! based on the resting depths and interface height perturbations,
     ! as long at this doesn't interfere with 1.
-    eta1D(nz+1) = -G%bathyT(i,j)
+    eta1D(nz+1) = -depth_tot(i,j)
     do k=nz,1,-1
       eta1D(K) = e0(K)
       if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
@@ -866,12 +878,14 @@ subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params)
 end subroutine initialize_thickness_uniform
 
 !> Initialize thickness from a 1D list
-subroutine initialize_thickness_list(h, G, GV, US, param_file, just_read_params)
+subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -926,7 +940,7 @@ subroutine initialize_thickness_list(h, G, GV, US, param_file, just_read_params)
     ! Angstrom thick, and 2.  the interfaces are where they should be
     ! based on the resting depths and interface height perturbations,
     ! as long at this doesn't interfere with 1.
-    eta1D(nz+1) = -G%bathyT(i,j)
+    eta1D(nz+1) = -depth_tot(i,j)
     do k=nz,1,-1
       eta1D(K) = e0(K)
       if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
@@ -1743,10 +1757,11 @@ end subroutine initialize_temp_salt_linear
 !! number of tracers should be restored within each sponge. The
 !! interface height is always subject to damping, and must always be
 !! the first registered field.
-subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_file, Layer_CSp, ALE_CSp, Time)
+subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, param_file, &
+                                   Layer_CSp, ALE_CSp, Time)
   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(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
   logical,                 intent(in) :: use_temperature !< If true, T & S are state variables.
   type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic
                                               !! variables.
@@ -1756,6 +1771,8 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f
   real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                            intent(in) :: v    !< The meridional velocity that is being
                                               !! initialized [L T-1 ~> m s-1]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
   type(sponge_CS),         pointer    :: Layer_CSp  !< A pointer that is set to point to the control
                                               !! structure for this module (in layered mode).
@@ -1916,7 +1933,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f
     call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z)
 
     do j=js,je ; do i=is,ie
-      eta(i,j,nz+1) = -G%bathyT(i,j)
+      eta(i,j,nz+1) = -depth_tot(i,j)
     enddo ; enddo
     do k=nz,1,-1 ; do j=js,je ; do i=is,ie
       if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) &
@@ -1977,7 +1994,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f
       allocate(h(isd:ied,jsd:jed,nz_data))
       call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z)
       do j=js,je ; do i=is,ie
-        eta(i,j,nz+1) = -G%bathyT(i,j)
+        eta(i,j,nz+1) = -depth_tot(i,j)
       enddo ; enddo
       do k=nz,1,-1 ; do j=js,je ; do i=is,ie
         if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) &
@@ -2256,13 +2273,15 @@ end subroutine set_velocity_depth_min
 !> This subroutine determines the isopycnal or other coordinate interfaces and
 !! layer potential temperatures and salinities directly from a z-space file on
 !! a latitude-longitude grid.
-subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params, frac_shelf_h)
+subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just_read_params, frac_shelf_h)
   type(ocean_grid_type),   intent(inout) :: 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(out)   :: h    !< Layer thicknesses being initialized [H ~> m or kg m-2]
   type(thermo_var_ptrs),   intent(inout) :: tv   !< A structure pointing to various thermodynamic
                                                  !! variables including temperature and salinity
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)    :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
   type(param_file_type),   intent(in)    :: PF   !< A structure indicating the open file
                                                  !! to parse for model parameter values.
@@ -2521,7 +2540,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
   call pass_var(rho_z,G%Domain)
 
   do j=js,je ; do i=is,ie
-    Z_bottom(i,j) = -G%bathyT(i,j)
+    Z_bottom(i,j) = -depth_tot(i,j)
   enddo ; enddo
 
   ! Done with horizontal interpolation.
diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90
index f632b95086..49c0a03235 100644
--- a/src/user/BFB_initialization.F90
+++ b/src/user/BFB_initialization.F90
@@ -76,13 +76,15 @@ end subroutine BFB_set_coord
 
 !> This subroutine sets up the sponges for the southern bouundary of the domain. Maximum damping occurs
 !! within 2 degrees lat of the boundary. The damping linearly decreases northward over the next 2 degrees.
-subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, param_file, CSp, h)
+subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, param_file, CSp, h)
   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
   logical,                 intent(in) :: use_temperature !< If true, temperature and salinity are used as
                                             !! state variables.
   type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic variables
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters
   type(sponge_CS),         pointer    :: CSp  !< A pointer to the sponge control structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
@@ -129,7 +131,7 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, para
   max_damping = 1.0  / (86400.0*US%s_to_T)
 
   do j=js,je ; do i=is,ie
-    if (G%bathyT(i,j) <= min_depth) then ; Idamp(i,j) = 0.0
+    if (depth_tot(i,j) <= min_depth) then ; Idamp(i,j) = 0.0
     elseif (G%geoLatT(i,j) < slat+2.0) then ; Idamp(i,j) = max_damping
     elseif (G%geoLatT(i,j) < slat+4.0) then
       Idamp(i,j) = max_damping * (slat+4.0-G%geoLatT(i,j))/2.0
diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90
index 293d601757..f99f0b8d5c 100644
--- a/src/user/DOME2d_initialization.F90
+++ b/src/user/DOME2d_initialization.F90
@@ -90,12 +90,14 @@ subroutine DOME2d_initialize_topography( D, G, param_file, max_depth )
 end subroutine DOME2d_initialize_topography
 
 !> Initialize thicknesses according to coordinate mode
-subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_params )
+subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params )
   type(ocean_grid_type),   intent(in)  :: G  !< Ocean grid structure
   type(verticalGrid_type), intent(in)  :: GV !< Vertical grid structure
   type(unit_scale_type),   intent(in)  :: US !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                            intent(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -150,7 +152,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
     case ( REGRIDDING_LAYER, REGRIDDING_RHO )
 
       do j=js,je ; do i=is,ie
-        eta1D(nz+1) = -G%bathyT(i,j)
+        eta1D(nz+1) = -depth_tot(i,j)
         do k=nz,1,-1
           eta1D(k) = e0(k)
           if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
@@ -172,7 +174,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
  !  case ( IC_RHO_C )
  !
  !    do j=js,je ; do i=is,ie
- !       eta1D(nz+1) = -G%bathyT(i,j)
+ !       eta1D(nz+1) = -depth_tot(i,j)
  !       do k=nz,1,-1
  !         eta1D(k) = e0(k)
  !         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
@@ -194,7 +196,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
     case ( REGRIDDING_ZSTAR )
 
       do j=js,je ; do i=is,ie
-        eta1D(nz+1) = -G%bathyT(i,j)
+        eta1D(nz+1) = -depth_tot(i,j)
         do k=nz,1,-1
           eta1D(k) = e0(k)
           if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
@@ -208,7 +210,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
 
     case ( REGRIDDING_SIGMA )
       do j=js,je ; do i=is,ie
-        h(i,j,:) = GV%Z_to_H*G%bathyT(i,j) / nz
+        h(i,j,:) = GV%Z_to_H*depth_tot(i,j) / nz
       enddo ; enddo
 
     case default
@@ -353,11 +355,13 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file,
 end subroutine DOME2d_initialize_temperature_salinity
 
 !> Set up sponges in 2d DOME configuration
-subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp)
+subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
   type(ocean_grid_type),   intent(in) :: G  !< Ocean grid structure
   type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
   type(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
   type(thermo_var_ptrs),   intent(in) :: tv !< Thermodynamics structure
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file !< Parameter file structure
   logical,                 intent(in) :: use_ALE !< If true, indicates model is in ALE mode
   type(sponge_CS),         pointer    :: CSp !< Layer-mode sponge structure
@@ -453,7 +457,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
     enddo
     e0(nz+1) = -G%max_depth
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) = e0(k)
         if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
@@ -470,7 +474,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
     ! Construct temperature and salinity on the arbitrary grid
     T(:,:,:) = 0.0 ; S(:,:,:) = 0.0
     do j=js,je ; do i=is,ie
-      z = -G%bathyT(i,j)
+      z = -depth_tot(i,j)
       do k = nz,1,-1
         z = z + 0.5 * GV%H_to_Z * h(i,j,k) ! Position of the center of layer k
         S(i,j,k) = 34.0 - 1.0 * (z / (G%max_depth))
@@ -491,7 +495,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
 
     ! Construct interface heights to restore toward
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(K) = -G%max_depth * real(k-1) / real(nz)
         if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
@@ -508,7 +512,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
         d_eta(nz) = dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_Z
       endif
 
-      eta(i,j,nz+1) = -G%bathyT(i,j)
+      eta(i,j,nz+1) = -depth_tot(i,j)
       do K=nz,1,-1
         eta(i,j,K) = eta(i,j,K+1) + d_eta(k)
       enddo
diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90
index 81444704b3..1f3d24e1c9 100644
--- a/src/user/DOME_initialization.F90
+++ b/src/user/DOME_initialization.F90
@@ -87,11 +87,13 @@ end subroutine DOME_initialize_topography
 
 ! -----------------------------------------------------------------------------
 !> This subroutine initializes layer thicknesses for the DOME experiment
-subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params)
+subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot   !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -124,7 +126,7 @@ subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params)
 !  Angstrom thick, and 2.  the interfaces are where they should be   !
 !  based on the resting depths and interface height perturbations,   !
 !  as long at this doesn't interfere with 1.                         !
-    eta1D(nz+1) = -G%bathyT(i,j)
+    eta1D(nz+1) = -depth_tot(i,j)
     do k=nz,1,-1
       eta1D(K) = e0(K)
       if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
@@ -145,17 +147,19 @@ end subroutine DOME_initialize_thickness
 !! number of tracers should be restored within each sponge. The       !
 !! interface height is always subject to damping, and must always be  !
 !! the first registered field.                                        !
-subroutine DOME_initialize_sponges(G, GV, US, tv, PF, CSp)
-  type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure.
+subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
+  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(in) :: tv   !< A structure containing pointers to any available
-                               !! thermodynamic fields, including potential temperature and
-                               !! salinity or mixed layer density. Absent fields have NULL ptrs.
-  type(param_file_type), intent(in) :: PF   !< A structure indicating the open file to
-                                            !! parse for model parameter values.
-  type(sponge_CS),       pointer    :: CSp  !< A pointer that is set to point to the control
-                                            !! structure for this module.
+  type(thermo_var_ptrs),   intent(in) :: tv   !< A structure containing pointers to any available
+                                !! thermodynamic fields, including potential temperature and
+                                !! salinity or mixed layer density. Absent fields have NULL ptrs.
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
+  type(param_file_type),   intent(in) :: PF   !< A structure indicating the open file to
+                                              !! parse for model parameter values.
+  type(sponge_CS),         pointer    :: CSp  !< A pointer that is set to point to the control
+                                              !! structure for this module.
 
   real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m].
   real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables. !
@@ -204,16 +208,16 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, PF, CSp)
     ! depth space for Boussinesq or non-Boussinesq models.
     eta(i,j,1) = 0.0
     do k=2,nz
-!     eta(i,j,K)=max(H0(k), -G%bathyT(i,j), GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j))
-      e_dense = -G%bathyT(i,j)
+!     eta(i,j,K) = max(H0(k), -depth_tot(i,j), GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j))
+      e_dense = -depth_tot(i,j)
       if (e_dense >= H0(k)) then ; eta(i,j,K) = e_dense
       else ; eta(i,j,K) = H0(k) ; endif
-      if (eta(i,j,K) < GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j)) &
-          eta(i,j,K) = GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j)
+      if (eta(i,j,K) < GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j)) &
+          eta(i,j,K) = GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j)
     enddo
-    eta(i,j,nz+1) = -G%bathyT(i,j)
+    eta(i,j,nz+1) = -depth_tot(i,j)
 
-    if (G%bathyT(i,j) > min_depth) then
+    if (depth_tot(i,j) > min_depth) then
       Idamp(i,j) = damp / 86400.0
     else ; Idamp(i,j) = 0.0 ; endif
   enddo ; enddo
diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90
index aa1c6cdfe6..76f60d9b99 100644
--- a/src/user/ISOMIP_initialization.F90
+++ b/src/user/ISOMIP_initialization.F90
@@ -128,12 +128,14 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US)
 end subroutine ISOMIP_initialize_topography
 
 !> Initialization of thicknesses
-subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read_params)
+subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot   !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   type(thermo_var_ptrs),   intent(in)  :: tv          !< A structure containing pointers to any
@@ -206,7 +208,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
 
     ! Calculate thicknesses
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) = e0(k)
         if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
@@ -221,7 +223,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
   case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR )   ! Initial thicknesses for z coordinates
     if (just_read) return ! All run-time parameters have been read, so return.
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) =  -G%max_depth * real(k-1) / real(nz)
         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
@@ -236,7 +238,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
   case ( REGRIDDING_SIGMA )             ! Initial thicknesses for sigma coordinates
     if (just_read) return ! All run-time parameters have been read, so return.
     do j=js,je ; do i=is,ie
-      h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz)
+      h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz)
     enddo ; enddo
 
   case default
@@ -248,7 +250,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
 end subroutine ISOMIP_initialize_thickness
 
 !> Initial values for temperature and salinity
-subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, &
+subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, US, param_file, &
                                                     eqn_of_state, just_read_params)
   type(ocean_grid_type),                     intent(in)  :: G  !< Ocean grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< Vertical grid structure
@@ -256,6 +258,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T  !< Potential temperature [degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< Layer thickness [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G)),          intent(in)  :: depth_tot  !< The nominal total bottom-to-top
+                                                               !! depth of the ocean [Z ~> m]
   type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
   type(EOS_type),                            pointer     :: eqn_of_state !< Equation of state structure
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -315,7 +319,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
       dS_dz = (s_sur - s_bot) / G%max_depth
       dT_dz = (t_sur - t_bot) / G%max_depth
       do j=js,je ; do i=is,ie
-        xi0 = -G%bathyT(i,j)
+        xi0 = -depth_tot(i,j)
         do k = nz,1,-1
           xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer
           S(i,j,k) = S_sur + dS_dz * xi0
@@ -420,7 +424,7 @@ end subroutine ISOMIP_initialize_temperature_salinity
 !> Sets up the the inverse restoration time (Idamp), and
 ! the values towards which the interface heights and an arbitrary
 ! number of tracers should be restored within each sponge.
-subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
+subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, ACSp)
   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
@@ -429,6 +433,8 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
                                             !! fields, potential temperature and
                                             !! salinity or mixed layer density.
                                             !! Absent fields have NULL ptrs.
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type), intent(in) :: PF   !< A structure indicating the
                                             !! open file to parse for model
                                             !! parameter values.
@@ -508,7 +514,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
   !  and mask2dT is 1.
 
   do j=js,je ; do i=is,ie
-    if (G%bathyT(i,j) <= min_depth) then
+    if (depth_tot(i,j) <= min_depth) then
       Idamp(i,j) = 0.0
     elseif (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then
       dummy1 = (G%geoLonT(i,j)-790.0)/(800.0-790.0)
@@ -549,7 +555,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
 
         ! Calculate thicknesses
         do j=js,je ; do i=is,ie
-          eta1D(nz+1) = -G%bathyT(i,j)
+          eta1D(nz+1) = -depth_tot(i,j)
           do k=nz,1,-1
             eta1D(k) = e0(k)
             if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
@@ -563,7 +569,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
 
       case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR )  ! Initial thicknesses for z coordinates
         do j=js,je ; do i=is,ie
-          eta1D(nz+1) = -G%bathyT(i,j)
+          eta1D(nz+1) = -depth_tot(i,j)
           do k=nz,1,-1
             eta1D(k) =  -G%max_depth * real(k-1) / real(nz)
             if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
@@ -577,7 +583,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
 
       case ( REGRIDDING_SIGMA )             ! Initial thicknesses for sigma coordinates
         do j=js,je ; do i=is,ie
-          h(i,j,:) = GV%Z_to_H * (G%bathyT(i,j) / dfloat(nz))
+          h(i,j,:) = GV%Z_to_H * (depth_tot(i,j) / dfloat(nz))
         enddo ; enddo
 
       case default
@@ -593,7 +599,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
     dS_dz = (s_sur - s_bot) / G%max_depth
     dT_dz = (t_sur - t_bot) / G%max_depth
     do j=js,je ; do i=is,ie
-      xi0 = -G%bathyT(i,j)
+      xi0 = -depth_tot(i,j)
       do k = nz,1,-1
         xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer
         S(i,j,k) = S_sur + dS_dz * xi0
diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90
index 4a136dd2db..ed944e5f0a 100644
--- a/src/user/Kelvin_initialization.F90
+++ b/src/user/Kelvin_initialization.F90
@@ -186,6 +186,7 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   real :: lambda       ! Offshore decay scale [L-1 ~> m-1]
   real :: omega        ! Wave frequency [T-1 ~> s-1]
   real :: PI
+  real :: depth_tot(SZI_(G),SZJ_(G))  ! The total depth of the ocean [Z ~> m]
   integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
   integer :: IsdB, IedB, JsdB, JedB
   real    :: mag_SSH ! An overall magnitude of the external wave sea surface height at the coastline [Z ~> m]
@@ -209,6 +210,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
   PI = 4.0*atan(1.0)
   km_to_L_scale = 1000.0*US%m_to_L
 
+  do j=jsd,jed ; do i=isd,ied
+    depth_tot(i,j) = G%bathyT(i,j)
+  enddo ; enddo
+  !### Instead this should be:
+  ! do j=jsd,jed ; do i=isd,ied
+  !   depth_tot(i,j) = 0.0
+  ! enddo ; enddo
+  ! do j=jsd,jed ; do i=isd,ied
+  !   depth_tot(i,j) = depth_tot(i,j) + GV%H_to_Z * h(i,j,k)
+  ! enddo ; enddo
+
   if (CS%mode == 0) then
     mag_SSH = 1.0*US%m_to_Z
     omega = 2.0 * PI / (12.42 * 3600.0*US%s_to_T)      ! M2 Tide period
@@ -245,20 +257,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
         y = -(x1 - CS%coast_offset1) * sina + y1 * cosa
         if (CS%mode == 0) then
           ! Use inside bathymetry
-          cff = sqrt(GV%g_Earth * G%bathyT(i+1,j) )
+          cff = sqrt(GV%g_Earth * depth_tot(i+1,j) )
           val2 = mag_SSH * exp(- CS%F_0 * y / cff)
           segment%eta(I,j) = GV%Z_to_H*val2 * cos(omega * time_sec)
-          segment%normal_vel_bt(I,j) = (val2 * (val1 * cff * cosa / &
-                 (G%bathyT(i+1,j) )) )
+          segment%normal_vel_bt(I,j) = val2 * (val1 * cff * cosa / depth_tot(i+1,j) )
           if (segment%nudged) then
             do k=1,nz
-              segment%nudged_normal_vel(I,j,k) = (val2 * (val1 * cff * cosa / &
-                     (G%bathyT(i+1,j))) )
+              segment%nudged_normal_vel(I,j,k) = val2 * (val1 * cff * cosa / depth_tot(i+1,j) )
             enddo
           elseif (segment%specified) then
             do k=1,nz
-              segment%normal_vel(I,j,k) = (val2 * (val1 * cff * cosa / &
-                     (G%bathyT(i+1,j) )) )
+              segment%normal_vel(I,j,k) = val2 * (val1 * cff * cosa / depth_tot(i+1,j) )
               segment%normal_trans(I,j,k) = segment%normal_vel(I,j,k) * h(i+1,j,k) * G%dyCu(I,j)
             enddo
           endif
@@ -288,11 +297,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
           y1 = km_to_L_scale * G%geoLatBu(I,J)
           x = (x1 - CS%coast_offset1) * cosa + y1 * sina
           y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
-          cff = sqrt(GV%g_Earth * G%bathyT(i+1,j) )
+          cff = sqrt(GV%g_Earth * depth_tot(i+1,j) )
           val2 = mag_SSH * exp(- CS%F_0 * y / cff)
           if (CS%mode == 0) then ; do k=1,nz
             segment%tangential_vel(I,J,k) = (val1 * val2 * cff * sina) / &
-               ( 0.5*(G%bathyT(i+1,j+1) +  G%bathyT(i+1,j) ) )
+               ( 0.5*(depth_tot(i+1,j+1) +  depth_tot(i+1,j) ) )
 
           enddo ; endif
         enddo ; enddo
@@ -306,20 +315,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
         x = (x1 - CS%coast_offset1) * cosa + y1 * sina
         y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
         if (CS%mode == 0) then
-          cff = sqrt(GV%g_Earth * G%bathyT(i,j+1) )
+          cff = sqrt(GV%g_Earth * depth_tot(i,j+1) )
           val2 = mag_SSH * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff)
           segment%eta(I,j) = GV%Z_to_H*val2 * cos(omega * time_sec)
-          segment%normal_vel_bt(I,j) = (val1 * cff * sina / &
-                 (G%bathyT(i,j+1) )) * val2
+          segment%normal_vel_bt(I,j) = (val1 * cff * sina / depth_tot(i,j+1) ) * val2
           if (segment%nudged) then
             do k=1,nz
-              segment%nudged_normal_vel(I,j,k) = (val1 * cff * sina / &
-                     (G%bathyT(i,j+1) )) * val2
+              segment%nudged_normal_vel(I,j,k) = (val1 * cff * sina / depth_tot(i,j+1)) * val2
             enddo
           elseif (segment%specified) then
             do k=1,nz
-              segment%normal_vel(I,j,k) = (val1 * cff * sina / &
-                     (G%bathyT(i,j+1) )) * val2
+              segment%normal_vel(I,j,k) = (val1 * cff * sina / depth_tot(i,j+1) ) * val2
               segment%normal_trans(i,J,k) = segment%normal_vel(i,J,k) * h(i,j+1,k) * G%dxCv(i,J)
             enddo
           endif
@@ -347,11 +353,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
           y1 = km_to_L_scale * G%geoLatBu(I,J)
           x = (x1 - CS%coast_offset1) * cosa + y1 * sina
           y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
-          cff = sqrt(GV%g_Earth * G%bathyT(i,j+1) )
+          cff = sqrt(GV%g_Earth * depth_tot(i,j+1) )
           val2 = mag_SSH * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff)
           if (CS%mode == 0) then ; do k=1,nz
-            segment%tangential_vel(I,J,k) = ((val1 * val2 * cff * sina) / &
-                ( 0.5*((G%bathyT(i+1,j+1)) + G%bathyT(i,j+1))) )
+            segment%tangential_vel(I,J,k) = (val1 * val2 * cff * sina) / &
+                ( 0.5*(depth_tot(i+1,j+1) + depth_tot(i,j+1)) )
           enddo ; endif
         enddo ; enddo
       endif
diff --git a/src/user/Neverworld_initialization.F90 b/src/user/Neverworld_initialization.F90
index 93a43e4a3e..3f5b8c8ab2 100644
--- a/src/user/Neverworld_initialization.F90
+++ b/src/user/Neverworld_initialization.F90
@@ -239,12 +239,14 @@ end function circ_ridge
 !! by finding the depths of interfaces in a specified latitude-dependent
 !! temperature profile with an exponentially decaying thermocline on top of a
 !! linear stratification.
-subroutine Neverworld_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, P_ref)
+subroutine Neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file, eqn_of_state, P_ref)
   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, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< The thickness that is being
                                                               !! initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file           !< A structure indicating the open
                                                               !! file to parse for model
                                                               !! parameter values.
@@ -283,7 +285,7 @@ subroutine Neverworld_initialize_thickness(h, G, GV, US, param_file, eqn_of_stat
   enddo
 
   do j=js,je ; do i=is,ie
-    e_interface = -G%bathyT(i,j)
+    e_interface = -depth_tot(i,j)
     do k=nz,2,-1
       h(i,j,k) = GV%Z_to_H * (e0(k) - e_interface) ! Nominal thickness
       x=(G%geoLonT(i,j)-G%west_lon)/G%len_lon
diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90
index dfa9c19460..448c86b5fb 100644
--- a/src/user/Phillips_initialization.F90
+++ b/src/user/Phillips_initialization.F90
@@ -35,12 +35,14 @@ module Phillips_initialization
 contains
 
 !> Initialize the thickness field for the Phillips model test case.
-subroutine Phillips_initialize_thickness(h, G, GV, US, param_file, just_read_params)
+subroutine Phillips_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read_params)
   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(out) :: h          !< The thickness that is being initialized [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file !< A structure indicating the open file
                                                      !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -113,7 +115,7 @@ subroutine Phillips_initialize_thickness(h, G, GV, US, param_file, just_read_par
     ! thicknesses are set to insure that: 1. each layer is at least an Angstrom thick, and
     ! 2. the interfaces are where they should be based on the resting depths and interface
     !    height perturbations, as long at this doesn't interfere with 1.
-    eta1D(nz+1) = -G%bathyT(i,j)
+    eta1D(nz+1) = -depth_tot(i,j)
     do k=nz,1,-1
       eta1D(K) = eta_im(j,K)
       if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
diff --git a/src/user/RGC_initialization.F90 b/src/user/RGC_initialization.F90
index 4df728c22a..d051bccc6c 100644
--- a/src/user/RGC_initialization.F90
+++ b/src/user/RGC_initialization.F90
@@ -46,7 +46,7 @@ module RGC_initialization
 
 !> Sets up the the inverse restoration time, and the values towards which the interface heights,
 !! velocities and tracers should be restored within the sponges for the RGC test case.
-subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp)
+subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, use_ALE, CSp, ACSp)
   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
@@ -59,6 +59,8 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp)
                  target, intent(in) :: u    !< Array with the u velocity [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                  target, intent(in) :: v    !< Array with the v velocity [L T-1 ~> m s-1]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                         intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type), intent(in) :: PF   !< A structure indicating the
                                             !! open file to parse for model
                                             !! parameter values.
@@ -93,10 +95,11 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp)
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
   iscB = G%iscB ; iecB = G%iecB; jscB = G%jscB ; jecB = G%jecB
 
-  call get_param(PF,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness',units='m',default=1.e-3)
+  call get_param(PF, mod,"MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
+                 units='m', default=1.e-3)
 
-  call get_param(PF, mod, "RGC_TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', &
-                 default=0.0, scale=86400.0*US%s_to_T)
+  call get_param(PF, mod, "RGC_TNUDG", TNUDG, 'Nudging time scale for sponge layers', &
+                 units='days', default=0.0, scale=86400.0*US%s_to_T)
 
   call get_param(PF, mod, "LENLAT", lenlat, &
                   "The latitudinal or y-direction length of the domain", &
@@ -114,7 +117,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp)
                  "Nudge velocities (u and v) towards zero in the sponge layer.", &
                  default=.false., do_not_log=.true.)
 
-  T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0; RHO(:,:,:) = 0.0
+  T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 ; RHO(:,:,:) = 0.0
 
   call get_param(PF, mod, "MINIMUM_DEPTH", min_depth, &
                  "The minimum depth of the ocean.", units="m", default=0.0)
@@ -130,7 +133,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp)
   !  and mask2dT is 1.
 
   do i=is,ie ; do j=js,je
-    if ((G%bathyT(i,j) <= min_depth) .or. (G%geoLonT(i,j) <= lensponge)) then
+    if ((depth_tot(i,j) <= min_depth) .or. (G%geoLonT(i,j) <= lensponge)) then
       Idamp(i,j) = 0.0
     elseif (G%geoLonT(i,j) >= (lenlon - lensponge) .AND. G%geoLonT(i,j) <= lenlon) then
       dummy1 = (G%geoLonT(i,j)-(lenlon - lensponge))/(lensponge)
diff --git a/src/user/adjustment_initialization.F90 b/src/user/adjustment_initialization.F90
index ad4eab33ff..b9f676dc55 100644
--- a/src/user/adjustment_initialization.F90
+++ b/src/user/adjustment_initialization.F90
@@ -57,8 +57,8 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read
   real    :: target_values(SZK_(GV)+1)  ! Target densities or density anomalies [R ~> kg m-3]
   logical :: just_read    ! If true, just read parameters but set nothing.
   character(len=20) :: verticalCoordinate
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   integer :: i, j, k, is, ie, js, je, nz
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
@@ -193,13 +193,15 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read
 end subroutine adjustment_initialize_thickness
 
 !> Initialization of temperature and salinity in the adjustment test case
-subroutine adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file, &
+subroutine adjustment_initialize_temperature_salinity(T, S, h, depth_tot, G, GV, param_file, &
                                                       eqn_of_state, just_read_params)
   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(out) :: T !< The temperature that is being initialized.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is being initialized.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h !< The model thicknesses [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file   !< A structure indicating the open file to
                                                       !! parse for model parameter values.
   type(EOS_type),                 pointer     :: eqn_of_state !< Equation of state.
@@ -260,7 +262,7 @@ subroutine adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file
     case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA )
       dSdz = -delta_S_strat / G%max_depth
       do j=js,je ; do i=is,ie
-        eta1d(nz+1) = -G%bathyT(i,j)
+        eta1d(nz+1) = -depth_tot(i,j)
         do k=nz,1,-1
           eta1d(k) = eta1d(k+1) + h(i,j,k)*GV%H_to_Z
         enddo
diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90
index b1c988e016..22f4d705a1 100644
--- a/src/user/baroclinic_zone_initialization.F90
+++ b/src/user/baroclinic_zone_initialization.F90
@@ -75,7 +75,7 @@ subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref,
 end subroutine bcz_params
 
 !> Initialization of temperature and salinity with the baroclinic zone initial conditions
-subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_file, &
+subroutine baroclinic_zone_init_temperature_salinity(T, S, h, depth_tot, G, GV, US, param_file, &
                                                      just_read_params)
   type(ocean_grid_type),                     intent(in)  :: G  !< Grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< The ocean's vertical grid structure.
@@ -83,6 +83,8 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_f
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T  !< Potential temperature [degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< Salinity [ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< The model thicknesses [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),                     intent(in)  :: param_file  !< Parameter file handle
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                       !! only read parameters without changing T & S.
@@ -109,7 +111,7 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_f
   PI = 4.*atan(1.)
 
   do j = G%jsc,G%jec ; do i = G%isc,G%iec
-    zi = -G%bathyT(i,j)
+    zi = -depth_tot(i,j)
     x = G%geoLonT(i,j) - (G%west_lon + 0.5*G%len_lon) ! Relative to center of domain
     xd = x / G%len_lon ! -1/2 < xd 1/2
     y = G%geoLatT(i,j) - (G%south_lat + 0.5*G%len_lat) ! Relative to center of domain
diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90
index ed0bbbf069..d077e0fa6f 100644
--- a/src/user/benchmark_initialization.F90
+++ b/src/user/benchmark_initialization.F90
@@ -82,13 +82,15 @@ end subroutine benchmark_initialize_topography
 !! by finding the depths of interfaces in a specified latitude-dependent
 !! temperature profile with an exponentially decaying thermocline on top of a
 !! linear stratification.
-subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, &
+subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, eqn_of_state, &
                                           P_Ref, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   type(EOS_type),          pointer     :: eqn_of_state !< Equation of state structure
@@ -181,7 +183,7 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state
     ! are set to insure that: 1. each layer is at least  Gv%Angstrom_m thick, and
     ! 2. the interfaces are where they should be based on the resting depths and interface
     ! height perturbations, as long at this doesn't interfere with 1.
-    eta1D(nz+1) = -G%bathyT(i,j)
+    eta1D(nz+1) = -depth_tot(i,j)
 
     do k=nz,2,-1
       T_int = 0.5*(T0(k) + T0(k-1))
diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90
index 4dd5a7c606..29fb6647b3 100644
--- a/src/user/circle_obcs_initialization.F90
+++ b/src/user/circle_obcs_initialization.F90
@@ -28,11 +28,13 @@ module circle_obcs_initialization
 contains
 
 !> This subroutine initializes layer thicknesses for the circle_obcs experiment.
-subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_params)
+subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot   !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -79,8 +81,8 @@ subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_para
   enddo
 
   ! Uniform thicknesses for base state
-  do j=js,je ; do i=is,ie                        !
-    eta1D(nz+1) = -G%bathyT(i,j)
+  do j=js,je ; do i=is,ie
+    eta1D(nz+1) = -depth_tot(i,j)
     do k=nz,1,-1
       eta1D(K) = e0(K)
       if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90
index e8fe345bb0..c1eb4fa2e7 100644
--- a/src/user/dense_water_initialization.F90
+++ b/src/user/dense_water_initialization.F90
@@ -152,11 +152,13 @@ subroutine dense_water_initialize_TS(G, GV, param_file, eqn_of_state, T, S, h, j
 end subroutine dense_water_initialize_TS
 
 !> Initialize the restoring sponges for the dense water experiment
-subroutine dense_water_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp)
+subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
   type(ocean_grid_type),   intent(in) :: G !< Horizontal grid control structure
   type(verticalGrid_type), intent(in) :: GV !< Vertical grid control structure
   type(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
   type(thermo_var_ptrs),   intent(in) :: tv !< Thermodynamic variables
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file !< Parameter file structure
   logical,                 intent(in) :: use_ALE !< ALE flag
   type(sponge_CS),         pointer    :: CSp !< Layered sponge control structure pointer
@@ -234,7 +236,7 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CS
 
     do j = G%jsc,G%jec
       do i = G%isc,G%iec
-        eta1D(nz+1) = -G%bathyT(i,j)
+        eta1D(nz+1) = -depth_tot(i,j)
         do k = nz,1,-1
           eta1D(k) = e0(k)
 
diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90
index f7b647dd27..463fe018b0 100644
--- a/src/user/dumbbell_initialization.F90
+++ b/src/user/dumbbell_initialization.F90
@@ -90,12 +90,14 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth )
 end subroutine dumbbell_initialize_topography
 
 !> Initializes the layer thicknesses to be uniform in the dumbbell test case
-subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
+subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -169,7 +171,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p
       e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom
     enddo
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) = e0(k)
         if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
@@ -184,7 +186,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p
   case ( REGRIDDING_ZSTAR )                       ! Initial thicknesses for z coordinates
     if (just_read) return ! All run-time parameters have been read, so return.
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) = -G%max_depth * real(k-1) / real(nz)
         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
@@ -199,7 +201,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p
   case ( REGRIDDING_SIGMA )             ! Initial thicknesses for sigma coordinates
     if (just_read) return ! All run-time parameters have been read, so return.
     do j=js,je ; do i=is,ie
-      h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz)
+      h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz)
     enddo ; enddo
 
 end select
@@ -284,11 +286,13 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file
 end subroutine dumbbell_initialize_temperature_salinity
 
 !> Initialize the restoring sponges for the dumbbell test case
-subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp)
+subroutine dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
   type(ocean_grid_type),   intent(in) :: G !< Horizontal grid control structure
   type(verticalGrid_type), intent(in) :: GV !< Vertical grid control structure
   type(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
   type(thermo_var_ptrs),   intent(in) :: tv !< Thermodynamic variables
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in) :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in) :: param_file !< Parameter file structure
   logical,                 intent(in) :: use_ALE !< ALE flag
   type(sponge_CS),         pointer    :: CSp !< Layered sponge control structure pointer
@@ -354,7 +358,7 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp,
   if (use_ALE) then
     ! construct a uniform grid for the sponge
     do j=G%jsc,G%jec ; do i=G%isc,G%iec
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) =  depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) = -G%max_depth * real(k-1) / real(nz)
         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90
index 9118133108..6bfaedc221 100644
--- a/src/user/seamount_initialization.F90
+++ b/src/user/seamount_initialization.F90
@@ -77,12 +77,14 @@ end subroutine seamount_initialize_topography
 
 !> Initialization of thicknesses.
 !! This subroutine initializes the layer thicknesses to be uniform.
-subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
+subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -152,7 +154,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p
       e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom
     enddo
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) = e0(k)
         if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
@@ -167,7 +169,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p
   case ( REGRIDDING_ZSTAR )                       ! Initial thicknesses for z coordinates
     if (just_read) return ! All run-time parameters have been read, so return.
     do j=js,je ; do i=is,ie
-      eta1D(nz+1) = -G%bathyT(i,j)
+      eta1D(nz+1) = -depth_tot(i,j)
       do k=nz,1,-1
         eta1D(k) =  -G%max_depth * real(k-1) / real(nz)
         if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
@@ -182,7 +184,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p
   case ( REGRIDDING_SIGMA )             ! Initial thicknesses for sigma coordinates
     if (just_read) return ! All run-time parameters have been read, so return.
     do j=js,je ; do i=is,ie
-      h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz)
+      h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz)
     enddo ; enddo
 
 end select
diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90
index e1c0a96d63..1c3334d8b0 100644
--- a/src/user/sloshing_initialization.F90
+++ b/src/user/sloshing_initialization.F90
@@ -53,12 +53,14 @@ end subroutine sloshing_initialize_topography
 !! same thickness but all interfaces (except bottom and sea surface) are
 !! displaced according to a half-period cosine, with maximum value on the
 !! left and minimum value on the right. This sets off a regular sloshing motion.
-subroutine sloshing_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
+subroutine sloshing_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params)
   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(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot  !< The nominal total depth of the ocean [Z ~> m]
   type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                       !! to parse for model parameter values.
   logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
@@ -152,7 +154,7 @@ subroutine sloshing_initialize_thickness ( h, G, GV, US, param_file, just_read_p
     enddo
 
     ! 3. The last interface must coincide with the seabed
-    z_inter(nz+1) = -G%bathyT(i,j)
+    z_inter(nz+1) = -depth_tot(i,j)
     ! Modify interface heights to make sure all thicknesses are strictly positive
     do k = nz,1,-1
       if ( z_inter(k) < (z_inter(k+1) + GV%Angstrom_Z) ) then
diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90
index ac6ec8c4bc..f62aa54f88 100644
--- a/src/user/soliton_initialization.F90
+++ b/src/user/soliton_initialization.F90
@@ -28,12 +28,14 @@ module soliton_initialization
 contains
 
 !> Initialization of thicknesses in Equatorial Rossby soliton test
-subroutine soliton_initialize_thickness(h, G, GV, US)
+subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US)
   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(out) :: h    !< The thickness that is being initialized [H ~> m or kg m-2].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)  :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
 
   integer :: i, j, k, is, ie, js, je, nz
   real    :: x, y, x0, y0
@@ -55,7 +57,7 @@ subroutine soliton_initialize_thickness(h, G, GV, US)
       y = G%geoLatT(i,j)-y0
       val3 = exp(-val1*x)
       val4 = val2 * ( 2.0*val3 / (1.0 + (val3*val3)) )**2
-      h(i,j,k) = GV%Z_to_H * (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + G%bathyT(i,j))
+      h(i,j,k) = GV%Z_to_H * (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j))
     enddo
   enddo ; enddo
 
diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90
index 793b87f149..3338121d9e 100644
--- a/src/user/user_initialization.F90
+++ b/src/user/user_initialization.F90
@@ -247,7 +247,7 @@ end subroutine write_user_log
 !!  - u - Zonal velocity [Z T-1 ~> m s-1].
 !!  - v - Meridional velocity [Z T-1 ~> m s-1].
 !!  - h - Layer thickness [H ~> m or kg m-2].  (Must be positive.)
-!!  - G%bathyT - Basin depth [Z ~> m].  (Must be positive.)
+!!  - G%bathyT - Basin depth [Z ~> m].
 !!  - G%CoriolisBu - The Coriolis parameter [T-1 ~> s-1].
 !!  - GV%g_prime - The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2].
 !!  - GV%Rlay - Layer potential density (coordinate variable) [R ~> kg m-3].