diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index c6017c7465..d0a7da8695 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -958,7 +958,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
             call check_redundant("Pre-ALE 1 ", u, v, G)
           endif
           call cpu_clock_begin(id_clock_ALE)
+          ! Switch thickness units from H to m for remapping
+          h = h*G%H_to_m
           call ALE_main(G, h, u, v, CS%tv, CS%ALE_CSp)
+          h = h*G%m_to_H
           call cpu_clock_end(id_clock_ALE)
           if (CS%debug) then
             call MOM_state_chksum("Post-ALE 1 ", u, v, h, CS%uh, CS%vh, G)
@@ -1257,7 +1260,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
             call check_redundant("Pre-ALE ", u, v, G)
           endif
           call cpu_clock_begin(id_clock_ALE)
+          ! Switch thickness units from H to m for remapping
+          h = h*G%H_to_m
           call ALE_main(G, h, u, v, CS%tv, CS%ALE_CSp)
+          h = h*G%m_to_H
           call cpu_clock_end(id_clock_ALE)
           if (CS%debug) then
             call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G)
@@ -1854,7 +1860,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   ! Allocate and initialize space for primary MOM variables.
   ALLOC_(CS%u(IsdB:IedB,jsd:jed,nz))   ; CS%u(:,:,:) = 0.0
   ALLOC_(CS%v(isd:ied,JsdB:JedB,nz))   ; CS%v(:,:,:) = 0.0
-  ALLOC_(CS%h(isd:ied,jsd:jed,nz))     ; CS%h(:,:,:) = G%Angstrom
+  ALLOC_(CS%h(isd:ied,jsd:jed,nz))     ; CS%h(:,:,:) = G%Angstrom_z
   ALLOC_(CS%uh(IsdB:IedB,jsd:jed,nz))  ; CS%uh(:,:,:) = 0.0
   ALLOC_(CS%vh(isd:ied,JsdB:JedB,nz))  ; CS%vh(:,:,:) = 0.0
   if (CS%use_temperature) then
@@ -1988,18 +1994,25 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
     if (CS%debug) then
       call uchksum(CS%u,"Pre initialize_ALE u", G, haloshift=1)
       call vchksum(CS%v,"Pre initialize_ALE v", G, haloshift=1)
-      call hchksum(CS%h,"Pre initialize_ALE h", G, haloshift=1)
+      call hchksum(CS%h*G%H_to_m,"Pre initialize_ALE h", G, haloshift=1)
     endif
+    ! Switch thickness units from H to m for regridding
+    CS%h = CS%h*G%H_to_m
+
     if (.not. query_initialized(CS%h,"h",CS%restart_CSp)) then
       ! This is a not a restart so we do the following...
       call adjustGridForIntegrity(CS%ALE_CSp, G, CS%h )
       call ALE_main( G, CS%h, CS%u, CS%v, CS%tv, CS%ALE_CSp )
     endif
+
+    ! Switch thickness units back to H
+    CS%h = CS%h*G%m_to_H
+
     call ALE_updateVerticalGridType( CS%ALE_CSp, G%GV )
     if (CS%debug) then
       call uchksum(CS%u,"Post initialize_ALE u", G, haloshift=1)
       call vchksum(CS%v,"Post initialize_ALE v", G, haloshift=1)
-      call hchksum(CS%h, "Post initialize_ALE h", G, haloshift=1)
+      call hchksum(CS%h*G%H_to_m, "Post initialize_ALE h", G, haloshift=1)
     endif
   endif
 
@@ -2652,8 +2665,8 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
       enddo
 
       do k=1,nz ; do i=is,ie
-        if (depth(i) + h(i,j,k) < depth_ml) then
-          dh = h(i,j,k)
+        if (depth(i) + h(i,j,k)*G%H_to_m < depth_ml) then
+          dh = h(i,j,k)*G%H_to_m
         elseif (depth(i) < depth_ml) then
           dh = depth_ml - depth(i)
         else
@@ -2669,7 +2682,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
       enddo ; enddo
   ! Calculate the average properties of the mixed layer depth.
       do i=is,ie
-        if (depth(i) < G%H_subroundoff) depth(i) = G%H_subroundoff
+        if (depth(i) < G%H_subroundoff*G%H_to_m) depth(i) = G%H_subroundoff*G%H_to_m
         if (CS%use_temperature) then
           state%SST(i,j) = state%SST(i,j) / depth(i)
           state%SSS(i,j) = state%SSS(i,j) / depth(i)
diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90
index bb862464f9..c7f1444f7b 100644
--- a/src/core/MOM_PressureForce_Montgomery.F90
+++ b/src/core/MOM_PressureForce_Montgomery.F90
@@ -516,7 +516,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
     do j=Jsq,Jeq+1 
       do i=Isq,Ieq+1 ; e(i,j,1) = -1.0*G%bathyT(i,j) ; enddo
       do k=1,nz ; do i=Isq,Ieq+1
-        e(i,j,1) = e(i,j,1) + h(i,j,k)
+        e(i,j,1) = e(i,j,1) + h(i,j,k)*G%H_to_m
       enddo ; enddo 
     enddo
     call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp)
@@ -537,7 +537,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
   endif
 !$OMP do
   do j=Jsq,Jeq+1 ; do k=nz,1,-1 ; do i=Isq,Ieq+1
-    e(i,j,K) = e(i,j,K+1) + h(i,j,k)
+    e(i,j,K) = e(i,j,K+1) + h(i,j,k)*G%H_to_m
   enddo ; enddo ; enddo
 !$OMP end parallel
   if (use_EOS) then
@@ -662,12 +662,12 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, CS, p_atm, pbce, eta)
     ! about 200 lines above.
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e,e_tidal)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-        eta(i,j) = e(i,j,1) + e_tidal(i,j)
+        eta(i,j) = e(i,j,1)*G%m_to_H + e_tidal(i,j)*G%m_to_H
       enddo ; enddo
     else
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-        eta(i,j) = e(i,j,1)
+        eta(i,j) = e(i,j,1)*G%m_to_H
       enddo ; enddo
     endif
   endif
@@ -729,7 +729,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
   Rho0xG = Rho0*g_Earth
   G_Rho0 = g_Earth/Rho0
   use_EOS = associated(tv%eqn_of_state)
-  h_neglect = G%H_subroundoff * G%H_to_m
+  h_neglect = G%H_subroundoff*G%H_to_m
 
   if (use_EOS) then
     if (present(rho_star)) then
@@ -737,8 +737,8 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
 !$OMP                          private(Ihtot)
       do j=Jsq,Jeq+1
         do i=Isq,Ieq+1
-          Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
-          pbce(i,j,1) = GFS_scale * rho_star(i,j,1)
+          Ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * G%m_to_H)
+          pbce(i,j,1) = GFS_scale * rho_star(i,j,1) * (1.0 / G%m_to_H)
         enddo
         do k=2,nz ; do i=Isq,Ieq+1
           pbce(i,j,k) = pbce(i,j,k-1) + (rho_star(i,j,k)-rho_star(i,j,k-1)) * &
@@ -750,13 +750,13 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
 !$OMP                          private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS)
       do j=Jsq,Jeq+1
         do i=Isq,Ieq+1
-          Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
+          Ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * G%m_to_H)
           press(i) = -Rho0xG*e(i,j,1)
         enddo
         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
                                Isq, Ieq-Isq+2, tv%eqn_of_state)
         do i=Isq,Ieq+1
-          pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i))
+          pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * (1.0 / G%m_to_H)
         enddo
         do k=2,nz
           do i=Isq,Ieq+1
@@ -779,8 +779,8 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,h_neglect,pbce) private(Ihtot)
     do j=Jsq,Jeq+1
       do i=Isq,Ieq+1
-        Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
-        pbce(i,j,1) = G%g_prime(1)
+        Ihtot(i) = 1.0 / (((e(i,j,1)-e(i,j,nz+1)) + h_neglect) * G%m_to_H)
+        pbce(i,j,1) = G%g_prime(1) * (1.0 / G%m_to_H)
       enddo
       do k=2,nz ; do i=Isq,Ieq+1
         pbce(i,j,k) = pbce(i,j,k-1) + &
diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90
index eff2c2c09f..026ed7c415 100644
--- a/src/core/MOM_PressureForce_analytic_FV.F90
+++ b/src/core/MOM_PressureForce_analytic_FV.F90
@@ -576,13 +576,13 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
     ! but that is not yet implemented, and the current form is correct for
     ! barotropic tides.
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,e,G,h)
-    do j=Jsq,Jeq+1 
+    do j=Jsq,Jeq+1
       do i=Isq,Ieq+1
         e(i,j,1) = -1.0*G%bathyT(i,j)
       enddo
       do k=1,nz ; do i=Isq,Ieq+1
-        e(i,j,1) = e(i,j,1) + G%H_to_m*h(i,j,k)
-      enddo ; enddo 
+        e(i,j,1) = e(i,j,1) + h(i,j,k)*G%H_to_m
+      enddo ; enddo
     enddo
     call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp)
   endif
@@ -602,7 +602,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
   endif
 !$OMP do 
   do j=Jsq,Jeq+1; do k=nz,1,-1 ; do i=Isq,Ieq+1
-    e(i,j,K) = e(i,j,K+1) + G%H_to_m*h(i,j,k)
+    e(i,j,K) = e(i,j,K+1) + h(i,j,k)*G%H_to_m
   enddo ; enddo ; enddo
 !$OMP end parallel
 
@@ -677,9 +677,9 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
   ! of freedeom needed to know the linear profile).
   if ( use_ALE ) then
     if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then
-      call pressure_gradient_plm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h );
+      call pressure_gradient_plm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h*G%H_to_m );
     elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then
-      call pressure_gradient_ppm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h );
+      call pressure_gradient_ppm (ALE_CSp, S_t, S_b, T_t, T_b, G, tv, h*G%H_to_m );
     endif
   endif
 
@@ -754,6 +754,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
                     rho_ref, CS%Rho0, G%g_Earth, G%Block(n), tv%eqn_of_state, &
                     dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk )
         endif
+        intz_dpa_bk = intz_dpa_bk*G%m_to_H
       else
         do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1
           i = ib+ioff_bk ; j = jb+joff_bk
@@ -775,7 +776,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
         PFu(I,j,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
                      (pa_bk(ib+1,jb)*h(i+1,j,k) + intz_dpa_bk(ib+1,jb))) + &
                      ((h(i+1,j,k) - h(i,j,k)) * intx_pa_bk(Ib,jb) - &
-                     (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa_bk(Ib,jb))) * &
+                     (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa_bk(Ib,jb) * G%m_to_H)) * &
                      ((2.0*I_Rho0*G%IdxCu(I,j)) / &
                      ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
         intx_pa_bk(Ib,jb) = intx_pa_bk(Ib,jb) + intx_dpa_bk(Ib,jb)
@@ -786,7 +787,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
         PFv(i,J,k) = (((pa_bk(ib,jb)*h(i,j,k) + intz_dpa_bk(ib,jb)) - &
                      (pa_bk(ib,jb+1)*h(i,j+1,k) + intz_dpa_bk(ib,jb+1))) + &
                      ((h(i,j+1,k) - h(i,j,k)) * inty_pa_bk(ib,Jb) - &
-                     (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa_bk(ib,Jb))) * &
+                     (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa_bk(ib,Jb) * G%m_to_H)) * &
                      ((2.0*I_Rho0*G%IdyCv(i,J)) / &
                      ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
         inty_pa_bk(ib,Jb) = inty_pa_bk(ib,Jb) + inty_dpa_bk(ib,Jb)
@@ -819,12 +820,12 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, CS, ALE_CSp, p_atm, pbce,
     ! about 200 lines above.
 !$OM parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e,e_tidal)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-        eta(i,j) = e(i,j,1) + e_tidal(i,j)
+        eta(i,j) = e(i,j,1)*G%m_to_H + e_tidal(i,j)*G%m_to_H
       enddo ; enddo
     else
 !$OM parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta,e)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-        eta(i,j) = e(i,j,1)
+        eta(i,j) = e(i,j,1)*G%m_to_H
       enddo ; enddo
     endif
   endif
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 08a7d62cb6..146c90f298 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1375,7 +1375,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, &
         ! Limit the sink (inward) correction to the amount of mass that is already
         ! inside the cell, plus any mass added by eta_source.
         Htot = eta(i,j)
-        if (G%Boussinesq) Htot = CS%bathyT(i,j) + eta(i,j)
+        if (G%Boussinesq) Htot = CS%bathyT(i,j)*G%m_to_H + eta(i,j)
 
         CS%eta_cor(i,j) = max(CS%eta_cor(i,j), -max(0.0,Htot + dt*CS%eta_source(i,j)))
       endif
@@ -2657,8 +2657,8 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, MS, halo, use_BT_cont, Datu, Datv,
         BT_OBC%Cg_u(I,j) = SQRT(G%g_prime(1)*(0.5* &
                                 (G%bathyT(i,j) + G%bathyT(i+1,j))))
         if (G%Boussinesq) then
-          BT_OBC%H_u(I,j) = 0.5*((G%bathyT(i,j) + eta(i,j)) + &
-                                 (G%bathyT(i+1,j) + eta(i+1,j)))
+          BT_OBC%H_u(I,j) = 0.5*((G%bathyT(i,j)*G%m_to_H + eta(i,j)) + &
+                                 (G%bathyT(i+1,j)*G%m_to_H + eta(i+1,j)))
         else
           BT_OBC%H_u(I,j) = 0.5*(eta(i,j) + eta(i+1,j))
         endif
@@ -2693,8 +2693,8 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, MS, halo, use_BT_cont, Datu, Datv,
         BT_OBC%Cg_v(i,J) = SQRT(G%g_prime(1)*(0.5* &
                                 (G%bathyT(i,j) + G%bathyT(i,j+1))))
         if (G%Boussinesq) then
-          BT_OBC%H_v(i,J) = 0.5*((G%bathyT(i,j) + eta(i,j)) + &
-                                 (G%bathyT(i,j+1) + eta(i,j+1)))
+          BT_OBC%H_v(i,J) = 0.5*((G%bathyT(i,j)*G%m_to_H + eta(i,j)) + &
+                                 (G%bathyT(i,j+1)*G%m_to_H + eta(i,j+1)))
         else
           BT_OBC%H_v(i,J) = 0.5*(eta(i,j) + eta(i,j+1))
         endif
@@ -3487,14 +3487,14 @@ subroutine find_face_areas(Datu, Datv, G, CS, MS, eta, halo, add_max)
     if (G%Boussinesq) then
 !$OMP do
       do j=js-hs,je+hs ; do I=is-1-hs,ie+hs
-        H1 = CS%bathyT(i,j) + eta(i,j) ; H2 = CS%bathyT(i+1,j) + eta(i+1,j)
+        H1 = CS%bathyT(i,j)*G%m_to_H + eta(i,j) ; H2 = CS%bathyT(i+1,j)*G%m_to_H + eta(i+1,j)
         Datu(I,j) = 0.0 ; if ((H1 > 0.0) .and. (H2 > 0.0)) &
         Datu(I,j) = CS%dy_Cu(I,j) * (2.0 * H1 * H2) / (H1 + H2)
 !       Datu(I,j) = CS%dy_Cu(I,j) * 0.5 * (H1 + H2)
       enddo ; enddo
 !$OMP do
       do J=js-1-hs,je+hs ; do i=is-hs,ie+hs
-        H1 = CS%bathyT(i,j) + eta(i,j) ; H2 = CS%bathyT(i,j+1) + eta(i,j+1)
+        H1 = CS%bathyT(i,j)*G%m_to_H + eta(i,j) ; H2 = CS%bathyT(i,j+1)*G%m_to_H + eta(i,j+1)
         Datv(i,J) = 0.0 ; if ((H1 > 0.0) .and. (H2 > 0.0)) &
         Datv(i,J) = CS%dx_Cv(i,J) * (2.0 * H1 * H2) / (H1 + H2)
 !       Datv(i,J) = CS%dy_v(i,J) * 0.5 * (H1 + H2)
@@ -3597,7 +3597,7 @@ subroutine bt_mass_source(h, eta, fluxes, set_cor, dt_therm, dt_since_therm, &
   do j=js,je
     do i=is,ie ; h_tot(i) = h(i,j,1) ; enddo
     if (G%Boussinesq) then
-      do i=is,ie ; eta_h(i) = h(i,j,1) - G%bathyT(i,j) ; enddo
+      do i=is,ie ; eta_h(i) = h(i,j,1) - G%bathyT(i,j)*G%m_to_H ; enddo
     else
       do i=is,ie ; eta_h(i) = h(i,j,1) ; enddo
     endif
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 570424e17f..d5d35161ee 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -496,12 +496,13 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   call cpu_clock_begin(id_clock_pres)
   call PressureForce(h, tv, CS%PFu, CS%PFv, G, CS%PressureForce_CSp, &
                      CS%ALE_CSp, p_surf, CS%pbce, CS%eta_PF)
+  if (CS%debug) then
+    call uchksum(CS%PFu,"After PressureForce CS%PFu",G,haloshift=0)
+    call vchksum(CS%PFv,"After PressureForce CS%PFv",G,haloshift=0)
+  endif
+
   if (dyn_p_surf) then
-    if (G%Boussinesq) then
-      Pa_to_eta = 1.0 / (G%Rho0*G%g_Earth)
-    else
-      Pa_to_eta = 1.0 / G%H_to_Pa
-    endif
+    Pa_to_eta = 1.0 / G%H_to_Pa
 !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,eta_PF_start,CS,Pa_to_eta,p_surf_begin,p_surf_end)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       eta_PF_start(i,j) = CS%eta_PF(i,j) - Pa_to_eta * &
@@ -1242,7 +1243,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, param_file, &
     ! dimensions as h, either m or kg m-3.  
     !   CS%eta(:,:) = 0.0 already from initialization.
     if (G%Boussinesq) then
-      do j=js,je ; do i=is,ie ; CS%eta(i,j) = -G%bathyT(i,j) ; enddo ; enddo
+      do j=js,je ; do i=is,ie ; CS%eta(i,j) = -G%bathyT(i,j) * G%m_to_H ; enddo ; enddo
     endif
     do k=1,nz ; do j=js,je ; do i=is,ie
        CS%eta(i,j) = CS%eta(i,j) + h(i,j,k)
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index 022dda86af..6556d1dc6b 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -687,7 +687,7 @@ subroutine calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, &
   useRiverHeatContent   = .False.    
   useCalvingHeatContent = .False.  
 
-  depthBeforeScalingFluxes = max( G%Angstrom, 1.e-30 )
+  depthBeforeScalingFluxes = max( G%Angstrom_z, 1.e-30 )
   pressure(:) = 0. ! Ignore atmospheric pressure
   GoRho       = G%g_Earth / G%Rho0
   start       = 1 + G%isc - G%isd
@@ -724,10 +724,10 @@ subroutine calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, &
 
   ! Convert to a buoyancy flux, excluding penetrating SW heating
   buoyancyFlux(G%isc:G%iec,1) = - GoRho * ( dRhodS(G%isc:G%iec) * netSalt(G%isc:G%iec) + &
-                                             dRhodT(G%isc:G%iec) * netHeat(G%isc:G%iec) ) ! m^2/s^3
+                                             dRhodT(G%isc:G%iec) * netHeat(G%isc:G%iec) ) * G%H_to_m ! m^2/s^3
   ! We also have a penetrative buoyancy flux associated with penetrative SW
   do k=2, G%ke+1
-    buoyancyFlux(G%isc:G%iec,k) = - GoRho * ( dRhodT(G%isc:G%iec) * netPen(G%isc:G%iec,k) ) ! m^2/s^3
+    buoyancyFlux(G%isc:G%iec,k) = - GoRho * ( dRhodT(G%isc:G%iec) * netPen(G%isc:G%iec,k) ) * G%H_to_m ! m^2/s^3
   enddo
 
 end subroutine calculateBuoyancyFlux1d
diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90
index eee7f3732c..6027d2c7a1 100644
--- a/src/core/MOM_grid.F90
+++ b/src/core/MOM_grid.F90
@@ -236,11 +236,16 @@ subroutine MOM_grid_init(G, param_file)
   call get_param(param_file, "MOM_grid", "ANGSTROM", G%Angstrom_z, &
                  "The minumum layer thickness, usually one-Angstrom.", &
                  units="m", default=1.0e-10)
-  if (.not.G%Boussinesq) &
+
+  if (.not.G%Boussinesq) then
     call get_param(param_file, "MOM_grid", "H_TO_KG_M2", G%H_to_kg_m2,&
                  "A constant that translates thicknesses from the model's \n"//&
                  "internal units of thickness to kg m-2.", units="kg m-2 H-1", &
                  default=1.0)
+  else
+    call get_param(param_file, "MOM_grid", "H_TO_M", G%H_to_m, default=1.0)
+  endif
+
   call get_param(param_file, "MOM_grid", "BATHYMETRY_AT_VEL", G%bathymetry_at_vel, &
                  "If true, there are separate values for the basin depths \n"//&
                  "at velocity points.  Otherwise the effects of of \n"//&
@@ -306,11 +311,10 @@ subroutine MOM_grid_init(G, param_file)
   endif
 
   if (G%Boussinesq) then
-    G%H_to_kg_m2 = G%Rho0
-    G%kg_m2_to_H = 1.0/G%Rho0
-    G%m_to_H = 1.0
-    G%H_to_m = 1.0
-    G%Angstrom = G%Angstrom_z
+    G%H_to_kg_m2 = G%Rho0 * G%H_to_m
+    G%kg_m2_to_H = 1.0 / G%H_to_kg_m2
+    G%m_to_H = 1.0 / G%H_to_m
+    G%Angstrom = G%Angstrom_z*G%m_to_H
   else
     G%kg_m2_to_H = 1.0 / G%H_to_kg_m2
     G%m_to_H = G%Rho0 * G%kg_m2_to_H
diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90
index e117cd4c2a..4d333b664a 100644
--- a/src/core/MOM_interface_heights.F90
+++ b/src/core/MOM_interface_heights.F90
@@ -114,12 +114,12 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size)
 !$OMP                               G_Earth,dz_geo,halo,I_gEarth) &
 !$OMP                       private(dilate,htot)
 !$OMP do
-  do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -G%bathyT(i,j) ; enddo ; enddo
+  do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -G%bathyT(i,j); enddo ; enddo
 
   if (G%Boussinesq) then
 !$OMP do
     do j=jsv,jev ; do k=nz,1,-1; do i=isv,iev
-      eta(i,j,K) = eta(i,j,K+1) + h(i,j,k)
+      eta(i,j,K) = eta(i,j,K+1) + h(i,j,k)*G%H_to_m
     enddo ; enddo ; enddo
     if (present(eta_bt)) then
       ! Dilate the water column to agree with the free surface height
@@ -127,10 +127,10 @@ subroutine find_eta_3d(h, tv, G_Earth, G, eta, eta_bt, halo_size)
 !$OMP do
       do j=jsv,jev
         do i=isv,iev
-          dilate(i) = (eta_bt(i,j)+G%bathyT(i,j)) / (eta(i,j,1)+G%bathyT(i,j))
+          dilate(i) = (eta_bt(i,j)+G%bathyT(i,j)*G%m_to_H) / (eta(i,j,1)+G%bathyT(i,j))
         enddo
         do k=1,nz ; do i=isv,iev
-          eta(i,j,K) = dilate(i) * (eta(i,j,K) + G%bathyT(i,j)) - G%bathyT(i,j)
+          eta(i,j,K) = dilate(i) * (eta(i,j,K) + G%bathyT(i,j)*G%m_to_H) - G%bathyT(i,j)
         enddo ; enddo
       enddo
     endif
@@ -234,7 +234,7 @@ subroutine find_eta_2d(h, tv, G_Earth, G, eta, eta_bt, halo_size)
     else
 !$OMP do
       do j=js,je ; do k=1,nz ; do i=is,ie
-        eta(i,j) = eta(i,j) + h(i,j,k)
+        eta(i,j) = eta(i,j) + h(i,j,k)*G%H_to_m
       enddo ; enddo ; enddo
     endif
   else
diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90
index 5cda60ab51..18bfb93398 100644
--- a/src/diagnostics/MOM_sum_output.F90
+++ b/src/diagnostics/MOM_sum_output.F90
@@ -638,7 +638,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, CS, tracer_CSp)
       do j=js,je ; do i=is,ie
         hbelow = 0.0
         do k=nz,1,-1
-          hbelow = hbelow + h(i,j,k)
+          hbelow = hbelow + h(i,j,k) * G%H_to_m
           hint = (H_0APE(K) + hbelow - G%bathyT(i,j))
           hbot = H_0APE(K) - G%bathyT(i,j)
           hbot = (hbot + ABS(hbot)) * 0.5
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 450da3ebdf..6174f2647f 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -325,7 +325,19 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, PF, dirs, &
                  "units of m to kg m-2 or vice versa, depending on whether \n"//&
                  "BOUSSINESQ is defined. This does not apply if a restart \n"//&
                  "file is read.", default=.false.)
-    if (convert) call convert_thickness(h, G, PF, tv)
+    if (convert .and. .not. G%Boussinesq) then
+      ! Convert h from m to kg m-2 then to thickness units (H)
+      call convert_thickness(h, G, PF, tv)
+    elseif (G%Boussinesq) then
+      ! Convert h from m to thickness units (H)
+      h = h*G%m_to_H
+    endif
+
+    if (debug) then
+      call hchksum(h*G%H_to_m, "MOM_initialize_state: h ", G, haloshift=1)
+      if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G, haloshift=1)
+      if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G, haloshift=1)
+    endif
 
 !  Remove the mass that would be displaced by an ice shelf or inverse barometer.
     call get_param(PF, mod, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90
index be40cf8151..a82dc56717 100644
--- a/src/parameterizations/vertical/MOM_KPP.F90
+++ b/src/parameterizations/vertical/MOM_KPP.F90
@@ -570,7 +570,7 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, &
         km1 = max(1, k-1)
         kk = 3*(k-1)
         deltaRho(k) = rho_1D(kk+2) - rho_1D(kk+1)
-        N2_1d(k)    = (GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / (0.5*(h(i,j,km1) + h(i,j,k))+G%H_subroundoff)
+        N2_1d(k)    = (GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / ((0.5*(h(i,j,km1) + h(i,j,k))+G%H_subroundoff)*G%H_to_m)
         N_1d(k)     = sqrt( max( N2_1d(k), 0.) )
       enddo
       N2_1d(G%ke+1 ) = 0.0
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 9a81cf5131..1bfeb343c0 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -2412,7 +2412,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, &
   ! To accommodate vanishing upper layers, we need to allow for an instantaneous
   ! distribution of forcing over some finite vertical extent. The bulk mixed layer
   ! code handles this issue properly. 
-  H_limit_fluxes = max(G%Angstrom, 1.E-30) 
+  H_limit_fluxes = max(G%Angstrom_z, 1.E-30) 
 
   ! The inverse scale, IforcingDepthScale, is a hack which 
   ! should not be tickled in Eulerian mode. It stops all of the forcing from 
@@ -2445,7 +2445,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, &
       h2d(i,k) = h(i,j,k)
       T2d(i,k) = tv%T(i,j,k)
       do n=1,nsw
-        opacityBand(n,i,k) = G%H_to_m*optics%opacity_band(n,i,j,k)
+        opacityBand(n,i,k) = (1.0 / G%m_to_H)*optics%opacity_band(n,i,j,k)
       enddo
     enddo ; enddo
 
@@ -2589,7 +2589,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, dt, fluxes, optics, ea, h, tv, &
           ! For layers thin relative to 1/IforcingDepthScale, then distribute 
           ! forcing into deeper layers. 
           ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
-          fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale)
+          fractionOfForcing = min(1.0, h2d(i,k)*G%H_to_m*IforcingDepthScale)
 
           ! In the case with (-1)*netMassOut greater than 0.8*h, then we limit 
           ! applied to the top cell, and distribute the fluxes downwards.
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index 44ea4cf499..c728f9f217 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -1654,7 +1654,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, G, CS,
       ! This is energy loss in addition to work done as mixing, apparently to Joule heating.
       TKE_remaining = exp(-Idecay*dh) * TKE_remaining
 
-      z = z + h(i,j,k) ! Distance between upper interface of layer and the bottom, in m.
+      z = z + h(i,j,k)*G%H_to_m ! Distance between upper interface of layer and the bottom, in m.
       D_minus_z = max(total_thickness - z, 0.) ! Thickness above layer, m.
 
       ! Diffusivity using law of the wall, limited by rotation, at height z, in m2/s.
diff --git a/src/parameterizations/vertical/MOM_shortwave_abs.F90 b/src/parameterizations/vertical/MOM_shortwave_abs.F90
index 80af047ff0..327d515b7f 100644
--- a/src/parameterizations/vertical/MOM_shortwave_abs.F90
+++ b/src/parameterizations/vertical/MOM_shortwave_abs.F90
@@ -385,7 +385,7 @@ subroutine sumSWoverBands(G, h, opacity_band, nsw, j, dt, &
       if (h(i,k) > 0.0) then
         do n=1,nsw ; if (Pen_SW_bnd(n,i) > 0.0) then
           ! SW_trans is the SW that is transmitted THROUGH the layer
-          opt_depth = h(i,k) * opacity_band(n,i,k)
+          opt_depth = h(i,k)*G%H_to_m * opacity_band(n,i,k)
           exp_OD = exp(-opt_depth)
           SW_trans = exp_OD
           ! Heating at a rate of less than 10-4 W m-2 = 10-3 K m / Century,